{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fe40a09b",
   "metadata": {},
   "source": [
    "# Read and plot Sigmet files available on AWS using Xradar \n",
    "This example shows how to access radar data from the Colombian national radar network public on Amazon Web Services. We will look at the bucket structure and plot a PPI using the Xradar library. Radar reflectivity is filtered using some polarimetric values and xarray functionality."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "924ac38f",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9d2cc01",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fsspec\n",
    "import xarray as xr\n",
    "import xradar as xd\n",
    "import boto3\n",
    "import botocore\n",
    "import numpy as np\n",
    "from pandas import to_datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import cmweather\n",
    "from botocore.client import Config\n",
    "from datetime import datetime\n",
    "from matplotlib import pyplot\n",
    "from pyart.graph import cm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59dbfdb5",
   "metadata": {},
   "source": [
    "## IDEAM AWS Bucket\n",
    "Instituto de Hidrología, Meteorología y Estudios Ambientales - IDEAM (Colombian National Weather Service) has made public the weather radar data. Data can be found [here](https://registry.opendata.aws/ideam-radares/), and documentation [here](http://www.pronosticosyalertas.gov.co/archivos-radar#:~:text=RED%20DE%20RADARES%20DE%20IDEAM%20EN%20AWS).\n",
    "\n",
    "The bucket structure is s3://s3-radaresideam/l2_data/YYYY/MM/DD/Radar_name/RRRAAMMDDTTTTTT.RAWXXXX where:\n",
    "* YYYY is the 4-digit year\n",
    "* MM is the 2-digit month\n",
    "* DD is the 2-digit day\n",
    "* Radar_name radar name. Options are Guaviare, Munchique, Barrancabermja, and Carimagua\n",
    "* RRRAAMMDDTTTTTT.RAWXXXX is the radar filename with the following:\n",
    "    - RRR three first letters of the radar name (e.g., GUA for Guaviare radar)\n",
    "    - YY is the 2-digit year\n",
    "    - MM is the 2-digit month\n",
    "    - DD is the 2-digit day\n",
    "    - TTTTTT is the time at which the scan was made (GTM)\n",
    "    - RAWXXXX Sigmet file format and unique code provided by IRIS software\n",
    "    \n",
    "This is too complicated! No worries. We created a function to help you list files within the bucket."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37b38b59",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_query(date, radar_site):\n",
    "    \"\"\"\n",
    "    Creates a string for quering the IDEAM radar files stored in AWS bucket\n",
    "    :param date: date to be queried. e.g datetime(2021, 10, 3, 12). Datetime python object\n",
    "    :param radar_site: radar site e.g. Guaviare\n",
    "    :return: string with a IDEAM radar bucket format\n",
    "    \"\"\"\n",
    "    if (date.hour != 0) and (date.minute != 0):\n",
    "        return f\"l2_data/{date:%Y}/{date:%m}/{date:%d}/{radar_site}/{radar_site[:3].upper()}{date:%y%m%d%H%M}\"\n",
    "    elif (date.hour != 0) and (date.minute == 0):\n",
    "        return f\"l2_data/{date:%Y}/{date:%m}/{date:%d}/{radar_site}/{radar_site[:3].upper()}{date:%y%m%d%H}\"\n",
    "    else:\n",
    "        return f\"l2_data/{date:%Y}/{date:%m}/{date:%d}/{radar_site}/{radar_site[:3].upper()}{date:%y%m%d}\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8e1c788",
   "metadata": {},
   "source": [
    "Let's suppose we want to check the radar files on **2022-10-6** from the Guaviare radar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee045227",
   "metadata": {},
   "outputs": [],
   "source": [
    "date_query = datetime(2022, 10, 6)\n",
    "radar_name = \"Guaviare\"\n",
    "query = create_query(date=date_query, radar_site=radar_name)\n",
    "query"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecfacdd5",
   "metadata": {},
   "source": [
    "### Connecting to the AWS bucket\n",
    "Once the query is defined, we can procced to list all the available files in the bucket using **boto3** and **botocore** libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "263d95a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "str_bucket = \"s3://s3-radaresideam/\"\n",
    "s3 = boto3.resource(\n",
    "    \"s3\",\n",
    "    config=Config(signature_version=botocore.UNSIGNED, user_agent_extra=\"Resource\"),\n",
    ")\n",
    "\n",
    "bucket = s3.Bucket(\"s3-radaresideam\")\n",
    "\n",
    "radar_files = [f\"{str_bucket}{i.key}\" for i in bucket.objects.filter(Prefix=f\"{query}\")]\n",
    "radar_files[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3993934-42d3-499d-b044-21a8cb8a608f",
   "metadata": {},
   "source": [
    "We can use the Filesystem interfaces for Python [fsspec](https://filesystem-spec.readthedocs.io/en/latest/) to access the data from the s3 bucket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc357556",
   "metadata": {},
   "outputs": [],
   "source": [
    "file = fsspec.open_local(\n",
    "    f\"simplecache::{radar_files[0]}\",\n",
    "    s3={\"anon\": True},\n",
    "    filecache={\"cache_storage\": \".\"},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97ed695c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset(file, engine=\"iris\", group=\"sweep_0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4ba117d",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe456548",
   "metadata": {},
   "source": [
    "## Reflectivity and Correlation coefficient plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed4fba5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(10, 4), dpi=120)\n",
    "ds.DBZH.plot(cmap=\"ChaseSpectral\", vmin=-10, vmax=50, ax=ax)\n",
    "ds.RHOHV.plot(cmap=\"ChaseSpectral\", vmin=0, vmax=1, ax=ax1)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "278eecc3",
   "metadata": {},
   "source": [
    "The dataset object has range and azimuth as coordinates. To create a polar plot, we need to add the georeference information using  `xd.georeference.get_x_y_z()` module from [xradar](https://docs.openradarscience.org/projects/xradar/en/stable/index.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3078aa6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xd.georeference.get_x_y_z(ds)\n",
    "display(ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d779053a",
   "metadata": {},
   "source": [
    "Now x, y, and z have been added to the dataset coordinates. Let's create the new plot using the georeference information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85bb2051",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(10, 4))\n",
    "ds.DBZH.plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=-10, vmax=50, ax=ax)\n",
    "ds.RHOHV.plot(x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=0, vmax=1, ax=ax1)\n",
    "ax.set_title(\"\")\n",
    "ax1.set_title(\"\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb18393",
   "metadata": {},
   "source": [
    "## Filtering data\n",
    "\n",
    "The blue background color indicates that the radar reflectivity is less than -10 dBZ. we can filter radar data using [xarray.where](https://docs.xarray.dev/en/stable/generated/xarray.where.html) module\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "834b156e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(10, 4))\n",
    "ds.DBZH.where(ds.DBZH >= -10).plot(\n",
    "    x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=-10, vmax=50, ax=ax\n",
    ")\n",
    "ds.RHOHV.where(ds.DBZH >= -10).plot(\n",
    "    x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=0, vmax=1, ax=ax1\n",
    ")\n",
    "ax.set_title(\"\")\n",
    "ax1.set_title(\"\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cb866a8",
   "metadata": {},
   "source": [
    "Polarimetric variables can also be used as indicators to remove different noises from different sources. For example, the $\\rho_{HV}$ measures the consistency of the shapes and sizes of targets within the radar beam. Thus, the greater the $\\rho_{HV}$, the more consistent the measurement. For this example we can use $\\rho_{HV} > 0.80$ as an acceptable threshold "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80242e59-746d-463f-a667-246e3f37d16f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(10, 4))\n",
    "ds.DBZH.where(ds.DBZH >= -10).where(ds.RHOHV >= 0.80).plot(\n",
    "    x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=-10, vmax=50, ax=ax\n",
    ")\n",
    "\n",
    "ds.DBZH.where(ds.DBZH >= -10).where(ds.RHOHV >= 0.85).plot(\n",
    "    x=\"x\", y=\"y\", cmap=\"ChaseSpectral\", vmin=0, vmax=1, ax=ax1\n",
    ")\n",
    "ax.set_title(\"\")\n",
    "ax1.set_title(\"\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad8060e0",
   "metadata": {},
   "source": [
    "##  Axis labels and titles\n",
    "We can change some axis labels as well as the colorbar label "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b48519af",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax, ax1) = plt.subplots(1, 2, figsize=(10, 4))\n",
    "ds.DBZH.where(ds.DBZH >= -10).where(ds.RHOHV >= 0.85).plot(\n",
    "    x=\"x\",\n",
    "    y=\"y\",\n",
    "    cmap=\"ChaseSpectral\",\n",
    "    vmin=-10,\n",
    "    vmax=50,\n",
    "    ax=ax,\n",
    "    cbar_kwargs={\"label\": \"$Reflectivity \\ [dBZ]$\"},\n",
    ")\n",
    "\n",
    "ds.RHOHV.where(ds.DBZH >= -10).where(ds.RHOHV >= 0.85).plot(\n",
    "    x=\"x\",\n",
    "    y=\"y\",\n",
    "    cmap=\"ChaseSpectral\",\n",
    "    vmin=0,\n",
    "    vmax=1,\n",
    "    ax=ax1,\n",
    "    cbar_kwargs={\"label\": \"$Corr. \\ Coef. \\  [unitless]$\"},\n",
    ")\n",
    "\n",
    "# lambda fucntion for unit trasformation m->km\n",
    "m2km = lambda x, _: f\"{x/1000:g}\"\n",
    "# set new ticks\n",
    "ax.xaxis.set_major_formatter(m2km)\n",
    "ax.yaxis.set_major_formatter(m2km)\n",
    "ax1.xaxis.set_major_formatter(m2km)\n",
    "ax1.yaxis.set_major_formatter(m2km)\n",
    "# removing the title in both plots\n",
    "ax.set_title(\"\")\n",
    "ax1.set_title(\"\")\n",
    "\n",
    "# renaming the axis\n",
    "ax.set_ylabel(\"$North - South \\ distance \\ [km]$\")\n",
    "ax.set_xlabel(\"$East - West \\ distance \\ [km]$\")\n",
    "ax1.set_ylabel(\"$North - South \\ distance \\ [km]$\")\n",
    "ax1.set_xlabel(\"$East - West \\ distance \\ [km]$\")\n",
    "\n",
    "# setting up the title\n",
    "ax.set_title(\n",
    "    f\"$Guaviare \\ radar$\"\n",
    "    + \"\\n\"\n",
    "    + f\"${to_datetime(ds.time.values[0]): %Y-%m-%d - %X}$\"\n",
    "    + \"$ UTC$\"\n",
    ")\n",
    "ax1.set_title(\n",
    "    f\"$Guaviare \\ radar$\"\n",
    "    + \"\\n\"\n",
    "    + f\"${to_datetime(ds.time.values[0]): %Y-%m-%d - %X}$\"\n",
    "    + \"$ UTC$\"\n",
    ")\n",
    "fig.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}