From 35aaff2fb5837d8210f0080688753d820f52d3be Mon Sep 17 00:00:00 2001 From: lauratomkins Date: Tue, 14 Nov 2023 14:05:19 -0500 Subject: [PATCH] ENH: updated feature detection function and example (#1487) * Rename functions and variables to remove "convective-stratiform" We want to make this function more generic so it can be used on a wide variety of fields and not just to detect convective and stratiform elements. I renamed most of the functions and some of the variables to reflect this and hopefully make it easier to use. * Add new reference * Incorporate binary closing option to functions * Add author to functions * Slight reformat for linting * More linting formatting * Add alternate option for getting under/over estimate field * Create plot_feature_detection.py * Delete plot_convective_stratiform.py * Update echo_class.py * Update plot_feature_detection.py * Update plot_feature_detection.py * Update test_echo_class.py * Keep original function and add warning * Revert "Delete plot_convective_stratiform.py" This reverts commit 2bff81ac0d8c33b3a51c4d6485d43994d52956ad. * Update __init__.py * Update echo_class.py * Update dictionary keys * Update plot_convective_stratiform.py * Fix issues * Update plot_convective_stratiform.py --- .../retrieve/plot_convective_stratiform.py | 42 +- examples/retrieve/plot_feature_detection.py | 810 ++++++++++++++++++ pyart/retrieve/__init__.py | 1 + pyart/retrieve/_echo_class.py | 457 +++++----- pyart/retrieve/echo_class.py | 316 +++++-- pyart/util/radar_utils.py | 2 + tests/retrieve/test_echo_class.py | 22 +- 7 files changed, 1328 insertions(+), 322 deletions(-) create mode 100644 examples/retrieve/plot_feature_detection.py diff --git a/examples/retrieve/plot_convective_stratiform.py b/examples/retrieve/plot_convective_stratiform.py index da0f22cc9f..2d905ee208 100644 --- a/examples/retrieve/plot_convective_stratiform.py +++ b/examples/retrieve/plot_convective_stratiform.py @@ -142,13 +142,13 @@ # add to grid object # mask zero values (no surface echo) -convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0) +convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0) # mask 3 values (weak echo) convsf_masked = np.ma.masked_equal(convsf_masked, 3) # add dimension to array to add to grid object -convsf_dict["convsf"]["data"] = convsf_masked[None, :, :] +convsf_dict["feature_detection"]["data"] = convsf_masked # add field -grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True) +grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True) # create plot using GridMapDisplay # plot variables @@ -195,23 +195,23 @@ # ``estimate_flag=False``), but we recommend keeping it turned on. # mask weak echo and no surface echo -convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0) +convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0) convsf_masked = np.ma.masked_equal(convsf_masked, 3) -convsf_dict["convsf"]["data"] = convsf_masked +convsf_dict["feature_detection"]["data"] = convsf_masked # underest. -convsf_masked = np.ma.masked_equal(convsf_dict["convsf_under"]["data"], 0) +convsf_masked = np.ma.masked_equal(convsf_dict["feature_under"]["data"], 0) convsf_masked = np.ma.masked_equal(convsf_masked, 3) -convsf_dict["convsf_under"]["data"] = convsf_masked +convsf_dict["feature_under"]["data"] = convsf_masked # overest. -convsf_masked = np.ma.masked_equal(convsf_dict["convsf_over"]["data"], 0) +convsf_masked = np.ma.masked_equal(convsf_dict["feature_over"]["data"], 0) convsf_masked = np.ma.masked_equal(convsf_masked, 3) -convsf_dict["convsf_over"]["data"] = convsf_masked +convsf_dict["feature_over"]["data"] = convsf_masked # Plot each estimation plt.figure(figsize=(10, 4)) ax1 = plt.subplot(131) ax1.pcolormesh( - convsf_dict["convsf"]["data"][0, :, :], + convsf_dict["feature_detection"]["data"][0, :, :], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3), @@ -220,13 +220,19 @@ ax1.set_aspect("equal") ax2 = plt.subplot(132) ax2.pcolormesh( - convsf_dict["convsf_under"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3) + convsf_dict["feature_under"]["data"][0, :, :], + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), ) ax2.set_title("Underestimate") ax2.set_aspect("equal") ax3 = plt.subplot(133) ax3.pcolormesh( - convsf_dict["convsf_over"]["data"], vmin=0, vmax=2, cmap=plt.get_cmap("viridis", 3) + convsf_dict["feature_over"]["data"][0, :, :], + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), ) ax3.set_title("Overestimate") ax3.set_aspect("equal") @@ -274,13 +280,13 @@ # add to grid object # mask zero values (no surface echo) -convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0) +convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0) # mask 3 values (weak echo) convsf_masked = np.ma.masked_equal(convsf_masked, 3) # add dimension to array to add to grid object -convsf_dict["convsf"]["data"] = convsf_masked[None, :, :] +convsf_dict["feature_detection"]["data"] = convsf_masked # add field -grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True) +grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True) # create plot using GridMapDisplay # plot variables @@ -393,13 +399,13 @@ # add to grid object # mask zero values (no surface echo) -convsf_masked = np.ma.masked_equal(convsf_dict["convsf"]["data"], 0) +convsf_masked = np.ma.masked_equal(convsf_dict["feature_detection"]["data"], 0) # mask 3 values (weak echo) convsf_masked = np.ma.masked_equal(convsf_masked, 3) # add dimension to array to add to grid object -convsf_dict["convsf"]["data"] = convsf_masked[None, :, :] +convsf_dict["feature_detection"]["data"] = convsf_masked # add field -grid.add_field("convsf", convsf_dict["convsf"], replace_existing=True) +grid.add_field("convsf", convsf_dict["feature_detection"], replace_existing=True) # create plot using GridMapDisplay # plot variables diff --git a/examples/retrieve/plot_feature_detection.py b/examples/retrieve/plot_feature_detection.py new file mode 100644 index 0000000000..5122ce6928 --- /dev/null +++ b/examples/retrieve/plot_feature_detection.py @@ -0,0 +1,810 @@ +""" +======================================= +Feature Detection classification +======================================= +This example shows how to use the feature detection algorithm. We show how to use the algorithm to detect convective +and stratiform regions in warm season events and how the algorithm can be used to detect both weak and strong +features in cool-season events. +""" + +print(__doc__) + +# Author: Laura Tomkins (lmtomkin@ncsu.edu) +# License: BSD 3 clause + + +import cartopy.crs as ccrs +import matplotlib.colors as mcolors +import matplotlib.pyplot as plt +import numpy as np +from open_radar_data import DATASETS + +import pyart + +###################################### +# How the algorithm works +# ---------- +# The feature detection algorithm works by identifying features that exceed the background value by an amount that +# varies with the background value. The algorithm is heavily customizable and is designed to work with a variety of +# datasets. Here, we show several examples of how to use the algorithm with different types of radar data. +# +# See Steiner et al. (1995), Yuter and Houze (1997), and Yuter et al. (2005) for full details on the algorithm. A +# manuscript (Tomkins et al. 2024) is in prep to describe feature detection in cool-season events (part 2). + +###################################### +# Part 1: Warm-season convective-stratiform classification +# ---------- +# **Classification of summer convective example** +# +# Our first example classifies echo from a summer convective event. + +# read in file +filename = pyart.testing.get_test_data("swx_20120520_0641.nc") +radar = pyart.io.read(filename) + +# extract the lowest sweep +radar = radar.extract_sweeps([0]) + +# interpolate to grid +grid = pyart.map.grid_from_radars( + (radar,), + grid_shape=(1, 201, 201), + grid_limits=((0, 10000), (-50000.0, 50000.0), (-50000.0, 50000.0)), + fields=["reflectivity_horizontal"], +) + +# get dx dy +dx = grid.x["data"][1] - grid.x["data"][0] +dy = grid.y["data"][1] - grid.y["data"][0] + +# feature detection +feature_dict = pyart.retrieve.feature_detection( + grid, + dx, + dy, + field="reflectivity_horizontal", + always_core_thres=40, + bkg_rad_km=20, + use_cosine=True, + max_diff=5, + zero_diff_cos_val=55, + weak_echo_thres=10, + max_rad_km=2, +) + +# add to grid object +# mask zero values (no surface echo) +feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0) +# mask 3 values (weak echo) +feature_masked = np.ma.masked_equal(feature_masked, 3) +# add dimension to array to add to grid object +feature_dict["feature_detection"]["data"] = feature_masked +# add field +grid.add_field( + "feature_detection", feature_dict["feature_detection"], replace_existing=True +) + +# create plot using GridMapDisplay +# plot variables +display = pyart.graph.GridMapDisplay(grid) +magma_r_cmap = plt.get_cmap("magma_r") +ref_cmap = mcolors.LinearSegmentedColormap.from_list( + "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N)) +) +projection = ccrs.AlbersEqualArea( + central_latitude=radar.latitude["data"][0], + central_longitude=radar.longitude["data"][0], +) + +# plot +plt.figure(figsize=(10, 4)) +ax1 = plt.subplot(1, 2, 1, projection=projection) +display.plot_grid( + "reflectivity_horizontal", + vmin=5, + vmax=45, + cmap=ref_cmap, + transform=ccrs.PlateCarree(), + ax=ax1, +) +ax2 = plt.subplot(1, 2, 2, projection=projection) +display.plot_grid( + "feature_detection", + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), + ax=ax2, + transform=ccrs.PlateCarree(), + ticks=[1 / 3, 1, 5 / 3], + ticklabs=["", "Stratiform", "Convective"], +) +plt.show() + +###################################### +# In addition to the default feature detection classification, the function also returns an underestimate +# (``feature_under``) and an overestimate (``feature_over``) to take into consideration the uncertainty when choosing +# classification parameters. The under and overestimate use the same parameters, but vary the input field by a +# certain value (default is 5 dBZ, can be changed with ``estimate_offset``). The estimation can be turned off ( +# ``estimate_flag=False``), but we recommend keeping it turned on. + +# mask weak echo and no surface echo +feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0) +feature_masked = np.ma.masked_equal(feature_masked, 3) +feature_dict["feature_detection"]["data"] = feature_masked +# underest. +feature_masked = np.ma.masked_equal(feature_dict["feature_under"]["data"], 0) +feature_masked = np.ma.masked_equal(feature_masked, 3) +feature_dict["feature_under"]["data"] = feature_masked +# overest. +feature_masked = np.ma.masked_equal(feature_dict["feature_over"]["data"], 0) +feature_masked = np.ma.masked_equal(feature_masked, 3) +feature_dict["feature_over"]["data"] = feature_masked + +# Plot each estimation +plt.figure(figsize=(10, 4)) +ax1 = plt.subplot(131) +ax1.pcolormesh( + feature_dict["feature_detection"]["data"][0, :, :], + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), +) +ax1.set_title("Best estimate") +ax1.set_aspect("equal") +ax2 = plt.subplot(132) +ax2.pcolormesh( + feature_dict["feature_under"]["data"][0, :, :], + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), +) +ax2.set_title("Underestimate") +ax2.set_aspect("equal") +ax3 = plt.subplot(133) +ax3.pcolormesh( + feature_dict["feature_over"]["data"][0, :, :], + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), +) +ax3.set_title("Overestimate") +ax3.set_aspect("equal") +plt.show() + +###################################### +# **Tropical example** +# +# Let's get a NEXRAD file from Hurricane Ian + +# Read in file +nexrad_file = "s3://noaa-nexrad-level2/2022/09/28/KTBW/KTBW20220928_190142_V06" +radar = pyart.io.read_nexrad_archive(nexrad_file) + +# extract the lowest sweep +radar = radar.extract_sweeps([0]) + +# interpolate to grid +grid = pyart.map.grid_from_radars( + (radar,), + grid_shape=(1, 201, 201), + grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)), + fields=["reflectivity"], +) + +# get dx dy +dx = grid.x["data"][1] - grid.x["data"][0] +dy = grid.y["data"][1] - grid.y["data"][0] + +# feature detection +feature_dict = pyart.retrieve.feature_detection( + grid, + dx, + dy, + field="reflectivity", + always_core_thres=40, + bkg_rad_km=20, + use_cosine=True, + max_diff=3, + zero_diff_cos_val=55, + weak_echo_thres=5, + max_rad_km=2, + estimate_flag=False, +) + +# add to grid object +# mask zero values (no surface echo) +feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0) +# mask 3 values (weak echo) +feature_masked = np.ma.masked_equal(feature_masked, 3) +# add dimension to array to add to grid object +feature_dict["feature_detection"]["data"] = feature_masked +# add field +grid.add_field( + "feature_detection", feature_dict["feature_detection"], replace_existing=True +) + +# create plot using GridMapDisplay +# plot variables +display = pyart.graph.GridMapDisplay(grid) +magma_r_cmap = plt.get_cmap("magma_r") +ref_cmap = mcolors.LinearSegmentedColormap.from_list( + "ref_cmap", magma_r_cmap(np.linspace(0, 0.9, magma_r_cmap.N)) +) +projection = ccrs.AlbersEqualArea( + central_latitude=radar.latitude["data"][0], + central_longitude=radar.longitude["data"][0], +) +# plot +plt.figure(figsize=(10, 4)) +ax1 = plt.subplot(1, 2, 1, projection=projection) +display.plot_grid( + "reflectivity", + vmin=5, + vmax=45, + cmap=ref_cmap, + transform=ccrs.PlateCarree(), + ax=ax1, + axislabels_flag=False, +) +ax2 = plt.subplot(1, 2, 2, projection=projection) +display.plot_grid( + "feature_detection", + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), + axislabels_flag=False, + transform=ccrs.PlateCarree(), + ticks=[1 / 3, 1, 5 / 3], + ticklabs=["", "Stratiform", "Convective"], + ax=ax2, +) +plt.show() + +###################################### +# **Comparison with original C++ algorithm** +# +# To ensure that our algorithm here matches the original algorithm developed in C++, we compare the results from an +# example during the KWAJEX field campaign. The file contains the original convective-stratiform classification from +# the C++ algorithm. We use the algorithm implemented in Py-ART with the same settings to produce identical results. + +# read in file +filename = DATASETS.fetch("convsf.19990811.221202.cdf") +grid = pyart.io.read_grid(filename) + +# colormap +convsf_cmap = mcolors.LinearSegmentedColormap.from_list( + "convsf", + [ + (33 / 255, 145 / 255, 140 / 255), + (253 / 255, 231 / 255, 37 / 255), + (210 / 255, 180 / 255, 140 / 255), + ], + N=3, +) + +# get original data from file processed in C++ +x = grid.x["data"] +y = grid.y["data"] +ref = grid.fields["maxdz"]["data"][0, :, :] +csb = grid.fields["convsf"]["data"][0, :, :] +csb_masked = np.ma.masked_less_equal(csb, 0) +csl = grid.fields["convsf_lo"]["data"][0, :, :] +csl_masked = np.ma.masked_less_equal(csl, 0) +csh = grid.fields["convsf_hi"]["data"][0, :, :] +csh_masked = np.ma.masked_less_equal(csh, 0) + +# plot +fig, axs = plt.subplots(2, 2, figsize=(12, 10)) +# reflectivity +rpm = axs[0, 0].pcolormesh(x / 1000, y / 1000, ref, vmin=0, vmax=45, cmap="magma_r") +axs[0, 0].set_aspect("equal") +axs[0, 0].set_title("Reflectivity [dBZ]") +plt.colorbar(rpm, ax=axs[0, 0]) +# convsf +csbpm = axs[0, 1].pcolormesh( + x / 1000, y / 1000, csb_masked, vmin=1, vmax=3, cmap=convsf_cmap +) +axs[0, 1].set_aspect("equal") +axs[0, 1].set_title("convsf") +cb = plt.colorbar(csbpm, ax=axs[0, 1], ticks=[4 / 3, 2, 8 / 3]) +cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"]) +# convsf lo +cslpm = axs[1, 0].pcolormesh( + x / 1000, y / 1000, csl_masked, vmin=1, vmax=3, cmap=convsf_cmap +) +axs[1, 0].set_aspect("equal") +axs[1, 0].set_title("convsf lo") +cb = plt.colorbar(cslpm, ax=axs[1, 0], ticks=[]) +# convsf hi +cshpm = axs[1, 1].pcolormesh( + x / 1000, y / 1000, csh_masked, vmin=1, vmax=3, cmap=convsf_cmap +) +axs[1, 1].set_aspect("equal") +axs[1, 1].set_title("convsf hi") +cb = plt.colorbar(cshpm, ax=axs[1, 1], ticks=[4 / 3, 2, 8 / 3]) +cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"]) +plt.show() + +# now let's compare with the Python algorithm +convsf_dict = pyart.retrieve.feature_detection( + grid, + dx=x[1] - x[0], + dy=y[1] - y[0], + level_m=None, + always_core_thres=40, + bkg_rad_km=11, + use_cosine=True, + max_diff=8, + zero_diff_cos_val=55, + scalar_diff=None, + use_addition=False, + calc_thres=0, + weak_echo_thres=15, + min_val_used=0, + dB_averaging=True, + remove_small_objects=False, + min_km2_size=None, + val_for_max_rad=30, + max_rad_km=5.0, + core_val=3, + nosfcecho=0, + weakecho=3, + bkgd_val=1, + feat_val=2, + field="maxdz", + estimate_flag=True, + estimate_offset=5, +) + +# new data +csb_lt = convsf_dict["feature_detection"]["data"][0, :, :] +csb_lt_masked = np.ma.masked_less_equal(csb_lt, 0) +csu_lt = convsf_dict["feature_under"]["data"][0, :, :] +csu_lt_masked = np.ma.masked_less_equal(csu_lt, 0) +cso_lt = convsf_dict["feature_over"]["data"][0, :, :] +cso_lt_masked = np.ma.masked_less_equal(cso_lt, 0) + +fig, axs = plt.subplots(2, 2, figsize=(12, 10)) +# reflectivity +rpm = axs[0, 0].pcolormesh(x / 1000, y / 1000, ref, vmin=0, vmax=45, cmap="magma_r") +axs[0, 0].set_aspect("equal") +axs[0, 0].set_title("Reflectivity [dBZ]") +plt.colorbar(rpm, ax=axs[0, 0]) +# convsf best +csbpm = axs[0, 1].pcolormesh( + x / 1000, y / 1000, csb_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap +) +axs[0, 1].set_aspect("equal") +axs[0, 1].set_title("convsf") +cb = plt.colorbar(csbpm, ax=axs[0, 1], ticks=[4 / 3, 2, 8 / 3]) +cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"]) +# convsf under +csupm = axs[1, 0].pcolormesh( + x / 1000, y / 1000, csu_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap +) +axs[1, 0].set_aspect("equal") +axs[1, 0].set_title("convsf under") +cb = plt.colorbar(csupm, ax=axs[1, 0], ticks=[]) +# convsf over +csopm = axs[1, 1].pcolormesh( + x / 1000, y / 1000, cso_lt_masked, vmin=1, vmax=3, cmap=convsf_cmap +) +axs[1, 1].set_aspect("equal") +axs[1, 1].set_title("convsf over") +cb = plt.colorbar(csopm, ax=axs[1, 1], ticks=[4 / 3, 2, 8 / 3]) +cb.ax.set_yticklabels(["Stratiform", "Convective", "Weak Echo"]) +plt.show() + +###################################### +# Part 2: Cool-season feature detection +# ---------- +###################################### +# **Winter storm example** +# +# In this example, we will show how to algorithm can be used to detect features (snow bands) in winter storms. Here, we +# rescale the reflectivity to snow rate (Rasumussen et al. 2003) in order to better represent the snow field. Note +# that we are not using the absolute values of the snow rate, we are only using the relative anomalies relative to +# the background average. We recommend using a rescaled reflectivity for this algorithm, but note that you need to +# change ``dB_averaging = False`` when using a rescaled field (set ``dB_averaging`` to True for reflectivity fields in +# dBZ units). Also note that since we are now using a snow rate field with different values, most of our parameters +# have changed as well. + +# Read in file +nexrad_file = "s3://noaa-nexrad-level2/2021/02/07/KOKX/KOKX20210207_161413_V06" +radar = pyart.io.read_nexrad_archive(nexrad_file) + +# extract the lowest sweep +radar = radar.extract_sweeps([0]) + +# interpolate to grid +grid = pyart.map.grid_from_radars( + (radar,), + grid_shape=(1, 201, 201), + grid_limits=((0, 10000), (-200000.0, 200000.0), (-200000.0, 200000.0)), + fields=["reflectivity", "cross_correlation_ratio"], +) + +# image mute grid object +grid = pyart.util.image_mute_radar( + grid, "reflectivity", "cross_correlation_ratio", 0.97, 20 +) + +# rescale reflectivity to snow rate +grid = pyart.retrieve.qpe.ZtoR( + grid, ref_field="reflectivity", a=57.3, b=1.67, save_name="snow_rate" +) + +# get dx dy +dx = grid.x["data"][1] - grid.x["data"][0] +dy = grid.y["data"][1] - grid.y["data"][0] + +# feature detection +feature_dict = pyart.retrieve.feature_detection( + grid, + dx, + dy, + field="snow_rate", + dB_averaging=False, + always_core_thres=4, + bkg_rad_km=40, + use_cosine=True, + max_diff=1.5, + zero_diff_cos_val=5, + weak_echo_thres=0, + min_val_used=0, + max_rad_km=1, + estimate_flag=False, +) + +# add to grid object +# mask zero values (no surface echo) +feature_masked = np.ma.masked_equal(feature_dict["feature_detection"]["data"], 0) +# mask 3 values (weak echo) +feature_masked = np.ma.masked_equal(feature_masked, 3) +# add dimension to array to add to grid object +feature_dict["feature_detection"]["data"] = feature_masked +# add field +grid.add_field( + "feature_detection", feature_dict["feature_detection"], replace_existing=True +) + +# create plot using GridMapDisplay +# plot variables +display = pyart.graph.GridMapDisplay(grid) + +projection = ccrs.AlbersEqualArea( + central_latitude=radar.latitude["data"][0], + central_longitude=radar.longitude["data"][0], +) +# plot +plt.figure(figsize=(10, 4)) +ax1 = plt.subplot(1, 2, 1, projection=projection) +display.plot_grid( + "snow_rate", + vmin=0, + vmax=10, + cmap=plt.get_cmap("viridis"), + transform=ccrs.PlateCarree(), + ax=ax1, + axislabels_flag=False, +) +ax2 = plt.subplot(1, 2, 2, projection=projection) +display.plot_grid( + "feature_detection", + vmin=0, + vmax=2, + cmap=plt.get_cmap("viridis", 3), + axislabels_flag=False, + transform=ccrs.PlateCarree(), + ticks=[1 / 3, 1, 5 / 3], + ticklabs=["", "Background", "Feature"], + ax=ax2, +) +plt.show() + +###################################### +# **Under and over estimate in snow** +# +# Here we show how to bound the snow rate with an under and over estimate + +# get reflectivity field and copy +ref_field = grid.fields["reflectivity"] +ref_field_over = ref_field.copy() +ref_field_under = ref_field.copy() + +# offset original field by +/- 2dB +ref_field_over["data"] = ref_field["data"] + 2 +ref_field_under["data"] = ref_field["data"] - 2 + +# adjust dictionary +ref_field_over["standard_name"] = ref_field["standard_name"] + "_overest" +ref_field_over["long_name"] = ref_field["long_name"] + " overestimate" +ref_field_under["standard_name"] = ref_field["standard_name"] + "_underest" +ref_field_under["long_name"] = ref_field["long_name"] + " underestimate" + +# add to grid object +grid.add_field("reflectivity_over", ref_field_over, replace_existing=True) +grid.add_field("reflectivity_under", ref_field_under, replace_existing=True) + +# convert to snow rate +grid = pyart.retrieve.qpe.ZtoR( + grid, ref_field="reflectivity_over", a=57.3, b=1.67, save_name="snow_rate_over" +) +grid = pyart.retrieve.qpe.ZtoR( + grid, ref_field="reflectivity_under", a=57.3, b=1.67, save_name="snow_rate_under" +) + +# now do feature detection with under and over estimate fields +feature_estimate_dict = pyart.retrieve.feature_detection( + grid, + dx, + dy, + field="snow_rate", + dB_averaging=False, + always_core_thres=4, + bkg_rad_km=40, + use_cosine=True, + max_diff=1.5, + zero_diff_cos_val=5, + weak_echo_thres=0, + min_val_used=0, + max_rad_km=1, + estimate_flag=True, + overest_field="snow_rate_over", + underest_field="snow_rate_under", +) + +# x and y data for plotting +x = grid.x["data"] +y = grid.y["data"] + +# now let's plot +fig, axs = plt.subplots(2, 2, figsize=(12, 10)) +# snow rate +rpm = axs[0, 0].pcolormesh( + x / 1000, + y / 1000, + grid.fields["snow_rate"]["data"][0, :, :], + vmin=0, + vmax=10, +) +axs[0, 0].set_aspect("equal") +axs[0, 0].set_title("Reflectivity [dBZ]") +plt.colorbar(rpm, ax=axs[0, 0]) +# features best +bpm = axs[0, 1].pcolormesh( + x / 1000, + y / 1000, + feature_estimate_dict["feature_detection"]["data"][0, :, :], + vmin=0, + vmax=3, + cmap=plt.get_cmap("viridis", 3), +) +axs[0, 1].set_aspect("equal") +axs[0, 1].set_title("Feature detection (best)") +cb = plt.colorbar(bpm, ax=axs[0, 1], ticks=[3 / 2, 5 / 2]) +cb.ax.set_yticklabels(["Background", "Feature"]) +# features underestimate +upm = axs[1, 0].pcolormesh( + x / 1000, + y / 1000, + feature_estimate_dict["feature_under"]["data"][0, :, :], + vmin=0, + vmax=3, + cmap=plt.get_cmap("viridis", 3), +) +axs[1, 0].set_aspect("equal") +axs[1, 0].set_title("Feature detection (underestimate)") +cb = plt.colorbar(upm, ax=axs[1, 0], ticks=[3 / 2, 5 / 2]) +cb.ax.set_yticklabels(["Background", "Feature"]) +# features overestimate +opm = axs[1, 1].pcolormesh( + x / 1000, + y / 1000, + feature_estimate_dict["feature_over"]["data"][0, :, :], + vmin=0, + vmax=3, + cmap=plt.get_cmap("viridis", 3), +) +axs[1, 1].set_aspect("equal") +axs[1, 1].set_title("Feature detection (overestimate)") +cb = plt.colorbar(opm, ax=axs[1, 1], ticks=[3 / 2, 5 / 2]) +cb.ax.set_yticklabels(["Background", "Feature"]) +plt.show() + +###################################### +# **Strong and Faint features** +# +# We have developed a new technique which runs the algorithm twice to find features that are very distinct from the +# background (strong features) and features that are not very distinct from the background (faint features). + +# run the algorithm again, but to find objects that are also faint +# feature detection +feature_dict = pyart.retrieve.feature_detection( + grid, + dx, + dy, + field="snow_rate", + dB_averaging=False, + always_core_thres=4, + bkg_rad_km=40, + use_cosine=False, + scalar_diff=1.5, + use_addition=False, + weak_echo_thres=0, + min_val_used=0, + max_rad_km=1, + estimate_flag=False, +) + +# separate strong vs. faint +# mask zero values (no surface echo) +scalar_features_masked = np.ma.masked_equal( + feature_dict["feature_detection"]["data"], 0 +) +# mask 3 values (weak echo) +scalar_features_masked = np.ma.masked_equal(scalar_features_masked, 3) +# get cosine features +cosine_features_masked = grid.fields["feature_detection"]["data"] +# isolate faint features only +faint_features = scalar_features_masked - cosine_features_masked +faint_features_masked = np.ma.masked_less_equal(faint_features, 0) +# add strong and faint features +features_faint_strong = 2 * faint_features_masked.filled( + 0 +) + cosine_features_masked.filled(0) +features_faint_strong_masked = np.ma.masked_where( + cosine_features_masked.mask, features_faint_strong +) + +# add to grid object +new_dict = { + "features_faint_strong": { + "data": features_faint_strong_masked, + "standard_name": "features_faint_strong", + "long name": "Faint and Strong Features", + "valid_min": 0, + "valid_max": 10, + "comment_1": "3: Faint features, 2: Strong features, 1: background", + } +} +# add field +grid.add_field( + "features_faint_strong", new_dict["features_faint_strong"], replace_existing=True +) + +# create plot using GridMapDisplay +# plot variables +display = pyart.graph.GridMapDisplay(grid) + +projection = ccrs.AlbersEqualArea( + central_latitude=radar.latitude["data"][0], + central_longitude=radar.longitude["data"][0], +) + +faint_strong_cmap = mcolors.LinearSegmentedColormap.from_list( + "faint_strong", + [ + (32 / 255, 144 / 255, 140 / 255), + (254 / 255, 238 / 255, 93 / 255), + (254 / 255, 152 / 255, 93 / 255), + ], + N=3, +) + +# plot +plt.figure(figsize=(10, 4)) +ax1 = plt.subplot(1, 2, 1, projection=projection) +display.plot_grid( + "snow_rate", + vmin=0, + vmax=10, + cmap=plt.get_cmap("viridis"), + transform=ccrs.PlateCarree(), + ax=ax1, + axislabels_flag=False, +) +ax2 = plt.subplot(1, 2, 2, projection=projection) +display.plot_grid( + "features_faint_strong", + vmin=1, + vmax=3, + cmap=faint_strong_cmap, + axislabels_flag=False, + transform=ccrs.PlateCarree(), + ticks=[1.33, 2, 2.66], + ticklabs=["Background", "Strong Feature", "Faint Feature"], + ax=ax2, +) +plt.show() + +###################################### +# **Image muted Strong and Faint features** +# +# This last section shows how we can image mute the fields (Tomkins et al. 2022) and reduce the visual prominence of +# melting and mixed precipitation. + + +# create a function to quickly apply image muting to other fields +def quick_image_mute(field, muted_ref): + nonmuted_field = np.ma.masked_where(~muted_ref.mask, field) + muted_field = np.ma.masked_where(muted_ref.mask, field) + + return nonmuted_field, muted_field + + +# get fields +snow_rate = grid.fields["snow_rate"]["data"][0, :, :] +muted_ref = grid.fields["muted_reflectivity"]["data"][0, :, :] +features = grid.fields["features_faint_strong"]["data"][0, :, :] +x = grid.x["data"] +y = grid.y["data"] + +# mute +snow_rate_nonmuted, snow_rate_muted = quick_image_mute(snow_rate, muted_ref) +features_nonmuted, features_muted = quick_image_mute(features, muted_ref) + +# muted colormap +grays_cmap = plt.get_cmap("gray_r") +muted_cmap = mcolors.LinearSegmentedColormap.from_list( + "muted_cmap", grays_cmap(np.linspace(0.2, 0.8, grays_cmap.N)), 3 +) + + +# plot +fig, axs = plt.subplots(1, 2, figsize=(10, 4)) +# snow rate +srpm = axs[0].pcolormesh(x / 1000, y / 1000, snow_rate_nonmuted, vmin=0, vmax=10) +srpmm = axs[0].pcolormesh( + x / 1000, y / 1000, snow_rate_muted, vmin=0, vmax=10, cmap=muted_cmap +) +axs[0].set_aspect("equal") +axs[0].set_title("Snow rate [mm hr$^{-1}$]") +plt.colorbar(srpm, ax=axs[0]) + +# features +csbpm = axs[1].pcolormesh( + x / 1000, y / 1000, features_nonmuted, vmin=1, vmax=3, cmap=faint_strong_cmap +) +csbpmm = axs[1].pcolormesh( + x / 1000, y / 1000, features_muted, vmin=1, vmax=3, cmap=muted_cmap +) + +axs[1].set_aspect("equal") +axs[1].set_title("Feature detection") +cb = plt.colorbar(csbpm, ax=axs[1], ticks=[1.33, 2, 2.66]) +cb.ax.set_yticklabels(["Background", "Strong Feature", "Faint Feature"]) + +plt.show() +###################################### +# Summary of recommendations and best practices +# ---------- +# * Tune your parameters to your specific purpose +# * Use a rescaled field if possible (i.e. linear reflectivity, rain or snow rate) +# * Keep ``estimate_flag=True`` to see uncertainty in classification +# +# References +# ---------- +# Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological +# Characterization of Three-Dimensional Storm Structure from Operational +# Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007. +# https://doi.org/10.1175/1520-0450(1995)034<1978:CCOTDS>2.0.CO;2. +# +# Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size +# distributions over the Pacific warm pool and implications for Z-R relations. +# J. Appl. Meteor., 36, 847-867. +# https://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2 +# +# Yuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser, +# 2005: Physical characterization of tropical oceanic convection observed in +# KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1 +# +# Rasmussen, R., M. Dixon, S. Vasiloff, F. Hage, S. Knight, J. Vivekanandan, +# and M. Xu, 2003: Snow Nowcasting Using a Real-Time Correlation of Radar +# Reflectivity with Snow Gauge Accumulation. J. Appl. Meteorol. Climatol., 42, 20–36. +# https://doi.org/10.1175/1520-0450(2003)042%3C0020:SNUART%3E2.0.CO;2 +# +# Tomkins, L. M., Yuter, S. E., Miller, M. A., and Allen, L. R., 2022: +# Image muting of mixed precipitation to improve identification of regions +# of heavy snow in radar data. Atmos. Meas. Tech., 15, 5515–5525, +# https://doi.org/10.5194/amt-15-5515-2022 diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py index 1276e35038..3b79ba5391 100644 --- a/pyart/retrieve/__init__.py +++ b/pyart/retrieve/__init__.py @@ -5,6 +5,7 @@ from .advection import grid_displacement_pc, grid_shift # noqa from .comp_z import composite_reflectivity # noqa +from .echo_class import feature_detection # noqa from .echo_class import conv_strat_yuter # noqa from .echo_class import get_freq_band # noqa from .echo_class import hydroclass_semisupervised # noqa diff --git a/pyart/retrieve/_echo_class.py b/pyart/retrieve/_echo_class.py index ca9ff8209a..362d351eca 100644 --- a/pyart/retrieve/_echo_class.py +++ b/pyart/retrieve/_echo_class.py @@ -239,8 +239,8 @@ def steiner_class_buff( return sclass -def _revised_conv_strat( - refl, +def _feature_detection( + field, dx, dy, always_core_thres=42, @@ -252,48 +252,48 @@ def _revised_conv_strat( use_addition=True, calc_thres=0.75, weak_echo_thres=5.0, - min_dBZ_used=5.0, + min_val_used=5.0, dB_averaging=True, remove_small_objects=True, min_km2_size=10, - val_for_max_conv_rad=30, - max_conv_rad_km=5.0, - cs_core=3, + binary_close=False, + val_for_max_rad=30, + max_rad_km=5.0, + core_val=3, nosfcecho=0, weakecho=3, - sf=1, - conv=2, + bkgd_val=1, + feat_val=2, ): """ - We perform the Yuter and Houze (1997) algorithm for echo classification - using only the reflectivity field in order to classify each grid point - as either convective, stratiform or undefined. Grid points are + These functions are used to detect features in fields based on how distinct they are from the background + average. Methodology described in Tomkins et al. (2023), originally based on convective-stratiform algorithms + developed by Steiner et al. (1995), Yuter and Houze (1997), and Yuter et al. (2005) Grid points are classified as follows, nosfcecho = No Surface Echo/ Undefined - sf = Stratiform - conv = Convective + bkgd_val = Background echo (e.g. Stratiform) + feat_val = Feature (e.g. Convective) weakecho = Weak Echo - refl : array - array of reflectivity values + field : array + array of values to find features x, y : array - x and y coordinates of reflectivity array, respectively + x and y coordinates of field array, respectively dx, dy : float The x- and y-dimension resolutions in meters, respectively. always_core_thres : float, optional - Threshold for points that are always convective. All values above the threshold are classifed as convective + Threshold for points that are always features. All values above the threshold are classified as features. bkg_rad_km : float, optional - Radius to compute background reflectivity in kilometers. Default is 11 km. Recommended to be at least 3 x - grid spacing + Radius to compute background field in kilometers. Default is 11 km. Recommended to be at least 3 x grid spacing use_cosine : bool, optional - Boolean used to determine if cosine scheme should be used for identifying convective cores (True) or a scalar - scheme (False) + Boolean used to determine if a cosine scheme (see Yuter and Houze (1997)) should be used for identifying + cores (True) or if a simpler scalar scheme (False) should be used. max_diff : float, optional - Maximum difference between background average and reflectivity in order to be classified as convective. + Maximum difference between background average and grid value in order to be classified as features. "a" value in Eqn. B1 in Yuter and Houze (1997) zero_diff_cos_val : float, optional - Value where difference between background average and reflectivity is zero in the cosine function + Value where difference between background average and grid value is zero in the cosine function "b" value in Eqn. B1 in Yuter and Houze (1997) scalar_diff : float, optional If using a scalar difference scheme, this value is the multiplier or addition to the background average @@ -303,50 +303,50 @@ def _revised_conv_strat( Minimum percentage of points needed to be considered in background average calculation weak_echo_thres : float, optional Threshold for determining weak echo. All values below this threshold will be considered weak echo - min_dBZ_used : float, optional - Minimum dBZ value used for classification. All values below this threshold will be considered no surface echo + min_val_used : float, optional + Minimum value used for classification. All values below this threshold will be considered no surface echo + See Yuter and Houze (1997) and Yuter et al. (2005) for more detail. Units based on input field dB_averaging : bool, optional True if using dBZ values that need to be converted to linear Z before averaging. False for other types of values remove_small_objects : bool, optional - Determines if small objects should be removed from convective core array. Default is True. + Determines if small objects should be removed from core array. Default is True. min_km2_size : float, optional - Minimum size of convective cores to be considered. Cores less than this size will be removed. Default is 10 - km^2. - val_for_max_conv_rad : float, optional - dBZ for maximum convective radius. Convective cores with values above this will have the maximum convective - radius - max_conv_rad_km : float, optional - Maximum radius around convective cores to classify as convective. Default is 5 km. - cs_core : int, optional - Value for points classified as convective cores + Minimum size of Cores to be considered. Cores less than this size will be removed. Default is 10 km^2. + binary_close : bool, optional + Determines if a binary closing should be performed on the cores. Default is False. + val_for_max_rad : float, optional + value used for maximum radius. Cores with values above this will have the maximum radius incorporated. + max_rad_km : float, optional + Maximum radius around cores to classify as feature. Default is 5 km + core_val : int, optional + Value for points classified as cores nosfcecho : int, optional - Value for points classified as no surface echo, based on min_dBZ_used + Value for points classified as no surface echo, based on min_val_used weakecho : int, optional - Value for points classified as weak echo, based on weak_echo_thres - sf : int, optional - Value for points classified as stratiform - conv : int, optional - Value for points classified as convective + Value for points classified as weak echo, based on weak_echo_thres. + bkgd_val : int, optional + Value for points classified as background echo. + feat_val : int, optional + Value for points classified as features. Returns ------- - refl_bkg : array + field_bkg : array Array of background values - conv_core_array : array - Array of initial convective cores (identified convective elements without convective radii applied) - conv_strat_array : array - Array of convective stratiform classifcation with convective radii applied + core_array : array + Array of initial cores (identified features without radii applied) + feature_array : array + Array of feature detection with radii applied """ - # Set up mask arrays for background average and - # prepare for convective mask arrays - # calculate maximum convective diameter from max. convective radius (input) - max_conv_diameter = int(np.floor((max_conv_rad_km / (dx / 1000)) * 2)) + # Set up mask arrays for background average and prepare for mask arrays + # calculate maximum diameter from max. radius (input) + max_diameter = int(np.floor((max_rad_km / (dx / 1000)) * 2)) # if diameter is even, make odd - if max_conv_diameter % 2 == 0: - max_conv_diameter = max_conv_diameter + 1 + if max_diameter % 2 == 0: + max_diameter = max_diameter + 1 # find center point - center_conv_mask_x = int(np.floor(max_conv_diameter / 2)) + center_mask_x = int(np.floor(max_diameter / 2)) # prepare background mask array for computing background average # calculate number of pixels for background array given requested background radius and dx @@ -370,42 +370,40 @@ def _revised_conv_strat( circular=True, ) - # Convective stratiform detection - # start by making reflectivity a masked array - refl = np.ma.masked_invalid(refl) + # Feature detection + # start by making field a masked array + field = np.ma.masked_invalid(field) # Compute background radius - refl_bkg = calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres) + field_bkg = calc_bkg_intensity(field, bkg_mask_array, dB_averaging, calc_thres) - # Get convective core array from cosine scheme, or scalar scheme + # Get core array from cosine scheme, or scalar scheme if use_cosine: - conv_core_array = convcore_cos_scheme( - refl, refl_bkg, max_diff, zero_diff_cos_val, always_core_thres, cs_core + core_array = core_cos_scheme( + field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, core_val ) else: - conv_core_array = convcore_scalar_scheme( - refl, - refl_bkg, + core_array = core_scalar_scheme( + field, + field_bkg, scalar_diff, always_core_thres, - cs_core, + core_val, use_addition=use_addition, ) - # Assign convective radii based on background reflectivity - conv_radius_km = assign_conv_radius_km( - refl_bkg, - val_for_max_conv_rad=val_for_max_conv_rad, - max_conv_rad=max_conv_rad_km, + # Assign radii based on background field + radius_array_km = assign_feature_radius_km( + field_bkg, val_for_max_rad=val_for_max_rad, max_rad=max_rad_km ) - # remove small objects in convective core array + # remove small objects in core array if remove_small_objects: # calculate minimum pixel size given dx and dy min_pix_size = min_km2_size / ((dx / 1000) * (dy / 1000)) - # label connected objects in convective core array - cc_labels, _ = scipy.ndimage.label(conv_core_array) - # mask labels where convective core array is masked - cc_labels = np.ma.masked_where(conv_core_array.mask, cc_labels) + # label connected objects in core array + cc_labels, _ = scipy.ndimage.label(core_array) + # mask labels where core array is masked + cc_labels = np.ma.masked_where(core_array.mask, cc_labels) # loop through each unique label for lab in np.unique(cc_labels): @@ -413,51 +411,58 @@ def _revised_conv_strat( size_lab = np.count_nonzero(cc_labels == lab) # if number of pixels is less than minimum, then remove core if size_lab < min_pix_size: - conv_core_array[cc_labels == lab] = 0 + core_array[cc_labels == lab] = 0 - # Incorporate convective radius using binary dilation + # perform binary closing + if binary_close: + # binary closing - returns binary array + close_core = scipy.ndimage.binary_closing(core_array).astype(int) + # set values to core values + core_array = close_core * core_val + + # Incorporate radius using binary dilation # Create empty array for assignment - temp_assignment = np.zeros_like(conv_core_array) + temp_assignment = np.zeros_like(core_array) # Loop through radii - for radius in np.arange(1, max_conv_rad_km + 1): + for radius in np.arange(1, max_rad_km + 1): # create mask array for radius incorporation - conv_mask_array = create_conv_radius_mask( - max_conv_diameter, radius, dx / 1000, dy / 1000, center_conv_mask_x + radius_mask_array = create_radius_mask( + max_diameter, radius, dx / 1000, dy / 1000, center_mask_x ) # find location of radius - temp = conv_radius_km == radius + temp = radius_array_km == radius # get cores for given radius - temp_core = np.ma.masked_where(~temp, conv_core_array) + temp_core = np.ma.masked_where(~temp, core_array) # dilate cores temp_dilated = scipy.ndimage.binary_dilation( - temp_core.filled(0), conv_mask_array + temp_core.filled(0), radius_mask_array ) # add to assignment array temp_assignment = temp_assignment + temp_dilated # add dilated cores to original array - conv_core_copy = np.ma.copy(conv_core_array) - conv_core_copy[temp_assignment >= 1] = cs_core - - # Now do convective stratiform classification - conv_strat_array = np.zeros_like(refl) - conv_strat_array = classify_conv_strat_array( - refl, - conv_strat_array, - conv_core_copy, + core_copy = np.ma.copy(core_array) + core_copy[temp_assignment >= 1] = core_val + + # Now do feature detection + feature_array = np.zeros_like(field) + feature_array = classify_feature_array( + field, + feature_array, + core_copy, nosfcecho, - conv, - sf, + feat_val, + bkgd_val, weakecho, - cs_core, - min_dBZ_used, + core_val, + min_val_used, weak_echo_thres, ) - # mask where reflectivity is masked - conv_strat_array = np.ma.masked_where(refl.mask, conv_strat_array) + # mask where field is masked + feature_array = np.ma.masked_where(field.mask, feature_array) - return refl_bkg, conv_core_array, conv_strat_array + return field_bkg, core_array, feature_array # functions @@ -524,15 +529,15 @@ def create_radial_mask( return mask_array -def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None): +def calc_bkg_intensity(field, bkg_mask_array, dB_averaging, calc_thres=None): """ - Calculate the background of the given refl array. The footprint used to + Calculate the background of the given field. The footprint used to calculate the average for each pixel is given by bkg_mask_array Parameters ---------- - refl : array - Reflectivity array to compute average + field : array + Field array to compute average bkg_mask_array : array Array of radial points to use for average dB_averaging : bool @@ -542,17 +547,17 @@ def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None): Returns ------- - refl_bkg : array + field_bkg : array Array of average values """ # if dBaverage is true, convert reflectivity to linear Z if dB_averaging: - refl = 10 ** (refl / 10) + field = 10 ** (field / 10) - # calculate background reflectivity with circular footprint - refl_bkg = scipy.ndimage.generic_filter( - refl.filled(np.nan), + # calculate background field with circular footprint + field_bkg = scipy.ndimage.generic_filter( + field.filled(np.nan), function=np.nanmean, mode="constant", footprint=bkg_mask_array.astype(bool), @@ -562,8 +567,8 @@ def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None): # if calc_thres is not none, then calculate the number of points used to calculate average if calc_thres is not None: # count valid points - refl_count = scipy.ndimage.generic_filter( - refl.filled(0), + field_count = scipy.ndimage.generic_filter( + field.filled(0), function=np.count_nonzero, mode="constant", footprint=bkg_mask_array.astype(bool), @@ -572,223 +577,221 @@ def calc_bkg_intensity(refl, bkg_mask_array, dB_averaging, calc_thres=None): # find threshold number of points val = calc_thres * np.count_nonzero(bkg_mask_array) # mask out values - refl_bkg = np.ma.masked_where(refl_count < val, refl_bkg) + field_bkg = np.ma.masked_where(field_count < val, field_bkg) - # mask where original reflectivity is invalid - refl_bkg = np.ma.masked_where(refl.mask, refl_bkg) + # mask where original field is invalid + field_bkg = np.ma.masked_where(field.mask, field_bkg) # if dBaveraging is true, convert background reflectivity to dBZ if dB_averaging: - refl_bkg = 10 * np.log10(refl_bkg) - # mask where original reflectivity is invalid - refl_bkg = np.ma.masked_where(refl.mask, refl_bkg) + field_bkg = 10 * np.log10(field_bkg) + # mask where original field is invalid + field_bkg = np.ma.masked_where(field.mask, field_bkg) - return refl_bkg + return field_bkg -def convcore_cos_scheme( - refl, refl_bkg, max_diff, zero_diff_cos_val, always_core_thres, CS_CORE +def core_cos_scheme( + field, field_bkg, max_diff, zero_diff_cos_val, always_core_thres, CS_CORE ): """ - Function for assigning convective cores based on a cosine function + Function for assigning cores based on a cosine function Parameters ---------- - refl : array - Reflectivity values - refl_bkg : array - Background average of reflectivity values + field : array + Field values + field_bkg : array + Background average of field values max_diff : float - Maximum difference between refl and refl_bkg needed for convective classification + Maximum difference between field and field_bkg needed for feature detection zero_diff_cos_val : float Value where the cosine function returns a zero difference always_core_thres : float - All values above this threshold considered to be convective + All values above this threshold considered to be features CS_CORE : int - Value assigned to convective pixels + Value assigned to features Returns ------- - conv_core_array : array - Array of booleans if point is convective (1) or not (0) + core_array : array + Array of booleans if point is feature (1) or not (0) """ - # initialize entire array to not a convective core - conv_core_array = np.zeros_like(refl) + # initialize entire array to not a core + core_array = np.zeros_like(field) # calculate zeDiff for entire array - zDiff = max_diff * np.cos((np.pi * refl_bkg) / (2 * zero_diff_cos_val)) + zDiff = max_diff * np.cos((np.pi * field_bkg) / (2 * zero_diff_cos_val)) zDiff[zDiff < 0] = 0 # where difference less than zero, set to zero - zDiff[refl_bkg < 0] = max_diff # where background less than zero, set to max. diff + zDiff[field_bkg < 0] = max_diff # where background less than zero, set to max. diff # set core values - # where refl >= always_core_thres and where difference exceeds minimum + # where field >= always_core_thres and where difference exceeds minimum core_elements = np.logical_or( - (refl >= always_core_thres), (refl - refl_bkg) >= zDiff + (field >= always_core_thres), (field - field_bkg) >= zDiff ) core_elements = core_elements.filled(0) - conv_core_array[core_elements] = CS_CORE + core_array[core_elements] = CS_CORE - # mask by refl array - conv_core_array = np.ma.masked_where(refl.mask, conv_core_array) + # mask by field array + core_array = np.ma.masked_where(field.mask, core_array) - return conv_core_array + return core_array -def convcore_scalar_scheme( - refl, refl_bkg, max_diff, always_core_thres, CS_CORE, use_addition=False +def core_scalar_scheme( + field, field_bkg, max_diff, always_core_thres, CORE, use_addition=False ): """ - Function for assigning convective cores based on a scalar difference + Function for assigning cores based on a scalar difference Parameters ---------- - refl : array - Reflectivity values - refl_bkg : array - Background average of reflectivity values + field : array + Field values + field_bkg : array + Background average of field values max_diff : float - Maximum difference between refl and refl_bkg needed for convective classification + Maximum difference between field and field_bkg needed for feature detection always_core_thres : float - All values above this threshold considered to be convective - CS_CORE : int - Value assigned to convective pixels + All values above this threshold considered to be a feature + CORE : int + Value assigned to features use_addition : bool Boolean to determine if scalar should be added (True) or multiplied (False) Returns ------- - conv_core_array : array - Array of booleans if point is convective (1) or not (0) + core_array : array + Array of booleans if point is feature (1) or not (0) """ - # initialize entire array to not a convective core - conv_core_array = np.zeros_like(refl) + # initialize entire array to not a core + core_array = np.zeros_like(field) # calculate zDiff for entire array # if addition, add difference. Else, multiply difference if use_addition: - zDiff = (max_diff + refl_bkg) - refl_bkg + zDiff = (max_diff + field_bkg) - field_bkg else: - zDiff = (max_diff * refl_bkg) - refl_bkg + zDiff = (max_diff * field_bkg) - field_bkg zDiff[zDiff < 0] = 0 # where difference less than zero, set to zero - zDiff[refl_bkg < 0] = 0 # where background less than zero, set to zero + zDiff[field_bkg < 0] = 0 # where background less than zero, set to zero # set core values - # where refl >= always_core_thres and where difference exceeds minimum + # where field >= always_core_thres and where difference exceeds minimum core_elements = np.logical_or( - (refl >= always_core_thres), (refl - refl_bkg) >= zDiff + (field >= always_core_thres), (field - field_bkg) >= zDiff ) core_elements = core_elements.filled(0) - conv_core_array[core_elements] = CS_CORE + core_array[core_elements] = CORE - # mask by refl array - conv_core_array = np.ma.masked_where(refl.mask, conv_core_array) + # mask by field array + core_array = np.ma.masked_where(field.mask, core_array) - return conv_core_array + return core_array -def create_conv_radius_mask( - max_conv_diameter, radius_km, x_spacing, y_spacing, center_conv_mask_x -): +def create_radius_mask(max_diameter, radius_km, x_spacing, y_spacing, center_mask_x): """ - Does and initial convective stratiform classification + Creates a circular mask based on input diameter Parameters ---------- - max_conv_diameter : int - maximum convective diameter in kilometers + max_diameter : int + maximum diameter in kilometers radius_km : int - convective radius in kilometers + radius in kilometers x_spacing, y_spacing : float x- and y-dimension pixel size in meters, respectively - center_conv_mask_x : int + center_mask_x : int index of center point Returns ------- - conv_mask_array : array - array masked based on distance of convective diameter + feature_mask_array : array + array masked based on distance of diameter to incorporate """ - conv_mask_array = np.zeros((max_conv_diameter, max_conv_diameter)) - conv_mask_array = create_radial_mask( - conv_mask_array, + feature_mask_array = np.zeros((max_diameter, max_diameter)) + feature_mask_array = create_radial_mask( + feature_mask_array, 0, radius_km, x_spacing, y_spacing, - center_conv_mask_x, - center_conv_mask_x, + center_mask_x, + center_mask_x, True, ) - return conv_mask_array + return feature_mask_array -def assign_conv_radius_km(refl_bkg, val_for_max_conv_rad, max_conv_rad=5): +def assign_feature_radius_km(field_bkg, val_for_max_rad, max_rad=5): """ - Assigns the convective radius in kilometers based on the background values + Assigns the radius in kilometers based on the background values Parameters ---------- - refl_bkg : array - array of background reflectivity values - val_for_max_conv_rad : float - reflectivity value for maximum convective radius (5 km) - max_conv_rad : float, optional - maximum convective radius in kilometers + field_bkg : array + array of background field values + val_for_max_rad : float + field value for maximum radius (5 km) + max_rad : float, optional + maximum radius in kilometers Returns ------- - convRadiuskm : array - array of convective radii based on background values and val for max. conv radius + radius_array_km : array + array of radii based on background values and val for max. radius """ - convRadiuskm = np.ones_like(refl_bkg) + radius_array_km = np.ones_like(field_bkg) - convRadiuskm[refl_bkg >= (val_for_max_conv_rad - 15)] = max_conv_rad - 3 - convRadiuskm[refl_bkg >= (val_for_max_conv_rad - 10)] = max_conv_rad - 2 - convRadiuskm[refl_bkg >= (val_for_max_conv_rad - 5)] = max_conv_rad - 1 - convRadiuskm[refl_bkg >= val_for_max_conv_rad] = max_conv_rad + radius_array_km[field_bkg >= (val_for_max_rad - 15)] = max_rad - 3 + radius_array_km[field_bkg >= (val_for_max_rad - 10)] = max_rad - 2 + radius_array_km[field_bkg >= (val_for_max_rad - 5)] = max_rad - 1 + radius_array_km[field_bkg >= val_for_max_rad] = max_rad - return convRadiuskm + return radius_array_km -def classify_conv_strat_array( - refl, - conv_strat_array, - conv_core_array, +def classify_feature_array( + field, + feature_array, + core_array, NOSFCECHO, - CONV, - SF, + FEAT_VAL, + BKGD_VAL, WEAKECHO, - CS_CORE, + CORE, MINDBZUSE, WEAKECHOTHRES, ): """ - Does and initial convective stratiform classification + Does an initial feature detection Parameters ---------- - refl : array - Array of reflectivity values - conv_strat_array : array - Array with convective stratiform classifications - conv_core_array : array - Array with convective cores + field : array + Array of field values + feature_array : array + Array with feature detection + core_array : array + Array with cores NOSFCECHO : int Value to assign points classified as no surface echo - CONV : int - Value to assign points classified as convective - SF : int - Value to assign points classified as stratiform + FEAT_VAL : int + Value to assign points classified as features + BKGD_VAL : int + Value to assign points classified as background echo WEAKECHO : int Value to assign points classfied as weak echo - CS_CORE : int - Value assigned to convective cores in conv_core_array + CORE : int + Value assigned to cores in core_array MINDBZUSE : float Minimum dBZ value to consider in classification, all values below this will be set to NOSFCECHO WEAKECHOTHRES : float @@ -796,20 +799,20 @@ def classify_conv_strat_array( Returns ------- - conv_strat_array : array + feature_array : array array of classifications """ # assuming order so that each point is only assigned one time, no overlapping assignment # initially, assign every point to stratiform - conv_strat_array[:] = SF - # where reflectivity is masked, set to no surface echo - conv_strat_array[refl.mask] = NOSFCECHO - # assign convective cores to CONV - conv_strat_array[conv_core_array == CS_CORE] = CONV - # assign reflectivity less than weakechothres to weak echo - conv_strat_array[refl < WEAKECHOTHRES] = WEAKECHO - # assign reflectivity less than minimum to no surface echo - conv_strat_array[refl < MINDBZUSE] = NOSFCECHO - - return conv_strat_array + feature_array[:] = BKGD_VAL + # where field is masked, set to no surface echo + feature_array[field.mask] = NOSFCECHO + # assign cores to FEAT + feature_array[core_array == CORE] = FEAT_VAL + # assign field values less than weakechothres to weak echo + feature_array[field < WEAKECHOTHRES] = WEAKECHO + # assign field values less than minimum to no surface echo + feature_array[field < MINDBZUSE] = NOSFCECHO + + return feature_array diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index f118a2e99a..7a15fa448c 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -8,7 +8,7 @@ import numpy as np from ..config import get_field_name, get_fillvalue, get_metadata -from ._echo_class import _revised_conv_strat, steiner_class_buff +from ._echo_class import _feature_detection, steiner_class_buff def steiner_conv_strat( @@ -245,14 +245,191 @@ def conv_strat_yuter( """ - # Maxmimum convective radius must be less than 5 km - if max_conv_rad_km > 5: - print("Max conv radius must be less than 5 km, exiting") + warn( + "This function will be deprecated in Py-ART 2.0." + " Please use feature_detection function.", + DeprecationWarning, + ) + + feature_dict = feature_detection( + grid, + dx=dx, + dy=dy, + level_m=level_m, + always_core_thres=always_core_thres, + bkg_rad_km=bkg_rad_km, + use_cosine=use_cosine, + max_diff=max_diff, + zero_diff_cos_val=zero_diff_cos_val, + scalar_diff=scalar_diff, + use_addition=use_addition, + calc_thres=calc_thres, + weak_echo_thres=weak_echo_thres, + min_val_used=min_dBZ_used, + dB_averaging=dB_averaging, + remove_small_objects=remove_small_objects, + min_km2_size=min_km2_size, + binary_close=False, + val_for_max_rad=val_for_max_conv_rad, + max_rad_km=max_conv_rad_km, + core_val=cs_core, + nosfcecho=nosfcecho, + weakecho=weakecho, + bkgd_val=sf, + feat_val=conv, + field=refl_field, + estimate_flag=estimate_flag, + estimate_offset=5, + overest_field=None, + underest_field=None, + ) + + return feature_dict + + +def feature_detection( + grid, + dx=None, + dy=None, + level_m=None, + always_core_thres=42, + bkg_rad_km=11, + use_cosine=True, + max_diff=5, + zero_diff_cos_val=55, + scalar_diff=1.5, + use_addition=True, + calc_thres=0.75, + weak_echo_thres=5.0, + min_val_used=5.0, + dB_averaging=True, + remove_small_objects=True, + min_km2_size=10, + binary_close=False, + val_for_max_rad=30, + max_rad_km=5.0, + core_val=3, + nosfcecho=0, + weakecho=3, + bkgd_val=1, + feat_val=2, + field=None, + estimate_flag=True, + estimate_offset=5, + overest_field=None, + underest_field=None, +): + """ + This function can be used to detect features in a field (e.g. reflectivity, rain rate, snow rate, + etc.) described by Tomkins et al. (2023) and based on original convective-stratiform algorithms developed by + Steiner et al. (1995), Yuter et al. (2005) and Yuter and Houze (1997) algorithm. + + Author: Laura Tomkins (@lauratomkins) + + Parameters + ---------- + grid : Grid + Grid containing reflectivity field to partition. + dx, dy : float, optional + The x- and y-dimension resolutions in meters, respectively. If None + the resolution is determined from the first two axes values parsed from grid object. + level_m : float, optional + Desired height in meters to run feature detection algorithm. + always_core_thres : float, optional + Threshold for points that are always features. All values above the threshold are classified as features. + bkg_rad_km : float, optional + Radius to compute background reflectivity in kilometers. Default is 11 km. Recommended to be at least 3 x + grid spacing + use_cosine : bool, optional + Boolean used to determine if a cosine scheme (see Yuter and Houze (1997)) should be used for identifying + cores (True) or if a simpler scalar scheme (False) should be used. + max_diff : float, optional + Maximum difference between background average and reflectivity in order to be classified as features. + "a" value in Eqn. B1 in Yuter and Houze (1997) + zero_diff_cos_val : float, optional + Value where difference between background average and reflectivity is zero in the cosine function + "b" value in Eqn. B1 in Yuter and Houze (1997) + scalar_diff : float, optional + If using a scalar difference scheme, this value is the multiplier or addition to the background average + use_addition : bool, optional + Determines if a multiplier (False) or addition (True) in the scalar difference scheme should be used + calc_thres : float, optional + Minimum percentage of points needed to be considered in background average calculation + weak_echo_thres : float, optional + Threshold for determining weak echo. All values below this threshold will be considered weak echo + min_val_used : float, optional + Minimum value used for classification. All values below this threshold will be considered no surface echo + See Yuter and Houze (1997) and Yuter et al. (2005) for more detail. Units based on input field + dB_averaging : bool, optional + True if using dBZ reflectivity values that need to be converted to linear Z before averaging. False for + other non-dBZ values (i.e. snow rate) + remove_small_objects : bool, optional + Determines if small objects should be removed from core array. Default is True. + min_km2_size : float, optional + Minimum size of Cores to be considered. Cores less than this size will be removed. Default is 10 km^2. + binary_close : bool, optional + Determines if a binary closing should be performed on the cores. Default is False. + val_for_max_rad : float, optional + value used for maximum radius. Cores with values above this will have the maximum radius incorporated. + max_rad_km : float, optional + Maximum radius around cores to classify as feature. Default is 5 km + core_val : int, optional + Value for points classified as cores + nosfcecho : int, optional + Value for points classified as no surface echo, based on min_val_used + weakecho : int, optional + Value for points classified as weak echo, based on weak_echo_thres. + bkgd_val : int, optional + Value for points classified as background echo. + feat_val : int, optional + Value for points classified as features. + field : str, optional + Field in grid to find objects in. None will use the default reflectivity field name from the Py-ART + configuration file. + estimate_flag : bool, optional + Determines if over/underestimation should be applied. If true, the algorithm will also be run on the same field + wih the estimate_offset added and the same field with the estimate_offset subtracted. + Default is True (recommended) + estimate_offset : float, optional + Value used to offset the field values by for the over/underestimation application. Default value is 5 dBZ. + overest_field : str, optional + Name of field in grid object used to calculate the overestimate if estimate_flag is True. + underest_field : str, optional + Name of field in grid object used to calculate the underestimate if estimate_flag is True. + + Returns + ------- + feature_dict : dict + Feature detection classification dictionary. + + References + ---------- + Steiner, M. R., R. A. Houze Jr., and S. E. Yuter, 1995: Climatological + Characterization of Three-Dimensional Storm Structure from Operational + Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007. + + Yuter, S. E., and R. A. Houze, Jr., 1997: Measurements of raindrop size + distributions over the Pacific warm pool and implications for Z-R relations. + J. Appl. Meteor., 36, 847-867. + https://doi.org/10.1175/1520-0450(1997)036%3C0847:MORSDO%3E2.0.CO;2 + + Yuter, S. E., R. A. Houze, Jr., E. A. Smith, T. T. Wilheit, and E. Zipser, + 2005: Physical characterization of tropical oceanic convection observed in + KWAJEX. J. Appl. Meteor., 44, 385-415. https://doi.org/10.1175/JAM2206.1 + + Tomkins, L. M., S. E. Yuter, and M. A. Miller, 2024: Objective identification + of faint and strong features in radar observations of winter storms. in prep. + + """ + + # Maxmimum radius must be less than 5 km + if max_rad_km > 5: + print("Max radius must be less than 5 km, exiting") raise # Parse field parameters - if refl_field is None: - refl_field = get_field_name("reflectivity") + if field is None: + field = get_field_name("reflectivity") dB_averaging = True # parse dx and dy if None @@ -274,15 +451,15 @@ def conv_strat_yuter( # Get reflectivity data at desired level if level_m is None: try: - ze = np.ma.copy(grid.fields[refl_field]["data"][0, :, :]) + ze = np.ma.copy(grid.fields[field]["data"][0, :, :]) except: - ze = np.ma.copy(grid.fields[refl_field]["data"][:, :]) + ze = np.ma.copy(grid.fields[field]["data"][:, :]) else: zslice = np.argmin(np.abs(z - level_m)) - ze = np.ma.copy(grid.fields[refl_field]["data"][zslice, :, :]) + ze = np.ma.copy(grid.fields[field]["data"][zslice, :, :]) - # run convective stratiform algorithm - _, _, convsf_best = _revised_conv_strat( + # run feature detection algorithm + _, _, feature_best = _feature_detection( ze, dx, dy, @@ -295,43 +472,47 @@ def conv_strat_yuter( use_addition=use_addition, calc_thres=calc_thres, weak_echo_thres=weak_echo_thres, - min_dBZ_used=min_dBZ_used, + min_val_used=min_val_used, dB_averaging=dB_averaging, remove_small_objects=remove_small_objects, min_km2_size=min_km2_size, - val_for_max_conv_rad=val_for_max_conv_rad, - max_conv_rad_km=max_conv_rad_km, - cs_core=cs_core, + binary_close=binary_close, + val_for_max_rad=val_for_max_rad, + max_rad_km=max_rad_km, + core_val=core_val, nosfcecho=nosfcecho, weakecho=weakecho, - sf=sf, - conv=conv, + bkgd_val=bkgd_val, + feat_val=feat_val, ) # put data into a dictionary to be added as a field - convsf_dict = { - "convsf": { - "data": convsf_best, - "standard_name": "convsf", - "long_name": "Convective stratiform classification", + feature_dict = { + "feature_detection": { + "data": feature_best[None, :, :], + "standard_name": "feature_detection", + "long_name": "Feature Detection", "valid_min": 0, "valid_max": 3, "comment_1": ( - "Convective-stratiform echo " - "classification based on " - "Yuter and Houze (1997) and Yuter et al. (2005)" - ), - "comment_2": ( - "0 = No surface echo/Undefined, 1 = Stratiform, " - "2 = Convective, 3 = weak echo" + "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format( + nosfcecho, bkgd_val, feat_val, weakecho + ) ), } } # If estimation is True, run the algorithm on the field with offset subtracted and the field with the offset added if estimate_flag: - _, _, convsf_under = _revised_conv_strat( - ze - estimate_offset, + # get underestimate field + if underest_field is None: + under_field = ze - estimate_offset + elif underest_field is not None: + under_field = np.ma.copy(grid.fields[underest_field]["data"][0, :, :]) + + # run algorithm to get feature detection underestimate + _, _, feature_under = _feature_detection( + under_field, dx, dy, always_core_thres=always_core_thres, @@ -343,21 +524,29 @@ def conv_strat_yuter( use_addition=use_addition, calc_thres=calc_thres, weak_echo_thres=weak_echo_thres, - min_dBZ_used=min_dBZ_used, + min_val_used=min_val_used, dB_averaging=dB_averaging, remove_small_objects=remove_small_objects, min_km2_size=min_km2_size, - val_for_max_conv_rad=val_for_max_conv_rad, - max_conv_rad_km=max_conv_rad_km, - cs_core=cs_core, + binary_close=binary_close, + val_for_max_rad=val_for_max_rad, + max_rad_km=max_rad_km, + core_val=core_val, nosfcecho=nosfcecho, weakecho=weakecho, - sf=sf, - conv=conv, + bkgd_val=bkgd_val, + feat_val=feat_val, ) - _, _, convsf_over = _revised_conv_strat( - ze + estimate_offset, + # get overestimate field + if overest_field is None: + over_field = ze + estimate_offset + elif overest_field is not None: + over_field = np.ma.copy(grid.fields[overest_field]["data"][0, :, :]) + + # run algorithm to get feature detection underestimate + _, _, feature_over = _feature_detection( + over_field, dx, dy, always_core_thres=always_core_thres, @@ -369,53 +558,48 @@ def conv_strat_yuter( use_addition=use_addition, calc_thres=calc_thres, weak_echo_thres=weak_echo_thres, - min_dBZ_used=min_dBZ_used, + min_val_used=min_val_used, dB_averaging=dB_averaging, remove_small_objects=remove_small_objects, min_km2_size=min_km2_size, - val_for_max_conv_rad=val_for_max_conv_rad, - max_conv_rad_km=max_conv_rad_km, - cs_core=cs_core, + binary_close=binary_close, + val_for_max_rad=val_for_max_rad, + max_rad_km=max_rad_km, + core_val=core_val, nosfcecho=nosfcecho, weakecho=weakecho, - sf=sf, - conv=conv, + bkgd_val=bkgd_val, + feat_val=feat_val, ) # save into dictionaries - convsf_dict["convsf_under"] = { - "data": convsf_under, - "standard_name": "convsf_under", - "long_name": "Convective stratiform classification (Underestimate)", + feature_dict["feature_under"] = { + "data": feature_under[None, :, :], + "standard_name": "feature_detection_under", + "long_name": "Feature Detection Underestimate", "valid_min": 0, "valid_max": 3, "comment_1": ( - "Convective-stratiform echo " - "classification based on " - "Yuter and Houze (1997) and Yuter et al. (2005)" - ), - "comment_2": ( - "0 = Undefined, 1 = Stratiform, " "2 = Convective, 3 = weak echo" + "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format( + nosfcecho, bkgd_val, feat_val, weakecho + ) ), } - convsf_dict["convsf_over"] = { - "data": convsf_over, - "standard_name": "convsf_under", - "long_name": "Convective stratiform classification (Overestimate)", + feature_dict["feature_over"] = { + "data": feature_over[None, :, :], + "standard_name": "feature_detection_over", + "long_name": "Feature Detection Overestimate", "valid_min": 0, "valid_max": 3, "comment_1": ( - "Convective-stratiform echo " - "classification based on " - "Yuter and Houze (1997)" - ), - "comment_2": ( - "0 = Undefined, 1 = Stratiform, " "2 = Convective, 3 = weak echo" + "{} = No surface echo/Undefined, {} = Background echo, {} = Features, {} = weak echo".format( + nosfcecho, bkgd_val, feat_val, weakecho + ) ), } - return convsf_dict + return feature_dict def hydroclass_semisupervised( diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index a2eacafc5a..1a0246f1d7 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -499,6 +499,8 @@ def image_mute_radar(radar, field, mute_field, mute_threshold, field_threshold=N the correlation coefficient is less than a certain threshold to discern melting precipitation. + Author: Laura Tomkins (@lauratomkins) + Parameters ---------- radar : Radar diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index ff11f74dd3..288ca32dd9 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -67,15 +67,15 @@ def test_steiner_conv_strat_modify_area(area_relation): assert eclass["data"].max() == 2 -def test_conv_strat__yuter_default(): +def test_feature_detection_default(): grid = pyart.testing.make_storm_grid() - dict = pyart.retrieve.conv_strat_yuter(grid, bkg_rad_km=50) + dict = pyart.retrieve.feature_detection(grid, bkg_rad_km=50) - assert "convsf" in dict.keys() - assert "convsf_under" in dict.keys() - assert "convsf_over" in dict.keys() + assert "feature_detection" in dict.keys() + assert "feature_under" in dict.keys() + assert "feature_over" in dict.keys() assert np.all( - dict["convsf"]["data"][25] + dict["feature_detection"]["data"][0, 25, :] == np.array( [ 0, @@ -123,13 +123,13 @@ def test_conv_strat__yuter_default(): ) -def test_conv_strat_yuter_noest(): +def test_feature_detection_noest(): grid = pyart.testing.make_storm_grid() - dict = pyart.retrieve.conv_strat_yuter(grid, bkg_rad_km=50, estimate_flag=False) + dict = pyart.retrieve.feature_detection(grid, bkg_rad_km=50, estimate_flag=False) - assert "convsf" in dict.keys() - assert "convsf_under" not in dict.keys() - assert "convsf_over" not in dict.keys() + assert "feature_detection" in dict.keys() + assert "feature_under" not in dict.keys() + assert "feature_over" not in dict.keys() def test_hydroclass_semisupervised():