{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Convective-Stratiform classification\nThis example shows how to use the updated convective stratiform classifcation algorithm. We show 3 examples,\na summer convective example, an example from Hurricane Ian, and an example from a winter storm.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(__doc__)\n\n# Author: Laura Tomkins (lmtomkin@ncsu.edu)\n# License: BSD 3 clause\n\n\nimport cartopy.crs as ccrs\nimport matplotlib.colors as mcolors\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport pyart"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## How the algorithm works\nThis first section describes how the convective-stratiform algorithm works (see references for full details). This\nalgorithm is a feature detection algorithm and classifies fields as \"convective\" or \"stratiform\". The algorithm is\ndesigned to detect features in a reflectivity field but can also detect features in fields such as rain rate or\nsnow rate. In this section we describe the steps of the convective stratiform algorithm and the variables used in\nthe function.\nThe first step of the algorithm calculates a background average of the field with a circular footprint using a radius\nprovided by ``bkg_rad_km``. A larger radius will yield a smoother field. The radius needs to be at least double the\ngrid spacing, but we recommend at least three times the grid spacing. If using reflectivity, ``dB_averaging`` should be set\nto True to convert reflectivity to linear Z before averaging, False for rescaled fields such as rain or snow rate.\n``calc_thres`` determines the minimum fraction of a circle that is considered in the background average calculation\n(default is 0.75, so the points along the edges where there is less than 75% of a full circle of data,\nthe algorithm is not run).\nOnce the background average has been calculated, the original field is compared to the background average.  In\norder for points to be considered \"convective cores\" they must exceed the background value by a certain value or\nsimply be greater than the ``always_core_thres``. This value is determined by either a cosine scheme, or a scalar\nvalue (i.e. the reflectivity value must be X times the background value (multiplier; ``use_addition=False``),\nor  X greater than the background value (``use_addition=True``) where X is the ``scalar_diff``).\n``use_cosine`` determines if a cosine scheme or scalar scheme is to be used. If ``use_cosine`` is True,\nthen the ``max_diff`` and ``zero_diff_cos_val`` come into use. These values define the cosine scheme that is used  to\ndetermine the minimum difference between the background and reflectivity value in order for a core to be\nidentified. ``max_diff`` is the maximum difference between the field and the background for a core to be identified,\nor where the cosine function crosses the y-axis. The ``zero_diff_cos_val`` is where the difference between the  field\nand the background is zero, or where the cosine function crosses the x-axis. Note, if\n``always_core_thres`` < ``zero_diff_cos_val``, ``zero_diff_cos_val`` only helps define the shape of the cosine curve and\nall values greater than ``always_core_thres`` will be considered a convective core. If\n``always_core_thres`` > ``zero_diff_cos_val`` then all values greater than ``zero_diff_cos_val`` will be considered a\nconvective core. We plot some examples of the schemes below:\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Example of the cosine scheme:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pyart.graph.plot_convstrat_scheme(\n    always_core_thres=30, use_cosine=True, max_diff=5, zero_diff_cos_val=45\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "when zero_diff_cos_val is greater than always_core_thres, the difference becomes zero at the zero_diff_cos_val\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pyart.graph.plot_convstrat_scheme(\n    always_core_thres=55, use_cosine=True, max_diff=5, zero_diff_cos_val=45\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "alternatively, we can use a simpler scalar difference instead of a cosine scheme\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pyart.graph.plot_convstrat_scheme(\n    always_core_thres=40,\n    use_cosine=False,\n    max_diff=None,\n    zero_diff_cos_val=None,\n    use_addition=True,\n    scalar_diff=2,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "if you are interested in picking up weak features, you can also use the scalar difference as a multiplier instead,\nso very weak features do not have to be that different from the background to be classified as convective.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pyart.graph.plot_convstrat_scheme(\n    always_core_thres=40,\n    use_cosine=False,\n    max_diff=None,\n    zero_diff_cos_val=None,\n    use_addition=False,\n    scalar_diff=2,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Once the cores are identified, there is an option to remove speckles (``remove_small_objects``) smaller than a  given\nsize (``min_km2_size``).\nAfter the convective cores are identified, We then incorporate convective radii using\n``val_for_max_conv_rad`` and ``max_conv_rad_km``. The convective radii act as a dilation and are used to classify\nadditional points around the cores as convective that may not have been identified previously.  The\n``val_for_max_conv_rad`` is the value where the maximum convective radius is applied and the ``max_conv_rad_km`` is the\nmaximum convective radius. Values less than the ``val_for_max_conv_rad`` are assigned a convective radius using a step\nfunction.\nFinally, the points are classified as NOSFCECHO (threshold set with ``min_dBZ_used``; 0), WEAKECHO (threshold set with\n``weak_echo_thres``; 3), SF (stratiform; 1), CONV (convective; 2).\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Examples\n**Classification of summer convective example**\n\nOur first example classifies echo from a summer convective event. We use a cosine scheme to classify the convective\npoints.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Now let's do a classification with our parameters\n# read in file\nfilename = pyart.testing.get_test_data(\"swx_20120520_0641.nc\")\nradar = pyart.io.read(filename)\n\n# extract the lowest sweep\nradar = radar.extract_sweeps([0])\n\n# interpolate to grid\ngrid = pyart.map.grid_from_radars(\n    (radar,),\n    grid_shape=(1, 201, 201),\n    grid_limits=((0, 10000), (-50000.0, 50000.0), (-50000.0, 50000.0)),\n    fields=[\"reflectivity_horizontal\"],\n)\n\n# get dx dy\ndx = grid.x[\"data\"][1] - grid.x[\"data\"][0]\ndy = grid.y[\"data\"][1] - grid.y[\"data\"][0]\n\n# convective stratiform classification\nconvsf_dict = pyart.retrieve.conv_strat_yuter(\n    grid,\n    dx,\n    dy,\n    refl_field=\"reflectivity_horizontal\",\n    always_core_thres=40,\n    bkg_rad_km=20,\n    use_cosine=True,\n    max_diff=5,\n    zero_diff_cos_val=55,\n    weak_echo_thres=10,\n    max_conv_rad_km=2,\n)\n\n# add to grid object\n# mask zero values (no surface echo)\nconvsf_masked = np.ma.masked_equal(convsf_dict[\"feature_detection\"][\"data\"], 0)\n# mask 3 values (weak echo)\nconvsf_masked = np.ma.masked_equal(convsf_masked, 3)\n# add dimension to array to add to grid object\nconvsf_dict[\"feature_detection\"][\"data\"] = convsf_masked\n# add field\ngrid.add_field(\"convsf\", convsf_dict[\"feature_detection\"], replace_existing=True)\n\n# create plot using GridMapDisplay\n# plot variables\ndisplay = pyart.graph.GridMapDisplay(grid)\nmagma_r_cmap = plt.get_cmap(\"magma_r\")\nref_cmap = mcolors.LinearSegmentedColormap.from_list(\n    \"ref_cmap\", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))\n)\nprojection = ccrs.AlbersEqualArea(\n    central_latitude=radar.latitude[\"data\"][0],\n    central_longitude=radar.longitude[\"data\"][0],\n)\n\n# plot\nplt.figure(figsize=(10, 4))\nax1 = plt.subplot(1, 2, 1, projection=projection)\ndisplay.plot_grid(\n    \"reflectivity_horizontal\",\n    vmin=5,\n    vmax=45,\n    cmap=ref_cmap,\n    transform=ccrs.PlateCarree(),\n    ax=ax1,\n)\nax2 = plt.subplot(1, 2, 2, projection=projection)\ndisplay.plot_grid(\n    \"convsf\",\n    vmin=0,\n    vmax=2,\n    cmap=plt.get_cmap(\"viridis\", 3),\n    ax=ax2,\n    transform=ccrs.PlateCarree(),\n    ticks=[1 / 3, 1, 5 / 3],\n    ticklabs=[\"\", \"Stratiform\", \"Convective\"],\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In addition to the default convective-stratiform classification, the function also returns an underestimate\n(convsf_under) and an overestimate (convsf_over) to take into consideration the uncertainty when choosing\nclassification parameters. The under and overestimate use the same parameters, but vary the input field by a\ncertain value (default is 5 dBZ, can be changed with ``estimate_offset``). The estimation can be turned off (\n``estimate_flag=False``), but we recommend keeping it turned on.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# mask weak echo and no surface echo\nconvsf_masked = np.ma.masked_equal(convsf_dict[\"feature_detection\"][\"data\"], 0)\nconvsf_masked = np.ma.masked_equal(convsf_masked, 3)\nconvsf_dict[\"feature_detection\"][\"data\"] = convsf_masked\n# underest.\nconvsf_masked = np.ma.masked_equal(convsf_dict[\"feature_under\"][\"data\"], 0)\nconvsf_masked = np.ma.masked_equal(convsf_masked, 3)\nconvsf_dict[\"feature_under\"][\"data\"] = convsf_masked\n# overest.\nconvsf_masked = np.ma.masked_equal(convsf_dict[\"feature_over\"][\"data\"], 0)\nconvsf_masked = np.ma.masked_equal(convsf_masked, 3)\nconvsf_dict[\"feature_over\"][\"data\"] = convsf_masked\n\n# Plot each estimation\nplt.figure(figsize=(10, 4))\nax1 = plt.subplot(131)\nax1.pcolormesh(\n    convsf_dict[\"feature_detection\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=2,\n    cmap=plt.get_cmap(\"viridis\", 3),\n)\nax1.set_title(\"Best estimate\")\nax1.set_aspect(\"equal\")\nax2 = plt.subplot(132)\nax2.pcolormesh(\n    convsf_dict[\"feature_under\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=2,\n    cmap=plt.get_cmap(\"viridis\", 3),\n)\nax2.set_title(\"Underestimate\")\nax2.set_aspect(\"equal\")\nax3 = plt.subplot(133)\nax3.pcolormesh(\n    convsf_dict[\"feature_over\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=2,\n    cmap=plt.get_cmap(\"viridis\", 3),\n)\nax3.set_title(\"Overestimate\")\nax3.set_aspect(\"equal\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Tropical example**\n\nLet's get a NEXRAD file from Hurricane Ian\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Read in file\nnexrad_file = \"s3://noaa-nexrad-level2/2022/09/28/KTBW/KTBW20220928_190142_V06\"\nradar = pyart.io.read_nexrad_archive(nexrad_file)\n\n# extract the lowest sweep\nradar = radar.extract_sweeps([0])\n\n# interpolate to grid\ngrid = pyart.map.grid_from_radars(\n    (radar,),\n    grid_shape=(1, 201, 201),\n    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),\n    fields=[\"reflectivity\"],\n)\n\n# get dx dy\ndx = grid.x[\"data\"][1] - grid.x[\"data\"][0]\ndy = grid.y[\"data\"][1] - grid.y[\"data\"][0]\n\n# convective stratiform classification\nconvsf_dict = pyart.retrieve.conv_strat_yuter(\n    grid,\n    dx,\n    dy,\n    refl_field=\"reflectivity\",\n    always_core_thres=40,\n    bkg_rad_km=20,\n    use_cosine=True,\n    max_diff=3,\n    zero_diff_cos_val=55,\n    weak_echo_thres=5,\n    max_conv_rad_km=2,\n    estimate_flag=False,\n)\n\n# add to grid object\n# mask zero values (no surface echo)\nconvsf_masked = np.ma.masked_equal(convsf_dict[\"feature_detection\"][\"data\"], 0)\n# mask 3 values (weak echo)\nconvsf_masked = np.ma.masked_equal(convsf_masked, 3)\n# add dimension to array to add to grid object\nconvsf_dict[\"feature_detection\"][\"data\"] = convsf_masked\n# add field\ngrid.add_field(\"convsf\", convsf_dict[\"feature_detection\"], replace_existing=True)\n\n# create plot using GridMapDisplay\n# plot variables\ndisplay = pyart.graph.GridMapDisplay(grid)\nmagma_r_cmap = plt.get_cmap(\"magma_r\")\nref_cmap = mcolors.LinearSegmentedColormap.from_list(\n    \"ref_cmap\", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))\n)\nprojection = ccrs.AlbersEqualArea(\n    central_latitude=radar.latitude[\"data\"][0],\n    central_longitude=radar.longitude[\"data\"][0],\n)\n# plot\nplt.figure(figsize=(10, 4))\nax1 = plt.subplot(1, 2, 1, projection=projection)\ndisplay.plot_grid(\n    \"reflectivity\",\n    vmin=5,\n    vmax=45,\n    cmap=ref_cmap,\n    transform=ccrs.PlateCarree(),\n    ax=ax1,\n    axislabels_flag=False,\n)\nax2 = plt.subplot(1, 2, 2, projection=projection)\ndisplay.plot_grid(\n    \"convsf\",\n    vmin=0,\n    vmax=2,\n    cmap=plt.get_cmap(\"viridis\", 3),\n    axislabels_flag=False,\n    transform=ccrs.PlateCarree(),\n    ticks=[1 / 3, 1, 5 / 3],\n    ticklabs=[\"\", \"Stratiform\", \"Convective\"],\n    ax=ax2,\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Winter storm example with image muting**\n\nHere is a final example of the convective stratiform classification using an example from a winter storm. Before\ndoing the classification, we image mute the reflectivity to remove regions with melting or mixed precipitation. We\nthen rescale the reflectivity to snow rate (Rasumussen et al. 2003). We recommend using a rescaled reflectivity\nto do the classification, but if you do make sure to changed dB_averaging to False because this parameter is  used\nto convert reflectivity to a linear value before averaging (set dB_averaging to True for reflectivity fields in\ndBZ units).\nIn this example, note how we change some of the other parameters since we are classifying snow rate instead of\nreflecitivity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Read in file\nnexrad_file = \"s3://noaa-nexrad-level2/2021/02/07/KOKX/KOKX20210207_161413_V06\"\nradar = pyart.io.read_nexrad_archive(nexrad_file)\n\n# extract the lowest sweep\nradar = radar.extract_sweeps([0])\n\n# interpolate to grid\ngrid = pyart.map.grid_from_radars(\n    (radar,),\n    grid_shape=(1, 201, 201),\n    grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)),\n    fields=[\"reflectivity\", \"cross_correlation_ratio\"],\n)\n\n# image mute grid object\ngrid = pyart.util.image_mute_radar(\n    grid, \"reflectivity\", \"cross_correlation_ratio\", 0.97, 20\n)\n\n# convect non-muted reflectivity to snow rate\nnonmuted_ref = grid.fields[\"nonmuted_reflectivity\"][\"data\"][0, :, :]\nnonmuted_ref = np.ma.masked_invalid(nonmuted_ref)\n\nnonmuted_ref_linear = 10 ** (nonmuted_ref / 10)  # mm6/m3\nsnow_rate = (nonmuted_ref_linear / 57.3) ** (1 / 1.67)  #\n\n# add to grid\nsnow_rate_dict = {\n    \"data\": snow_rate[None, :, :],\n    \"standard_name\": \"snow_rate\",\n    \"long_name\": \"Snow rate converted from linear reflectivity\",\n    \"units\": \"mm/hr\",\n    \"valid_min\": 0,\n    \"valid_max\": 40500,\n}\ngrid.add_field(\"snow_rate\", snow_rate_dict, replace_existing=True)\n\n# get dx dy\ndx = grid.x[\"data\"][1] - grid.x[\"data\"][0]\ndy = grid.y[\"data\"][1] - grid.y[\"data\"][0]\n\n# convective stratiform classification\nconvsf_dict = pyart.retrieve.conv_strat_yuter(\n    grid,\n    dx,\n    dy,\n    refl_field=\"snow_rate\",\n    dB_averaging=False,\n    always_core_thres=4,\n    bkg_rad_km=40,\n    use_cosine=True,\n    max_diff=1.5,\n    zero_diff_cos_val=5,\n    weak_echo_thres=0,\n    min_dBZ_used=0,\n    max_conv_rad_km=1,\n    estimate_flag=False,\n)\n\n# add to grid object\n# mask zero values (no surface echo)\nconvsf_masked = np.ma.masked_equal(convsf_dict[\"feature_detection\"][\"data\"], 0)\n# mask 3 values (weak echo)\nconvsf_masked = np.ma.masked_equal(convsf_masked, 3)\n# add dimension to array to add to grid object\nconvsf_dict[\"feature_detection\"][\"data\"] = convsf_masked\n# add field\ngrid.add_field(\"convsf\", convsf_dict[\"feature_detection\"], replace_existing=True)\n\n# create plot using GridMapDisplay\n# plot variables\ndisplay = pyart.graph.GridMapDisplay(grid)\nmagma_r_cmap = plt.get_cmap(\"magma_r\")\nref_cmap = mcolors.LinearSegmentedColormap.from_list(\n    \"ref_cmap\", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N))\n)\nprojection = ccrs.AlbersEqualArea(\n    central_latitude=radar.latitude[\"data\"][0],\n    central_longitude=radar.longitude[\"data\"][0],\n)\n# plot\nplt.figure(figsize=(10, 4))\nax1 = plt.subplot(1, 2, 1, projection=projection)\ndisplay.plot_grid(\n    \"snow_rate\",\n    vmin=0,\n    vmax=10,\n    cmap=plt.get_cmap(\"viridis\"),\n    transform=ccrs.PlateCarree(),\n    ax=ax1,\n    axislabels_flag=False,\n)\nax2 = plt.subplot(1, 2, 2, projection=projection)\ndisplay.plot_grid(\n    \"convsf\",\n    vmin=0,\n    vmax=2,\n    cmap=plt.get_cmap(\"viridis\", 3),\n    axislabels_flag=False,\n    transform=ccrs.PlateCarree(),\n    ticks=[1 / 3, 1, 5 / 3],\n    ticklabs=[\"\", \"Stratiform\", \"Convective\"],\n    ax=ax2,\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Summary of recommendations and best practices\n* Tune your parameters to your specific purpose\n* Use a rescaled field if possible (i.e. linear reflectivity, rain or snow rate)\n* Keep ``estimate_flag=True`` to see uncertainty in classification\n\n## References\nSteiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological\nCharacterization of Three-Dimensional Storm Structure from Operational\nRadar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007.\nhttps://doi.org/10.1175/1520-0450(1995)034<1978:CCOTDS>2.0.CO;2.\n\nYuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size\ndistributions over the Pacific warm pool and implications for Z-R relations.\nJ. Appl. Meteor., 36, 847-867.\nhttps://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2\n\nYuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser,\n2005: Physical characterization of tropical oceanic convection observed in\nKWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1\n\nRasmussen, R., M. Dixon, S. Vasiloff, F. Hage, S. Knight, J. Vivekanandan,\nand M. Xu, 2003: Snow Nowcasting Using a Real-Time Correlation of Radar\nReflectivity with Snow Gauge Accumulation. J. Appl. Meteorol. Climatol., 42, 20\u201336.\nhttps://doi.org/10.1175/1520-0450(2003)042%3C0020:SNUART%3E2.0.CO;2\n\n"
      ]
    }
  ],
  "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.13.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}