{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Feature Detection classification\nThis example shows how to use the feature detection algorithm. We show how to use the algorithm to detect convective\nand stratiform regions in warm season events and how the algorithm can be used to detect both weak and strong\nfeatures in cool-season events.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(__doc__)\n\n# Author: Laura Tomkins (lauramtomkins@gmail.com)\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\nfrom open_radar_data import DATASETS\n\nimport pyart"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## How the algorithm works\nThe feature detection algorithm works by identifying features that exceed the background value by an amount that\nvaries with the background value. The algorithm is heavily customizable and is designed to work with a variety of\ndatasets. Here, we show several examples of how to use the algorithm with different types of radar data.\n\nSee Steiner et al. (1995), Yuter and Houze (1997), and Yuter et al. (2005) for full details on the algorithm. Tomkins\net al. 2024 builds on this work to describe feature detection in cool-season events (part 2).\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Part 1: Warm-season convective-stratiform classification\n**Classification of summer convective example**\n\nOur first example classifies echo from a summer convective event.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# 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    nb=1.5,\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# feature detection\nfeature_dict = pyart.retrieve.feature_detection(\n    grid,\n    dx,\n    dy,\n    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_rad_km=2,\n)\n\n# add to grid object\n# mask zero values (no surface echo)\nfeature_masked = np.ma.masked_equal(feature_dict[\"feature_detection\"][\"data\"], 0)\n# mask 3 values (weak echo)\nfeature_masked = np.ma.masked_equal(feature_masked, 3)\n# add dimension to array to add to grid object\nfeature_dict[\"feature_detection\"][\"data\"] = feature_masked\n# add field\ngrid.add_field(\n    \"feature_detection\", feature_dict[\"feature_detection\"], replace_existing=True\n)\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    \"feature_detection\",\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 feature detection classification, the function also returns an underestimate\n(``feature_under``) and an overestimate (``feature_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\nfeature_masked = np.ma.masked_equal(feature_dict[\"feature_detection\"][\"data\"], 0)\nfeature_masked = np.ma.masked_equal(feature_masked, 3)\nfeature_dict[\"feature_detection\"][\"data\"] = feature_masked\n# underest.\nfeature_masked = np.ma.masked_equal(feature_dict[\"feature_under\"][\"data\"], 0)\nfeature_masked = np.ma.masked_equal(feature_masked, 3)\nfeature_dict[\"feature_under\"][\"data\"] = feature_masked\n# overest.\nfeature_masked = np.ma.masked_equal(feature_dict[\"feature_over\"][\"data\"], 0)\nfeature_masked = np.ma.masked_equal(feature_masked, 3)\nfeature_dict[\"feature_over\"][\"data\"] = feature_masked\n\n# Plot each estimation\nplt.figure(figsize=(10, 4))\nax1 = plt.subplot(131)\nax1.pcolormesh(\n    feature_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    feature_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    feature_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://unidata-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    nb=1.5,\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# feature detection\nfeature_dict = pyart.retrieve.feature_detection(\n    grid,\n    dx,\n    dy,\n    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_rad_km=2,\n    estimate_flag=False,\n)\n\n# add to grid object\n# mask zero values (no surface echo)\nfeature_masked = np.ma.masked_equal(feature_dict[\"feature_detection\"][\"data\"], 0)\n# mask 3 values (weak echo)\nfeature_masked = np.ma.masked_equal(feature_masked, 3)\n# add dimension to array to add to grid object\nfeature_dict[\"feature_detection\"][\"data\"] = feature_masked\n# add field\ngrid.add_field(\n    \"feature_detection\", feature_dict[\"feature_detection\"], replace_existing=True\n)\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    \"feature_detection\",\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": [
        "**Comparison with original C++ algorithm**\n\nTo ensure that our algorithm here matches the original algorithm developed in C++, we compare the results from an\nexample during the KWAJEX field campaign. The file contains the original convective-stratiform classification from\nthe C++ algorithm. We use the algorithm implemented in Py-ART with the same settings to produce identical results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# read in file\nfilename = DATASETS.fetch(\"convsf.19990811.221202.cdf\")\ngrid = pyart.io.read_grid(filename)\n\n# colormap\nconvsf_cmap = mcolors.LinearSegmentedColormap.from_list(\n    \"convsf\",\n    [\n        (33 / 255, 145 / 255, 140 / 255),\n        (253 / 255, 231 / 255, 37 / 255),\n        (210 / 255, 180 / 255, 140 / 255),\n    ],\n    N=3,\n)\n\n# get original data from file processed in C++\nx = grid.x[\"data\"]\ny = grid.y[\"data\"]\nref = grid.fields[\"maxdz\"][\"data\"][0, :, :]\ncsb = grid.fields[\"convsf\"][\"data\"][0, :, :]\ncsb_masked = np.ma.masked_less_equal(csb, 0)\ncsl = grid.fields[\"convsf_lo\"][\"data\"][0, :, :]\ncsl_masked = np.ma.masked_less_equal(csl, 0)\ncsh = grid.fields[\"convsf_hi\"][\"data\"][0, :, :]\ncsh_masked = np.ma.masked_less_equal(csh, 0)\n\n# plot\nfig, axs = plt.subplots(2, 2, figsize=(12, 10))\n# reflectivity\nrpm = axs[0, 0].pcolormesh(x / 1000, y / 1000, ref, vmin=0, vmax=45, cmap=\"magma_r\")\naxs[0, 0].set_aspect(\"equal\")\naxs[0, 0].set_title(\"Reflectivity [dBZ]\")\nplt.colorbar(rpm, ax=axs[0, 0])\n# convsf\ncsbpm = axs[0, 1].pcolormesh(\n    x / 1000, y / 1000, csb_masked, vmin=1, vmax=3, cmap=convsf_cmap\n)\naxs[0, 1].set_aspect(\"equal\")\naxs[0, 1].set_title(\"convsf\")\ncb = plt.colorbar(csbpm, ax=axs[0, 1], ticks=[4 / 3, 2, 8 / 3])\ncb.ax.set_yticklabels([\"Stratiform\", \"Convective\", \"Weak Echo\"])\n# convsf lo\ncslpm = axs[1, 0].pcolormesh(\n    x / 1000, y / 1000, csl_masked, vmin=1, vmax=3, cmap=convsf_cmap\n)\naxs[1, 0].set_aspect(\"equal\")\naxs[1, 0].set_title(\"convsf lo\")\ncb = plt.colorbar(cslpm, ax=axs[1, 0], ticks=[])\n# convsf hi\ncshpm = axs[1, 1].pcolormesh(\n    x / 1000, y / 1000, csh_masked, vmin=1, vmax=3, cmap=convsf_cmap\n)\naxs[1, 1].set_aspect(\"equal\")\naxs[1, 1].set_title(\"convsf hi\")\ncb = plt.colorbar(cshpm, ax=axs[1, 1], ticks=[4 / 3, 2, 8 / 3])\ncb.ax.set_yticklabels([\"Stratiform\", \"Convective\", \"Weak Echo\"])\nplt.show()\n\n# now let's compare with the Python algorithm\nconvsf_dict = pyart.retrieve.feature_detection(\n    grid,\n    dx=x[1] - x[0],\n    dy=y[1] - y[0],\n    level_m=None,\n    always_core_thres=40,\n    bkg_rad_km=11,\n    use_cosine=True,\n    max_diff=8,\n    zero_diff_cos_val=55,\n    scalar_diff=None,\n    use_addition=False,\n    calc_thres=0,\n    weak_echo_thres=15,\n    min_val_used=0,\n    dB_averaging=True,\n    remove_small_objects=False,\n    min_km2_size=None,\n    val_for_max_rad=30,\n    max_rad_km=5.0,\n    core_val=3,\n    nosfcecho=0,\n    weakecho=3,\n    bkgd_val=1,\n    feat_val=2,\n    field=\"maxdz\",\n    estimate_flag=True,\n    estimate_offset=5,\n)\n\n# new data\ncsb_lt = convsf_dict[\"feature_detection\"][\"data\"][0, :, :]\ncsb_lt_masked = np.ma.masked_less_equal(csb_lt, 0)\ncsu_lt = convsf_dict[\"feature_under\"][\"data\"][0, :, :]\ncsu_lt_masked = np.ma.masked_less_equal(csu_lt, 0)\ncso_lt = convsf_dict[\"feature_over\"][\"data\"][0, :, :]\ncso_lt_masked = np.ma.masked_less_equal(cso_lt, 0)\n\nfig, axs = plt.subplots(2, 2, figsize=(12, 10))\n# reflectivity\nrpm = axs[0, 0].pcolormesh(x / 1000, y / 1000, ref, vmin=0, vmax=45, cmap=\"magma_r\")\naxs[0, 0].set_aspect(\"equal\")\naxs[0, 0].set_title(\"Reflectivity [dBZ]\")\nplt.colorbar(rpm, ax=axs[0, 0])\n# convsf best\ncsbpm = axs[0, 1].pcolormesh(\n    x / 1000, y / 1000, csb_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap\n)\naxs[0, 1].set_aspect(\"equal\")\naxs[0, 1].set_title(\"convsf\")\ncb = plt.colorbar(csbpm, ax=axs[0, 1], ticks=[4 / 3, 2, 8 / 3])\ncb.ax.set_yticklabels([\"Stratiform\", \"Convective\", \"Weak Echo\"])\n# convsf under\ncsupm = axs[1, 0].pcolormesh(\n    x / 1000, y / 1000, csu_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap\n)\naxs[1, 0].set_aspect(\"equal\")\naxs[1, 0].set_title(\"convsf under\")\ncb = plt.colorbar(csupm, ax=axs[1, 0], ticks=[])\n# convsf over\ncsopm = axs[1, 1].pcolormesh(\n    x / 1000, y / 1000, cso_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap\n)\naxs[1, 1].set_aspect(\"equal\")\naxs[1, 1].set_title(\"convsf over\")\ncb = plt.colorbar(csopm, ax=axs[1, 1], ticks=[4 / 3, 2, 8 / 3])\ncb.ax.set_yticklabels([\"Stratiform\", \"Convective\", \"Weak Echo\"])\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Part 2: Cool-season feature detection\n**Winter storm example**\n\nIn this example, we will show how to algorithm can be used to detect features (snow bands) in winter storms. Here, we\nrescale the reflectivity to snow rate (Rasumussen et al. 2003) in order to better represent the snow field. Note\nthat we are not using the absolute values of the snow rate, we are only using the relative anomalies relative to\nthe background average. We recommend using a rescaled reflectivity for this algorithm, but note that you need to\nchange ``dB_averaging = False`` when using a rescaled field (set ``dB_averaging`` to True for reflectivity fields in\ndBZ units). Also note that since we are now using a snow rate field with different values, most of our parameters\nhave changed as well.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Read in file\nnexrad_file = \"s3://unidata-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    nb=1.5,\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# rescale reflectivity to snow rate\ngrid = pyart.retrieve.qpe.ZtoR(\n    grid, ref_field=\"reflectivity\", a=57.3, b=1.67, save_name=\"snow_rate\"\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# feature detection\nfeature_dict = pyart.retrieve.feature_detection(\n    grid,\n    dx,\n    dy,\n    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_val_used=0,\n    max_rad_km=1,\n    estimate_flag=False,\n)\n\n# add to grid object\n# mask zero values (no surface echo)\nfeature_masked = np.ma.masked_equal(feature_dict[\"feature_detection\"][\"data\"], 0)\n# mask 3 values (weak echo)\nfeature_masked = np.ma.masked_equal(feature_masked, 3)\n# add dimension to array to add to grid object\nfeature_dict[\"feature_detection\"][\"data\"] = feature_masked\n# add field\ngrid.add_field(\n    \"feature_detection\", feature_dict[\"feature_detection\"], replace_existing=True\n)\n\n# create plot using GridMapDisplay\n# plot variables\ndisplay = pyart.graph.GridMapDisplay(grid)\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    \"feature_detection\",\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=[\"\", \"Background\", \"Feature\"],\n    ax=ax2,\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Under and over estimate in snow**\n\nHere we show how to bound the snow rate with an under and over estimate\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get reflectivity field and copy\nref_field = grid.fields[\"reflectivity\"]\nref_field_over = ref_field.copy()\nref_field_under = ref_field.copy()\n\n# offset original field by +/- 2dB\nref_field_over[\"data\"] = ref_field[\"data\"] + 2\nref_field_under[\"data\"] = ref_field[\"data\"] - 2\n\n# adjust dictionary\nref_field_over[\"standard_name\"] = ref_field[\"standard_name\"] + \"_overest\"\nref_field_over[\"long_name\"] = ref_field[\"long_name\"] + \" overestimate\"\nref_field_under[\"standard_name\"] = ref_field[\"standard_name\"] + \"_underest\"\nref_field_under[\"long_name\"] = ref_field[\"long_name\"] + \" underestimate\"\n\n# add to grid object\ngrid.add_field(\"reflectivity_over\", ref_field_over, replace_existing=True)\ngrid.add_field(\"reflectivity_under\", ref_field_under, replace_existing=True)\n\n# convert to snow rate\ngrid = pyart.retrieve.qpe.ZtoR(\n    grid, ref_field=\"reflectivity_over\", a=57.3, b=1.67, save_name=\"snow_rate_over\"\n)\ngrid = pyart.retrieve.qpe.ZtoR(\n    grid, ref_field=\"reflectivity_under\", a=57.3, b=1.67, save_name=\"snow_rate_under\"\n)\n\n# now do feature detection with under and over estimate fields\nfeature_estimate_dict = pyart.retrieve.feature_detection(\n    grid,\n    dx,\n    dy,\n    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_val_used=0,\n    max_rad_km=1,\n    estimate_flag=True,\n    overest_field=\"snow_rate_over\",\n    underest_field=\"snow_rate_under\",\n)\n\n# x and y data for plotting\nx = grid.x[\"data\"]\ny = grid.y[\"data\"]\n\n# now let's plot\nfig, axs = plt.subplots(2, 2, figsize=(12, 10))\n# snow rate\nrpm = axs[0, 0].pcolormesh(\n    x / 1000,\n    y / 1000,\n    grid.fields[\"snow_rate\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=10,\n)\naxs[0, 0].set_aspect(\"equal\")\naxs[0, 0].set_title(\"Reflectivity [dBZ]\")\nplt.colorbar(rpm, ax=axs[0, 0])\n# features best\nbpm = axs[0, 1].pcolormesh(\n    x / 1000,\n    y / 1000,\n    feature_estimate_dict[\"feature_detection\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=3,\n    cmap=plt.get_cmap(\"viridis\", 3),\n)\naxs[0, 1].set_aspect(\"equal\")\naxs[0, 1].set_title(\"Feature detection (best)\")\ncb = plt.colorbar(bpm, ax=axs[0, 1], ticks=[3 / 2, 5 / 2])\ncb.ax.set_yticklabels([\"Background\", \"Feature\"])\n# features underestimate\nupm = axs[1, 0].pcolormesh(\n    x / 1000,\n    y / 1000,\n    feature_estimate_dict[\"feature_under\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=3,\n    cmap=plt.get_cmap(\"viridis\", 3),\n)\naxs[1, 0].set_aspect(\"equal\")\naxs[1, 0].set_title(\"Feature detection (underestimate)\")\ncb = plt.colorbar(upm, ax=axs[1, 0], ticks=[3 / 2, 5 / 2])\ncb.ax.set_yticklabels([\"Background\", \"Feature\"])\n# features overestimate\nopm = axs[1, 1].pcolormesh(\n    x / 1000,\n    y / 1000,\n    feature_estimate_dict[\"feature_over\"][\"data\"][0, :, :],\n    vmin=0,\n    vmax=3,\n    cmap=plt.get_cmap(\"viridis\", 3),\n)\naxs[1, 1].set_aspect(\"equal\")\naxs[1, 1].set_title(\"Feature detection (overestimate)\")\ncb = plt.colorbar(opm, ax=axs[1, 1], ticks=[3 / 2, 5 / 2])\ncb.ax.set_yticklabels([\"Background\", \"Feature\"])\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Strong and Faint features**\n\nWe have developed a new technique which runs the algorithm twice to find features that are very distinct from the\nbackground (strong features) and features that are not very distinct from the background (faint features).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# run the algorithm again, but to find objects that are also faint\n# feature detection\nfeature_dict = pyart.retrieve.feature_detection(\n    grid,\n    dx,\n    dy,\n    field=\"snow_rate\",\n    dB_averaging=False,\n    always_core_thres=4,\n    bkg_rad_km=40,\n    use_cosine=False,\n    scalar_diff=1.5,\n    use_addition=False,\n    weak_echo_thres=0,\n    min_val_used=0,\n    max_rad_km=1,\n    estimate_flag=False,\n)\n\n# separate strong vs. faint\n# mask zero values (no surface echo)\nscalar_features_masked = np.ma.masked_equal(\n    feature_dict[\"feature_detection\"][\"data\"], 0\n)\n# mask 3 values (weak echo)\nscalar_features_masked = np.ma.masked_equal(scalar_features_masked, 3)\n# get cosine features\ncosine_features_masked = grid.fields[\"feature_detection\"][\"data\"]\n# isolate faint features only\nfaint_features = scalar_features_masked - cosine_features_masked\nfaint_features_masked = np.ma.masked_less_equal(faint_features, 0)\n# add strong and faint features\nfeatures_faint_strong = 2 * faint_features_masked.filled(\n    0\n) + cosine_features_masked.filled(0)\nfeatures_faint_strong_masked = np.ma.masked_where(\n    cosine_features_masked.mask, features_faint_strong\n)\n\n# add to grid object\nnew_dict = {\n    \"features_faint_strong\": {\n        \"data\": features_faint_strong_masked,\n        \"standard_name\": \"features_faint_strong\",\n        \"long name\": \"Faint and Strong Features\",\n        \"valid_min\": 0,\n        \"valid_max\": 10,\n        \"comment_1\": \"3: Faint features, 2: Strong features, 1: background\",\n    }\n}\n# add field\ngrid.add_field(\n    \"features_faint_strong\", new_dict[\"features_faint_strong\"], replace_existing=True\n)\n\n# create plot using GridMapDisplay\n# plot variables\ndisplay = pyart.graph.GridMapDisplay(grid)\n\nprojection = ccrs.AlbersEqualArea(\n    central_latitude=radar.latitude[\"data\"][0],\n    central_longitude=radar.longitude[\"data\"][0],\n)\n\nfaint_strong_cmap = mcolors.LinearSegmentedColormap.from_list(\n    \"faint_strong\",\n    [\n        (32 / 255, 144 / 255, 140 / 255),\n        (254 / 255, 238 / 255, 93 / 255),\n        (254 / 255, 152 / 255, 93 / 255),\n    ],\n    N=3,\n)\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    \"features_faint_strong\",\n    vmin=1,\n    vmax=3,\n    cmap=faint_strong_cmap,\n    axislabels_flag=False,\n    transform=ccrs.PlateCarree(),\n    ticks=[1.33, 2, 2.66],\n    ticklabs=[\"Background\", \"Strong Feature\", \"Faint Feature\"],\n    ax=ax2,\n)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Image muted Strong and Faint features**\n\nThis last section shows how we can image mute the fields (Tomkins et al. 2022) and reduce the visual prominence of\nmelting and mixed precipitation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# create a function to quickly apply image muting to other fields\ndef quick_image_mute(field, muted_ref):\n    nonmuted_field = np.ma.masked_where(~muted_ref.mask, field)\n    muted_field = np.ma.masked_where(muted_ref.mask, field)\n\n    return nonmuted_field, muted_field\n\n\n# get fields\nsnow_rate = grid.fields[\"snow_rate\"][\"data\"][0, :, :]\nmuted_ref = grid.fields[\"muted_reflectivity\"][\"data\"][0, :, :]\nfeatures = grid.fields[\"features_faint_strong\"][\"data\"][0, :, :]\nx = grid.x[\"data\"]\ny = grid.y[\"data\"]\n\n# mute\nsnow_rate_nonmuted, snow_rate_muted = quick_image_mute(snow_rate, muted_ref)\nfeatures_nonmuted, features_muted = quick_image_mute(features, muted_ref)\n\n# muted colormap\ngrays_cmap = plt.get_cmap(\"gray_r\")\nmuted_cmap = mcolors.LinearSegmentedColormap.from_list(\n    \"muted_cmap\", grays_cmap(np.linspace(0.2, 0.8, grays_cmap.N)), 3\n)\n\n\n# plot\nfig, axs = plt.subplots(1, 2, figsize=(10, 4))\n# snow rate\nsrpm = axs[0].pcolormesh(x / 1000, y / 1000, snow_rate_nonmuted, vmin=0, vmax=10)\nsrpmm = axs[0].pcolormesh(\n    x / 1000, y / 1000, snow_rate_muted, vmin=0, vmax=10, cmap=muted_cmap\n)\naxs[0].set_aspect(\"equal\")\naxs[0].set_title(\"Snow rate [mm hr$^{-1}$]\")\nplt.colorbar(srpm, ax=axs[0])\n\n# features\ncsbpm = axs[1].pcolormesh(\n    x / 1000, y / 1000, features_nonmuted, vmin=1, vmax=3, cmap=faint_strong_cmap\n)\ncsbpmm = axs[1].pcolormesh(\n    x / 1000, y / 1000, features_muted, vmin=1, vmax=3, cmap=muted_cmap\n)\n\naxs[1].set_aspect(\"equal\")\naxs[1].set_title(\"Feature detection\")\ncb = plt.colorbar(csbpm, ax=axs[1], ticks=[1.33, 2, 2.66])\ncb.ax.set_yticklabels([\"Background\", \"Strong Feature\", \"Faint Feature\"])\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\nTomkins, L. M., Yuter, S. E., Miller, M. A., and Allen, L. R., 2022:\nImage muting of mixed precipitation to improve identification of regions\nof heavy snow in radar data. Atmos. Meas. Tech., 15, 5515\u20135525,\nhttps://doi.org/10.5194/amt-15-5515-2022\n\nTomkins, L. M., Yuter, S. E., and Miller, M. A., 2024: Dual adaptive differential\nthreshold method for automated detection of faint and strong echo features\nin radar observations of winter storms. Atmos. Meas. Tech., 17, 3377\u20133399,\nhttps://doi.org/10.5194/amt-17-3377-2024\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.14.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}