{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Creating a CFAD diagram\nThis example shows how to create a contoured frequency by altitude (CFAD) diagram\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 matplotlib.pyplot as plt\nimport numpy as np\nfrom open_radar_data import DATASETS\n\nimport pyart"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Description of a CFAD\n\nA contoured frequency by altitude diagram (CFAD) is essentially a 2D histogram used to depict the vertical\ndistribution of a particular variable, such as radar reflectivity. The x-axis represents the frequency of the\nfield and the y-axis represents the frequency of the altitude. A key feature that distinguishes a CFAD from a\nregular 2D histogram is that it is normalized by altitude and altitudes where there is insufficient data are\nremoved. See Yuter and Houze (1995) for full details of the diagram.\n\n**Interpretation**\nIn a CFAD diagram, for a given altitude, the normalized frequency values will sum to a value of 1. When\ninterpreting a CFAD diagram, the normalized frequency for a given bin will be the frequency of that value for the\nassociated altitude (i.e. each bin can be interpreted as the fraction of data points at each altitude).\n\n**Minimum fraction threshold**\nIn a CFAD diagram, altitudes where there is insufficient data are removed. This feature is controlled with the\n`min_frac_thres` value. The default value is 0.1, meaning that altitudes where the number of observations is less\nthan one tenth of the maximum number of observations across all altitudes are removed. Increasing the value will\nact to be more aggressive in removing altitudes and decreasing the value will act to include more data. Setting\nthis value to zero will include all altitudes with data, but use caution when interpreting CFADs with all data\navailable as those altitudes with less data will not be as representative as other altitudes with sufficient data.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Basic example with reflectivity from RHI\n\n**Example with RHI**\nFirst, we will show an example from an RHI scan. We will get the reflectivity data and mask outside -15 and 35 dBZ\nto remove any noisy data. For the best results, we recommend cleaning up the field as much as possible (i.e.\ndespeckling, filtering, etc.)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get test data\nfilename = pyart.testing.get_test_data(\"sgpxsaprrhicmacI5.c0.20110524.015604_NC4.nc\")\nradar = pyart.io.read_cfradial(filename)\n\n# compute CFAD\n# get reflectivity data and mask extremes\n# extract first sweep\nradar = radar.extract_sweeps([0])\n# get mask array\nref_data = radar.fields[\"reflectivity_horizontal\"][\"data\"][:]\nref_data_masked = np.ma.masked_outside(ref_data, -15, 35)\nfield_mask = ref_data_masked.mask\n\n# get CFAD\nfreq_norm, height_edges, field_edges = pyart.retrieve.create_cfad(\n    radar,\n    field_bins=np.linspace(-15, 35, 100),\n    altitude_bins=np.linspace(0, 15000, 50),\n    field=\"reflectivity_horizontal\",\n    field_mask=field_mask,\n    min_frac_thres=0.1,\n)\n\n# plot CFAD\n# mask frequency values less than zero\nfreq_norm_masked = np.ma.masked_less_equal(freq_norm, 0)\n# plot CFAD\nplt.figure()\nax = plt.axes()\ncfad_pm = ax.pcolormesh(\n    field_edges, height_edges / 1000, freq_norm_masked, cmap=\"plasma\", vmin=0, vmax=0.10\n)\nplt.colorbar(cfad_pm)\nax.set_xlabel(\"Reflectivity [dBZ]\")\nax.set_ylabel(\"Height [km]\")\nax.grid(ls=\"--\", color=\"gray\", lw=0.5, alpha=0.7)\nplt.show()\n\n# plot RHI data\ndisplay = pyart.graph.RadarDisplay(radar)\n\nplt.figure(figsize=[12, 3])\nax = plt.axes()\nplt.tight_layout()\ndisplay.plot(\n    \"reflectivity_horizontal\",\n    0,\n    vmin=-15,\n    vmax=35,\n    mask_outside=True,\n    cmap=\"magma_r\",\n    ax=ax,\n)\ndisplay.set_limits(ylim=[0, 15], ax=ax)\nax.set_aspect(\"equal\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can see a general increase in reflectivity values from the echo top to around 8km. The maximum frequency value in\neach altitude row represents the mode of the reflectivity distribution which we can see also increases from echo\ntop to 8km. Below 8km the distribution of reflectivity values widens, likely associated with some of the noise in\nthe RHI.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Minimum fraction threshold example\n# ----------------------------------\n# Previously, we used the default `min_frac_thres` of 0.1. Next, we will increase the threshold and set the threshold\n# to 0 (not recommended) so show how the CFAD changes.\n\n# Let's see the effect of changing the minimum fraction threshold:\nfreq_norm2, height_edges, field_edges = pyart.retrieve.create_cfad(\n    radar,\n    field_bins=np.linspace(-15, 35, 100),\n    altitude_bins=np.linspace(0, 15000, 50),\n    field=\"reflectivity_horizontal\",\n    field_mask=field_mask,\n    min_frac_thres=0.2,\n)\nfreq_norm0, height_edges, field_edges = pyart.retrieve.create_cfad(\n    radar,\n    field_bins=np.linspace(-15, 35, 100),\n    altitude_bins=np.linspace(0, 15000, 50),\n    field=\"reflectivity_horizontal\",\n    field_mask=field_mask,\n    min_frac_thres=0,\n)\n\n# plot CFAD\n# mask zero values\nfreq_norm2_masked = np.ma.masked_less_equal(freq_norm2, 0)\nfreq_norm0_masked = np.ma.masked_less_equal(freq_norm0, 0)\n# plot\nplt.figure(figsize=(12, 3))\nax1 = plt.subplot(1, 3, 3)\ncfad_pm = ax1.pcolormesh(\n    field_edges,\n    height_edges / 1000,\n    freq_norm2_masked,\n    cmap=\"plasma\",\n    vmin=0,\n    vmax=0.10,\n)\nplt.colorbar(cfad_pm, ax=ax1)\nax1.set_xlabel(\"Reflectivity [dBZ]\")\nax1.grid(ls=\"--\", color=\"gray\", lw=0.5, alpha=0.7)\nax1.set_title(\"min_frac_thres = 0.2\")\n\nax2 = plt.subplot(1, 3, 2)\ncfad_pm = ax2.pcolormesh(\n    field_edges,\n    height_edges / 1000,\n    freq_norm_masked,\n    cmap=\"plasma\",\n    vmin=0,\n    vmax=0.10,\n)\nplt.colorbar(cfad_pm, ax=ax2)\nax2.set_xlabel(\"Reflectivity [dBZ]\")\nax2.grid(ls=\"--\", color=\"gray\", lw=0.5, alpha=0.7)\nax2.set_title(\"min_frac_thres = 0.1\")\n\nax3 = plt.subplot(1, 3, 1)\ncfad_pm = ax3.pcolormesh(\n    field_edges,\n    height_edges / 1000,\n    freq_norm0_masked,\n    cmap=\"plasma\",\n    vmin=0,\n    vmax=0.10,\n)\nplt.colorbar(cfad_pm, ax=ax3)\nax3.set_xlabel(\"Reflectivity [dBZ]\")\nax3.set_ylabel(\"Height [km]\")\nax3.grid(ls=\"--\", color=\"gray\", lw=0.5, alpha=0.7)\nax3.set_title(\"min_frac_thres = 0\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Setting the `min_frac_thres` to 0 (left panel) shows all data, even near the top of the chart (14km) where there is\nlimited echo. Setting the `min_frac_thres` higher to 0.2 (right panel) removed altitudes between 1 and 5 km where\nthere is less echo than between 6 and 12km where there is a consistent swath of reflectivity throughout the entire\ncross section.\n\n\n## Velocity example\nNext, we will show a CFAD for the doppler velocity from the above example. First, we have to dealias the velocity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# create a gatefilter\ngatefilter = pyart.filters.GateFilter(radar)\ngatefilter.exclude_invalid(\"reflectivity_horizontal\")\ngatefilter.exclude_outside(\"reflectivity_horizontal\", -15, 30)\ngatefilter.exclude_invalid(\"mean_doppler_velocity\")\n\n# dealias velocity\nvelocity_dealias = pyart.correct.dealias_region_based(\n    radar,\n    gatefilter=gatefilter,\n    vel_field=\"mean_doppler_velocity\",\n    nyquist_vel=13,\n)\nradar.add_field(\"corrected_velocity\", velocity_dealias)\n\n# plot RHI data\ndisplay = pyart.graph.RadarDisplay(radar)\n\nplt.figure(figsize=[12, 3])\nax = plt.axes()\nplt.tight_layout()\ndisplay.plot(\n    \"corrected_velocity\",\n    0,\n    vmin=-20,\n    vmax=20,\n    mask_outside=True,\n    cmap=\"RdBu_r\",\n    ax=ax,\n)\ndisplay.set_limits(ylim=[0, 15], ax=ax)\nax.set_aspect(\"equal\")\nplt.show()\n\nfreq_norm, height_edges, field_edges = pyart.retrieve.create_cfad(\n    radar,\n    field_bins=np.linspace(-30, 30, 50),\n    altitude_bins=np.linspace(0, 15000, 50),\n    field=\"corrected_velocity\",\n    field_mask=field_mask,\n    min_frac_thres=0.2,\n)\n\n# plot CFAD\n# mask zero values\nfreq_norm_masked = np.ma.masked_less_equal(freq_norm, 0)\n\n# plot\nplt.figure()\nax = plt.axes()\ncfad_pm = ax.pcolormesh(\n    field_edges,\n    height_edges / 1000,\n    freq_norm_masked,\n    cmap=\"plasma\",\n    vmin=0,\n    vmax=0.20,\n)\nplt.colorbar(cfad_pm)\nax.set_xlabel(\"Velocity [m s$^{-1}$]\")\nax.set_ylabel(\"Height [km]\")\nax.grid(ls=\"--\", color=\"gray\", lw=0.5, alpha=0.7)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The velocity CFAD is very different from the reflectivity CFAD. In most altitudes, there is more of a bimodal\npattern associated with the changing sign of the velocity values on either side of the radar. In general,\nthe distribution of velocity values is consistently wide throughout the profile compared to the reflectivity CFAD.\n\n## Validation\nFinally, we wanted to compare this function with the original method, so here we reproduce Fig. 2c from Yuter and\nHouze (1995) to demonstrate that it works the same. Instead of using the `pcolormesh` function, we are using\ncontour lines.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get test data\nfilename = DATASETS.fetch(\"ddop.910815.213931.cdf\")\ngrid = pyart.io.read_grid(filename)\n\n# make corrections to altitude field\naltitude_data = grid.point_z[\"data\"]\ngrid.point_z[\"data\"] = (altitude_data - 800) / 1000\n\n# get fields to create a mask\nfield_data = grid.fields[\"maxdz\"][\"data\"][:]\nvvel_data = grid.fields[\"w_wind\"][\"data\"][:]\nvvel_masked = np.ma.masked_invalid(vvel_data)\nfield_data_masked = np.ma.masked_less_equal(field_data, -15)\nfield_data_masked = np.ma.masked_where(vvel_masked.mask, field_data_masked)\n\n# define histogram bins\nfield_bins = np.arange(-20, 65, 5)\naltitude_bins = np.arange(-0.2, 18.5, 0.4)\n\nfreq_norm, height_edges, field_edges = pyart.retrieve.create_cfad(\n    grid,\n    field_bins=field_bins,\n    altitude_bins=altitude_bins,\n    field=\"maxdz\",\n    field_mask=field_data_masked.mask,\n    min_frac_thres=0.1,\n)\n\n# plot CFAD with contour plot\nfreq_norm_masked = np.ma.masked_less_equal(freq_norm, 0)\nh, f = np.meshgrid(height_edges, field_edges)\n\nplt.figure(figsize=(6, 3))\nax = plt.axes()\ncont = ax.contour(\n    f[:-1, :-1] + 2.5,\n    h[:-1, :-1] + 0.2,\n    freq_norm.T,\n    levels=np.arange(0.05, 0.3, 0.05),\n    colors=\"black\",\n)\nax.set_yticks([0, 5, 10, 15])\nax.set_xticks([-10, 0, 10, 20, 30, 40, 50])\nax.set_xlim([-12, 58])\nax.set_ylim([-0.5, 16])\nax.set_ylabel(\"Height [km]\")\nax.set_xlabel(\"Reflectivity [dBZ]\")\nax.axhline(8, ls=\"--\", lw=0.75, color=\"black\")\nax.set_title(\"Yuter and Houze (1995) Fig. 2c\")\nplt.show()\n\n# References\n# ----------\n# Yuter, S. E., and R. A. Houze, 1995: Three-Dimensional Kinematic and Microphysical Evolution of Florida\n# Cumulonimbus. Part II: Frequency Distributions of Vertical Velocity, Reflectivity, and Differential Reflectivity.\n# Mon. Wea. Rev. 123, 1941-1963. https://doi.org/10.1175/1520-0493(1995)123%3C1941:TDKAME%3E2.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
}