{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Radar-Based Quantitative Precipitation Estimation (QPE)\nThis script processes radar data to compute vertical column profiles at a specific point location\nand estimate rain rates using categorized reflectivity. The rain rates are then saved as a\nNetCDF file.\n\n- The `column_profile` function extracts a vertical profile at a given latitude/longitude.\n- The `rain_rate_categorized` function computes rain rates based on categorized reflectivity.\n- The script reads multiple radar files, processes them, and visualizes the results.\n\n**Note:** If you want to increase the area considered for the vertical profile, you can increase\n`azimuth_spread` and `spatial_spread`. Each increment corresponds to the radar's azimuth resolution\nand range resolution, respectively.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(__doc__)\n\n# Author: Hamid Ali Syed (syed@rixvi.com)\n# License: BSD 3-Clause\n\nimport logging\nfrom datetime import datetime\n\nimport fsspec\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pytz\nimport xarray as xr\n\nimport pyart\n\n\ndef column_profile(\n    radar: pyart.core.Radar,\n    latitude: float = 33.8,\n    longitude: float = -88.3,\n    azimuth_spread: float = 3,\n    spatial_spread: float = 3,\n    v_res: float = 100,\n    min_alt: float = None,\n    max_alt: float = None,\n) -> xr.Dataset:\n    \"\"\"\n    Extracts a vertical column profile from radar data and interpolates it onto\n    a uniform height grid.\n\n    **Expanding the Area:**\n    - Increase `azimuth_spread` to expand the azimuthal range.\n    - Increase `spatial_spread` to include more horizontal coverage.\n\n    Parameters\n    ----------\n    radar : pyart.core.Radar\n        Py-ART radar object containing volume scan data.\n    latitude : float, optional\n        Latitude of the point of interest (default is 33.8).\n    longitude : float, optional\n        Longitude of the point of interest (default is -88.3).\n    azimuth_spread : float, optional\n        Azimuthal spread in degrees around the point of interest (default is 3).\n    spatial_spread : float, optional\n        Horizontal spatial spread in kilometers for averaging (default is 3).\n    v_res : float, optional\n        Vertical resolution in meters (default is 100).\n    min_alt : float\n        Minimum altitude in meters.\n    max_alt : float\n        Maximum altitude in meters.\n\n    Returns\n    -------\n    xr.Dataset\n        Interpolated columnar vertical profile on a uniform height grid.\n    \"\"\"\n    if min_alt is None or max_alt is None:\n        raise ValueError(\"Both min_alt and max_alt must be specified.\")\n\n    col_prof = pyart.util.column_vertical_profile(\n        radar,\n        latitude=latitude,\n        longitude=longitude,\n        azimuth_spread=azimuth_spread,\n        spatial_spread=spatial_spread,\n    )\n\n    new_heights = xr.Dataset(\n        coords={\"height\": (\"height\", np.arange(min_alt, max_alt + v_res, v_res))}\n    )\n\n    return col_prof.interp_like(new_heights)\n\n\ndef rain_rate_categorized(\n    dbz, conv_threshold=40.0, a_conv=300.0, b_conv=1.4, a_strat=200.0, b_strat=1.6\n):\n    \"\"\"\n    Computes rain rate from reflectivity, categorizing into convective and stratiform regions.\n    \"\"\"\n    Z = 10.0 ** (dbz / 10.0)\n    is_conv = dbz >= conv_threshold\n    rr = xr.zeros_like(Z)\n    rr = xr.where(is_conv, (Z / a_conv) ** (1.0 / b_conv), rr)\n    rr = xr.where(~is_conv, (Z / a_strat) ** (1.0 / b_strat), rr)\n    rr.attrs.update(\n        {\n            \"standard_name\": \"Rain rate\",\n            \"units\": \"mm/h\",\n            \"description\": \"Rain rate calculated using categorized reflectivity\",\n        }\n    )\n    return rr\n\n\ndef download_nexrad(timezone, date, site, local_date=False):\n    \"\"\"Download NEXRAD radar data from an S3 bucket.\"\"\"\n    try:\n        utc_date = (\n            pytz.timezone(timezone).localize(date).astimezone(pytz.utc)\n            if local_date\n            else date\n        )\n        logging.info(f\"Time: {utc_date}\")\n        fs = fsspec.filesystem(\"s3\", anon=True)\n        nexrad_path = utc_date.strftime(\n            f\"s3://unidata-nexrad-level2/%Y/%m/%d/{site}/{site}%Y%m%d_%H*\"\n        )\n        files = sorted(fs.glob(nexrad_path))\n        return [file for file in files if not file.endswith(\"_MDM\")]\n    except Exception as e:\n        logging.error(\"Error in processing: %s\", e)\n        return []\n\n\n# Load NEXRAD data\nsite = \"KGWX\"\ntimezone = \"UTC\"\ndate = datetime(2022, 3, 31, 0, 0)\nfiles = download_nexrad(timezone, date, site, local_date=False)[:5]\n\nrain_rate_list = []\n\nfor file in files:\n    radar = pyart.io.read_nexrad_archive(\"s3://\" + file)\n    col_prof_interp = column_profile(\n        radar,\n        latitude=33.5,\n        longitude=-88.3,\n        azimuth_spread=3,\n        spatial_spread=3,\n        v_res=100,\n        min_alt=200,\n        max_alt=8000,\n    )\n    rain_rate = rain_rate_categorized(col_prof_interp[\"reflectivity\"])\n    rain_rate[\"time\"] = col_prof_interp[\"base_time\"]\n    rain_rate_list.append(rain_rate)\n\nda_rain_rate = xr.concat(rain_rate_list, dim=\"time\")\n\nda_rain_rate.plot.contourf(\n    x=\"time\", levels=range(0, 26, 2), cmap=\"RRate11\", figsize=[10, 3]\n)\nplt.title(\"QPE\")\nplt.xlabel(\"Time\")\nplt.ylabel(\"Height ASL [m]\")\nplt.show()\n\nda_rain_rate.sel(height=1000, method=\"nearest\").plot.line(\"k-\", figsize=[10, 3])\nplt.title(\"QPE\")\nplt.xlabel(\"Time\")\nplt.show()\n\n# References\n# Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological\n# Characterization of Three-Dimensional Storm Structure from Operational\n# Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007.\n# https://doi.org/10.1175/1520-0450(1995)034<1978:CCOTDS>2.0.CO;2."
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "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.14.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}