Spaceborne-Ground Radar Matching (Methodology)#

In this tutorial, we demonstrate how to exploit GPM-API along with other radar software such xradar and wradlib to match reflectivities measurements of spaceborne (SR) and ground (GR) radars.

We guide you step-by-step through the process of obtaining spatially and temporally coincident radar samples.

The procedure, based on Schwaller and Morris (2011) and adapted by Warren, et. al. (2018), involves:

averaging SR reflectivities vertically along the SR beam between the half-power points of the GR sweep.
averaging GR reflectivities horizontally within the SR beam’s footprint.

The basic principle is illustrated in Fig. 2 of the original paper of Schwaller and Morris (2011).

Warren et al. (2018) describe the method as follows: “[…] intersections between individual SR beams and GR elevation sweeps are identified and the reflectivity values from both instruments are averaged within a spatial neighborhood around the intersection. Specifically, SR data are averaged in range over the width of the GR beam at the GR range of the intersection, while GR data are averaged in the range-azimuth plane within the footprint of the SR beam. The result is a pair of reflectivity measurements corresponding to approximately the same volume of atmosphere. […]”.

The procedure should become clearer in Fig. 3:

Relevant software#

To run this tutorial, it is necessary to install additional libraries which are not automatically required by GPM-API. The libraries are xradar, wradlib, geopandas, fsspec, s3fs and folium. You can install them with terminal command conda install -c conda-forge xradar wradlib geopandas fsspec s3fs folium

If you have such libraries installed, you can start importing the required packages:

[ ]:

import warnings

warnings.filterwarnings("ignore")
import datetime
from functools import reduce

import cartopy.crs as ccrs
import fsspec
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import wradlib as wrl
import xarray as xr
import xradar as xd
from IPython.display import display
from shapely import Point

import gpm
from gpm.gv.routines import (
    convert_s_to_ku_band,
    retrieve_gates_projection_coordinates,
    xyz_to_antenna_coordinates,
)
from gpm.utils.remapping import reproject_coords
from gpm.utils.zonal_stats import PolyAggregator

np.set_printoptions(suppress=True)

Retrieve GR gate resolution#

In this subsection we derive the horizontal and vertical resolution of GR gates.

The GR horizontal and vertical resolution change as a function of the range distance.

[ ]:

# Compute horizontal and vertical GR gate resolution
h_res_gr, v_res_gr = ds_gr.xradar_dev.resolution_at_range(
    azimuth_beamwidth=azimuth_beamwidth_gr,
    elevation_beamwidth=elevation_beamwidth_gr,
)
h_res_gr.xradar_dev.plot_map()
v_res_gr.xradar_dev.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe7d8c652d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_42_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_42_2.png

Define GR range distance interval#

Now let’s define the minimum range distance at which the GR gate vertical resolution is larger than 150/250 m (the range resolution of GPM) and at which GR distance the GR gate vertical resolution exceed 1500 m (which correspond to 12/6 SR gates included in a GR gate)

[ ]:

mask_gr_gate_depth_above_250 = v_res_gr > 250
mask_gr_gate_depth_above_250.xradar_dev.plot_map()
min_gr_range = np.max(
    mask_gr_gate_depth_above_250["range"].data[mask_gr_gate_depth_above_250.sum(dim="azimuth") == 0],
)  # min_range_distance

mask_gr_gate_depth_below_2000 = v_res_gr > 1500
mask_gr_gate_depth_below_2000.xradar_dev.plot_map()
max_gr_range = np.max(
    mask_gr_gate_depth_below_2000["range"].data[mask_gr_gate_depth_below_2000.sum(dim="azimuth") == 0],
)  # min_range_distance

print("Minimum GR distance:", min_gr_range, "m")
print("Maximum GR distance:", max_gr_range, "m")

Minimum GR distance: 14125.0 m
Maximum GR distance: 85875.0 m

../_images/tutorials_tutorial_03_SR_GR_Matching_45_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_45_2.png

Now let’s specify the minimum and maximum GR range distance to consider when matching SR/GR gates.

It’s good practice to avoid restricting too much at this stage. Matched SR/GR gates can be filtered afterwards.

[ ]:

# Specify GR range distance interval over which to perform SR matching
min_gr_range = 7_000  # about 15_000 km
max_gr_range = 150_000

Display SR slice and GR sweep#

[ ]:

# Slice SR volume
da_sr = ds_sr["zFactorFinal"].isel(range=-5)
da_sr = ds_sr["zFactorFinal"].gpm.slice_range_at_bin(ds_sr["binClutterFreeBottom"])
da_sr = ds_sr["zFactorFinal"].gpm.slice_range_at_height(value=2000)

[ ]:

# Define SR and GR mask
mask_gr = ds_gr["DBZH"] >= z_min_threshold_gr
# mask_gr = np.logical_and(ds_gr["DBZH"] >= z_min_threshold_gr, ds_gr["RHOHV"] >= 0.80)

mask_sr = True  # aka do not mask

# Retrieve SR extent
sr_extent = ds_sr.gpm.extent()

# Define cmap
cmap = "Spectral_r"

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), dpi=120, subplot_kw={"projection": ccrs.PlateCarree()})
# - Plot SR
p_sr = da_sr.where(mask_sr).gpm.plot_map(
    ax=ax1,
    cmap=cmap,
    vmin=0,
    vmax=40,
    cbar_kwargs={"label": "SR Reflectivity (dBz)"},
    add_colorbar=True,
)
# - Plot GR
ds_gr["DBZH"].where(mask_gr).xradar_dev.plot_map(
    ax=ax2,
    cmap=cmap,
    vmin=0,
    vmax=40,
    cbar_kwargs={"label": "GR Reflectivity (dBz)", "extend": "both"},
    add_colorbar=True,
)
# - Add the SR swath boundary
da_sr.gpm.plot_swath_lines(ax=ax2)
# - Restrict the extent to the SR overpass
ax2.set_extent(sr_extent)
# - Display GR range distances
for ax in [ax1, ax2]:
    ds_gr.xradar_dev.plot_range_distance(
        distance=min_gr_range,
        ax=ax,
        add_background=True,
        linestyle="dashed",
        edgecolor="black",
    )
    ds_gr.xradar_dev.plot_range_distance(
        distance=max_gr_range,
        ax=ax,
        add_background=True,
        linestyle="dashed",
        edgecolor="black",
    )
    ds_gr.xradar_dev.plot_range_distance(
        distance=250_000,
        ax=ax,
        add_background=True,
        linestyle="dashed",
        edgecolor="black",
    )
# - Add GR location
ax1.scatter(ds_gr["longitude"], ds_gr["latitude"], c="black", marker="X", s=4)
ax2.scatter(ds_gr["longitude"], ds_gr["latitude"], c="black", marker="X", s=4)
# - Set title
ax1.set_title("GPM", fontsize=12, loc="left")
ax2.set_title("Ground Radar", fontsize=12, loc="left")
# - Improve layout and display
plt.tight_layout()
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_50_0.png

Retrieve GR (lon, lat) coordinates#

Here we retrieve the 2D non-dimensional longitude and latitude coordinates in the CRS of the SR. These 2D coordinates will have dimensions (azimuth, range)

[ ]:

lon_gr, lat_gr, height_gr = reproject_coords(x=ds_gr["x"], y=ds_gr["y"], z=ds_gr["z"], src_crs=crs_gr, dst_crs=crs_sr)

# If crs_gr and crs_sr share same datum/ellispoid, z does not vary
print(crs_gr.datum)
print(crs_sr.datum)

World Geodetic System 1984
Unknown based on WGS 84 ellipsoid

Retrieve SR (x,y,z) coordinates#

Here we retrieve the 3D (x,y,z) coordinates of SR gates in the GR CRS projection.

The GR CRS is a azimuthal equidistant (AZEQ) projection centered on the ground radar. In the GR CRS, the (x,y) coordinates (0,0) corresponds to the location of the ground radar. Therefore, the retrieved SR (x,y) coordinates are relative to the ground radar site, while the z coordinates is relative to the ellipsoid.

[ ]:

x_sr, y_sr, z_sr = retrieve_gates_projection_coordinates(ds_sr, dst_crs=crs_gr)

Retrieve SR (range, azimuth, elevation) coordinates#

Now we inverse the SR (x,y,z) 3D GR CRS coordinates to derive the SR (range, azimuth, elevation) coordinates with respect to GR sweep. These coordinates allow to identify the SR gates intersecting the GR sweep.

[ ]:

range_sr, azimuth_sr, elevation_sr = xyz_to_antenna_coordinates(
    x=x_sr,
    y=y_sr,
    z=z_sr,
    site_altitude=ds_gr["altitude"],
    crs=crs_gr,
    effective_radius_fraction=None,
)
print(range_sr.isel(cross_track=0).min().item(), "m")  # far from the radar
print(range_sr.isel(cross_track=-1).min().item(), "m")  # closer to radar

252837.9704755049 m
9164.066906023296 m

[ ]:

# Show elevation angle at different SR scan angle
elevation_sr.isel(cross_track=0).gpm.plot_cross_section(vmin=0)  # far from the radar
elevation_sr.isel(cross_track=-1).gpm.plot_cross_section(vmin=0)  # closer to radar

<matplotlib.collections.QuadMesh at 0x7fe7d84fae90>

../_images/tutorials_tutorial_03_SR_GR_Matching_58_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_58_2.png

Retrieve SR range distance#

Here we retrieve the range distance from the satellite for each SR radar gate. We can accurately retrieve it from the GPM L1B product using the following code:

[ ]:

# Retrieve L1B accurate range distance
ds_l1b_sr["range_distance_from_satellite"] = ds_l1b_sr.gpm.retrieve("range_distance_from_satellite")
ds_sr["range_distance_from_satellite"] = ds_l1b_sr.gpm.extract_l2_dataset()["range_distance_from_satellite"]

The range_distance_from_satellite can also be derived from the GPM L2 product if we assume a cross-track scan angle. However, uncertainties in cross-tracks scan angle can erroneously impacts the estimates. As we show here below, an uncertainty of 0.1° translates to a range distance difference of up to 400 m which corresponds to 3/4 radar gates !

[ ]:

# Show impact of scan angle on range distance
# - Variation of 0.01 impacts range distance up to 30 m
# - Variation of 0.05 impacts range distance up to 100 m
# - Variation of 0.1 impacts range distance up to 400 m
scan_angle_template = np.abs(np.linspace(-17.04, 17.04, 49))  # assuming Ku swath of size 49
scan_angle = xr.DataArray(scan_angle_template[ds_sr["gpm_cross_track_id"].data], dims="cross_track")
scan_angle = scan_angle + 0.1
l2_distance = ds_sr.gpm.retrieve("range_distance_from_satellite", scan_angle=scan_angle)

diff = l2_distance - ds_sr["range_distance_from_satellite"]
diff.isel(along_track=0).gpm.plot_cross_section()  # difference > 500 m

<matplotlib.collections.QuadMesh at 0x7fe7d82f0650>

../_images/tutorials_tutorial_03_SR_GR_Matching_62_1.png

Retrieve radar gates volumes#

Here we compute SR and GR gate volumes.

An accurate beam width is necessary to derive accurate gate volume estimates. This is demonstrated here below !

[ ]:

# Compute GR gates volumes
vol_gr = wrl.qual.pulse_volume(
    ds_gr["range"],  # range distance
    h=ds_gr["range"].diff("range").median(),  # range resolution
    theta=elevation_beamwidth_gr,
)  # beam width
vol_gr = vol_gr.broadcast_like(ds_gr["DBZH"])
vol_gr = vol_gr.assign_coords({"lon": ds_gr["lon"], "lat": ds_gr["lat"]})

[ ]:

# Compute SR gates volumes
range_distance = ds_sr["range_distance_from_satellite"]
vol_sr = ds_sr.gpm.retrieve("gate_volume", beam_width=beam_width_sr, range_distance=range_distance)
vol_sr.isel(along_track=0).gpm.plot_cross_section()
vol_sr.isel(range=-1).gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe7963055d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_65_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_65_2.png

[ ]:

# Assuming fixed beamwidth with SR leads to volume difference up to 0.35 km3
vol_sr1 = ds_sr.gpm.retrieve("gate_volume", beam_width=0.71, range_distance=range_distance)
diff = vol_sr - vol_sr1
diff.isel(along_track=0).gpm.plot_cross_section()
diff.isel(range=-1).gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe796286710>

../_images/tutorials_tutorial_03_SR_GR_Matching_66_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_66_2.png

Retrieve SR gate resolution#

In this subsection we derive the horizontal and vertical SR gate resolution.

The SR horizontal resolution change for each radar gate, while the vertical resolution varies only across SR footprints. Note that the vertical resolution here derived is equivalent to the GPM L2 product variable height.

Assuming a fixed beamwidth of 0.71 can leads to variations in horizontal resolution estimates up to 400 m on the outer-track.

[ ]:

# Compute horizontal and vertical SR gate resolution
range_distance = ds_sr["range_distance_from_satellite"]
h_res_sr, v_res_sr = ds_sr.gpm.retrieve("gate_resolution", beam_width=beam_width_sr, range_distance=range_distance)
h_radius = h_res_sr / 2

[ ]:

# Validate v_res_sr [OK !]
diff = -ds_sr["height"].diff(dim="range") - v_res_sr
diff.isel(along_track=0).gpm.plot_cross_section(y="range")
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_69_0.png

[ ]:

# Validate h_res [PROBLEMATIC WITH FIXED BEAM WIDTH]
h_res1_sr, _ = ds_sr.gpm.retrieve("gate_resolution", beam_width=None, range_distance=range_distance)  # Use 0.71
diff = h_res_sr - h_res1_sr
diff.isel(along_track=0).gpm.plot_cross_section(y="range")
diff.isel(range=-1).gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe7960e5350>

../_images/tutorials_tutorial_03_SR_GR_Matching_70_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_70_2.png

Retrieve Bright Band (BB) Ratio#

Here we derive the bright band ratio, which is defined as follow:

A BB ratio of < 0 indicates that a bin is located below the melting layer (ML).
A BB ratio of > 0 indicates that a bin is located above the ML.
A BB ratio with values in between 0 and 1 indicates tha the radar is inside the ML.

[ ]:

da_bb_ratio, da_bb_mask = ds_sr.gpm.retrieve("bright_band_ratio", return_bb_mask=True)
da_bb_ratio.isel(cross_track=24).gpm.plot_cross_section(vmin=0, vmax=1)
da_bb_ratio.isel(range=0).gpm.plot_map()  # top
da_bb_ratio.isel(range=-1).gpm.plot_map()  # bottom

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe796223510>

../_images/tutorials_tutorial_03_SR_GR_Matching_72_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_72_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_72_3.png

[ ]:

# Show difference between L2 qualityBB and custom derived BB mask
ds_sr["qualityBB"].gpm.plot_map()
da_bb_mask.gpm.plot_map()
ds_sr["qualityBB"].attrs

{'gpm_api_product': '2A-DPR',
 'flag_values': [0, 1],
 'flag_meanings': ['BB not detected', 'BB detected'],
 'description': 'Quality of Bright-Band Detection',
 'gpm_api_decoded': 'yes',
 'grid_mapping': 'crsWGS84'}

../_images/tutorials_tutorial_03_SR_GR_Matching_73_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_73_2.png

Retrieve Precipitation and Hydrometeors Types#

[ ]:

# Retrieve precipitation type classification
ds_sr["flagPrecipitationType"] = ds_sr.gpm.retrieve("flagPrecipitationType", method="major_rain_type")
ds_sr["flagPrecipitationType"].gpm.plot_map()
ds_sr["flagPrecipitationType"].attrs

{'flag_values': [0, 1, 2, 3],
 'flag_meanings': ['no rain', 'stratiform', 'convective', 'other'],
 'description': 'Precipitation type'}

../_images/tutorials_tutorial_03_SR_GR_Matching_75_1.png

[ ]:

# Retrieve hydrometeors classification
ds_sr["flagHydroClass"] = ds_sr.gpm.retrieve("flagHydroClass")
ds_sr["flagHydroClass"].isel(cross_track=24).gpm.plot_cross_section()
ds_sr["flagHydroClass"].attrs

{'flag_values': [0, 1, 2, 3, 4, 5],
 'flag_meanings': ['no precipitation',
  'rain',
  'snow',
  'hail',
  'melting layer',
  'clutter']}

../_images/tutorials_tutorial_03_SR_GR_Matching_76_1.png

[ ]:

# Visualize cross-section
precip_type3d = xr.ones_like(ds_sr["flagHydroClass"]) * ds_sr["flagPrecipitationType"]
precip_type3d = precip_type3d.where(ds_sr["zFactorFinal"] > 10)
precip_type3d.name = ds_sr["flagPrecipitationType"].name

precip_type3d.isel(cross_track=24).gpm.plot_cross_section()
ds_sr["flagHydroClass"].where(ds_sr["flagHydroClass"] > 0).isel(cross_track=24).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe79d880c50>

../_images/tutorials_tutorial_03_SR_GR_Matching_77_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_77_2.png

Convert SR Ku-band to S-band#

Here we convert reflectivities from Ku-band to S-band based on Cao et al., (2013).

When filtering the SR/GR matched samples, plase consider to select:

only reflectivities below 35 dBZ for S-band GR radars to avoid non-rayleigh scattering
only reflectivities below 25 dBZ for C/X-band GR radars

[ ]:

# With attenuation corrected reflectivity
da_z_s_final = ds_sr.gpm.retrieve("s_band_cao2013", reflectivity="zFactorFinal", bb_ratio=None, precip_type=None)
da_z_ku_final = ds_sr["zFactorFinal"]
dfr_s_ku_final = da_z_s_final - da_z_ku_final
dfr_s_ku_final.name = "DRF_S_Ku"

dfr_s_ku_final.isel(cross_track=24).gpm.plot_cross_section()
da_z_ku_final.isel(cross_track=24).gpm.plot_cross_section()
da_z_s_final.isel(cross_track=24).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe79dd6e7d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_80_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_80_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_80_3.png

[ ]:

# With measured reflectivity
da_z_s_measured = ds_sr.gpm.retrieve("s_band_cao2013", reflectivity="zFactorMeasured", bb_ratio=None, precip_type=None)
da_z_ku_measured = ds_sr["zFactorMeasured"]
dfr_s_ku_measured = da_z_s_measured - da_z_ku_measured
dfr_s_ku_measured.name = "DRF_S_Ku"

dfr_s_ku_measured.isel(cross_track=24).gpm.plot_cross_section()
da_z_ku_measured.isel(cross_track=24).gpm.plot_cross_section()
da_z_s_measured.isel(cross_track=24).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe79e35c3d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_81_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_81_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_81_3.png

[ ]:

# Compute attenuation correction in S-band
da_z_s_correction = da_z_s_final - da_z_s_measured
da_z_s_correction.isel(cross_track=24).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe79e6793d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_82_1.png

Convert SR Ku-band to X and C-band#

Here we convert reflectivities from Ku-band to X and C-band based on XXX et al., (XXX).

We welcome contributions to improve this portion of code !

[ ]:

da_z_c = ds_sr.gpm.retrieve("c_band_tan", bb_ratio=None, precip_type=None)
da_z_x = ds_sr.gpm.retrieve("x_band_tan", bb_ratio=None, precip_type=None)

da_z_c.isel(cross_track=24).gpm.plot_cross_section()
da_z_x.isel(cross_track=24).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe81c737c50>

../_images/tutorials_tutorial_03_SR_GR_Matching_84_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_84_2.png

Convert GR S-band to Ku-band#

Here we provide some code to convert GR S-band reflectivities to Ku-band.

However, the resulting reflectivities are not used in the rest of the tutorial.

The methodology we adopt currently use only an average bright band height to discern between liquid and solid phase precipitation, and it does not account for mixed-phase precipitation and the phase transition.

We welcome contributions to improve this portion of code ;)

[ ]:

# Estimate bright band height
bright_band_height = np.nanmean(ds_sr["heightBB"])
# Retrieve Ku-band reflectivities from GR S-band reflectivities
da_gr_ku = convert_s_to_ku_band(ds_gr=ds_gr, bright_band_height=bright_band_height)

da_gr_s = ds_gr["DBZH"]
da_gr_s.xradar_dev.plot_map(vmin=0, vmax=50)
plt.show()

da_gr_ku.xradar_dev.plot_map(vmin=0, vmax=50)
plt.show()

diff = da_gr_ku - da_gr_s
diff.rename("diff").xradar_dev.plot_map()
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_86_0.png

../_images/tutorials_tutorial_03_SR_GR_Matching_86_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_86_2.png

Identify SR gates intersecting GR sweep#

Here we define the mask used to select SR footprints intersecting the GR sweep gates. Using other terms, we can describe this task as finding SR gates that, once aggregated, “match” with the GR sweep plane-position-indicator (PPI) display.

To this end, we need to filter SR gates by GR range distance and elevation angle.

Since the GR CRS is centered on the GR radar location, we can define the SR range mask using the “GR” range coordinate. So here below, we define the SR range mask with GR range between min_gr_range_lb and max_gr_range_ub. We suggest to specify min_gr_range_lb and max_gr_range_ub not to restrictively, and rather perform more aggressive filtering once the SR/GR matched database is obtained.

The SR elevation mask is obtained by selecting SR gates with elevation angle between the GR sweep elevation angle minus/plus its half beamwidth.

Important note:

Here we currently use the GR range computed at the SR gate centroids (instead of at the SR gate lower and upper bound)
Here we currently use the GR elevation computed at the SR gate centroids (instead of at the SR gate lower and upper bound)

[ ]:

# Define mask of SR footprints within GR range
r = np.sqrt(x_sr**2 + y_sr**2)
mask_sr_within_gr_range = r < ds_gr["range"].max().item()
mask_sr_within_gr_range.isel(range=-1).gpm.plot_map()

mask_sr_within_gr_range_interval = (r >= min_gr_range_lb) & (r <= max_gr_range_ub)

../_images/tutorials_tutorial_03_SR_GR_Matching_88_0.png

[ ]:

# Show how elevation angle (of GR) varies along SR scan angles
elevation_sr.isel(cross_track=0).gpm.plot_cross_section(vmin=0, y="range")
elevation_sr.isel(cross_track=24).gpm.plot_cross_section(vmin=0, y="range")
elevation_sr.isel(cross_track=-1).gpm.plot_cross_section(vmin=0, y="range")

<matplotlib.image.AxesImage at 0x7fe7d835d790>

../_images/tutorials_tutorial_03_SR_GR_Matching_89_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_89_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_89_3.png

[ ]:

# Show SR/GR matched elevation
mask_sr_matched_elevation = (elevation_sr >= (ds_gr["sweep_fixed_angle"] - elevation_beamwidth_gr / 2.0)) & (
    elevation_sr <= (ds_gr["sweep_fixed_angle"] + elevation_beamwidth_gr / 2.0)
)
mask_sr_matched_elevation.isel(cross_track=0).gpm.plot_cross_section()
mask_sr_matched_elevation.isel(cross_track=24).gpm.plot_cross_section()
mask_sr_matched_elevation.isel(cross_track=-1).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe79e247290>

../_images/tutorials_tutorial_03_SR_GR_Matching_90_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_90_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_90_3.png

[ ]:

# Define mask of SR footprints matching GR gates
mask_ppi = mask_sr_matched_elevation & mask_sr_within_gr_range_interval

[ ]:

# Show mask (as SR map)
mask_ppi.any(dim="range").gpm.plot_map(vmin=0)
mask_ppi.isel(range=-10).gpm.plot_map(vmin=0)
mask_ppi.isel(range=-30).gpm.plot_map(vmin=0)

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe79dc59ed0>

../_images/tutorials_tutorial_03_SR_GR_Matching_92_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_92_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_92_3.png

[ ]:

# Show mask (as SR cross-section)
mask_ppi.isel(cross_track=24).gpm.plot_cross_section()
mask_ppi.isel(cross_track=-1).gpm.plot_cross_section()
mask_ppi.isel(cross_track=-1, along_track=slice(10, 60), range=slice(-40, None)).gpm.plot_cross_section()

<matplotlib.collections.QuadMesh at 0x7fe795dab1d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_93_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_93_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_93_3.png

Define SR Masks#

Here we start by defining and analysing various SR masks that can be used to subset SR footprints to be matched with the GR sweep.

[ ]:

# Select scan with "normal" dataQuality (for entire cross-track scan)
mask_sr_quality = ds_sr["dataQuality"] == 0

[ ]:

# Select beams with detected precipitation
ds_sr["flagPrecip"].gpm.plot_map()
mask_sr_precip = ds_sr["flagPrecip"] > 0
mask_sr_precip.name = "mask_precip"
mask_sr_precip.gpm.plot_map()
ds_sr["flagPrecip"].attrs

{'gpm_api_product': '2A-DPR',
 'flag_values': [0, 1, 2, 10, 11, 12, 20, 21, 22],
 'flag_meanings': ['not detected by both Ku and Ka',
  'detected by Ka only',
  'detected by Ka only',
  'detected by Ku only',
  'detected by both Ku and Ka',
  'detected by both Ku and Ka',
  'detected by both Ku and Ka',
  'detected by both Ku and Ka',
  'detected by both Ku and Ka'],
 'description': 'Flag for precipitation detection',
 'gpm_api_decoded': 'yes',
 'grid_mapping': 'crsWGS84'}

../_images/tutorials_tutorial_03_SR_GR_Matching_96_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_96_2.png

[ ]:

# Select above minimum SR reflectivity
mask_sr_minimum_z = ds_sr["zFactorFinal"] > z_min_threshold_sr
mask_sr_minimum_z.name = "Minimum_Reflectitvity"
mask_sr_minimum_z.any(dim="range").gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe795b60b50>

../_images/tutorials_tutorial_03_SR_GR_Matching_97_1.png

[ ]:

# Select only 'normal' data
# mask_sr_quality_data = ds_sr["qualityData"] == 0
# mask_sr_quality_data.gpm.plot_map()

[ ]:

# Select only 'high quality' data
# - qualityFlag derived from qualityData.
# - qualityFlag == 1 indicates low quality retrievals
# - qualityFlag == 2 indicates bad/missing retrievals
mask_sr_quality_data = ds_sr["qualityFlag"] == 0

[ ]:

# Select only beams with detected bright band
ds_sr["qualityBB"].gpm.plot_map()
mask_sr_quality_bb = ds_sr["qualityBB"] == 1
mask_sr_quality_bb.name = "mask_bb"
mask_sr_quality_bb.gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe795962950>

../_images/tutorials_tutorial_03_SR_GR_Matching_100_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_100_2.png

[ ]:

# Select only beams with confident precipitation type
ds_sr["qualityTypePrecip"].gpm.plot_map()
mask_sr_quality_precip = ds_sr["qualityTypePrecip"] == 1
mask_sr_quality_precip.name = "mask_precip_type"
mask_sr_quality_precip.gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe7958b5150>

../_images/tutorials_tutorial_03_SR_GR_Matching_101_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_101_2.png

[ ]:

# Select only stratiform precipitation
ds_sr["flagPrecipitationType"].gpm.plot_map()
mask_sr_stratiform = ds_sr["flagPrecipitationType"] == 1
mask_sr_stratiform.name = "mask_stratiform"
mask_sr_stratiform.gpm.plot_map()

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe795625890>

../_images/tutorials_tutorial_03_SR_GR_Matching_102_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_102_2.png

[ ]:

# Select only beams with reduced path attenuation
# ...
ds_sr["piaFinal"].gpm.plot_map()  # PIA from DSD module
ds_sr["pathAtten"].gpm.plot_map()  # PIAeff
ds_sr["reliabFlag"].gpm.plot_map()  # reliability of PIAeff

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe7953fcd90>

../_images/tutorials_tutorial_03_SR_GR_Matching_103_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_103_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_103_3.png

[ ]:

# Select only beams with limited attenuation correction
# ...
ds_sr["zFactorCorrection"].isel(cross_track=24).gpm.plot_cross_section(vmin=0, vmax=3, zoom=False)
ds_sr["zFactorCorrection"].isel(range=-20).gpm.plot_map(vmin=0, vmax=6)

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7fe79581e250>

../_images/tutorials_tutorial_03_SR_GR_Matching_104_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_104_2.png

Here we define the actual SR mask used to select the SR footprints on which the SR/GR matching will be performed. Once again, we suggest to avoid, at this stage, being too restrictive in the selection of SR footprints.

[ ]:

# Define 3D mask of SR gates matching GR PPI gates
# - TIP: do not mask eccessively here ... but rather at final stage
mask_matched_ppi_3d = mask_ppi & mask_sr_precip & mask_sr_quality

# Define 2D mask of SR beams matching the GR PPI
mask_matched_ppi_2d = mask_matched_ppi_3d.any(dim="range")
mask_matched_ppi_2d.gpm.plot_map()

# Number of matching SR gates per beam
n_matching_sr_gates = mask_matched_ppi_3d.sum(dim="range")
n_matching_sr_gates.where(n_matching_sr_gates > 0).gpm.plot_map()

# Beam indices with at least one matching bin
cross_track_indices, along_track_indices = np.where(mask_matched_ppi_3d.sum(dim="range"))

# Number of beams with at least one valid bin
len(along_track_indices)

../_images/tutorials_tutorial_03_SR_GR_Matching_106_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_106_2.png

Aggregate SR radar gates#

In this subsection, we aggregate SR radar gates intersecting the GR sweep.

SR reflectivities are averaged along the SR beam (approximately vertically) between the half-power points of the GR sweep.

To select SR radar gates for aggregation, the following criteria must be met:

The SR footprints must contain precipitation.
The SR radar gates must overlap with GR gates.
The SR radar gates must fall within a specified GR range distance, defined by min_gr_range and max_gr_range.

The min_gr_range ensures that the GR beamwidth is not smaller than the SR gate spacing.

Important Notes:

Reflectivity must not be aggregated in dBZ, but rather in mm6 m−3

[ ]:

# Add variables to SR dataset
ds_sr["zFactorMeasured_Ku"] = ds_sr["zFactorMeasured"]
ds_sr["zFactorFinal_Ku"] = ds_sr["zFactorFinal"]
ds_sr["zFactorFinal_S"] = da_z_s_final
ds_sr["zFactorMeasured_S"] = da_z_s_measured
ds_sr["zFactorCorrection_Ku"] = ds_sr["zFactorCorrection"]
ds_sr["zFactorCorrection_S"] = da_z_s_correction
ds_sr["hres"] = h_res_sr
ds_sr["vres"] = v_res_sr
ds_sr["gate_volume"] = vol_sr
ds_sr["x"] = x_sr
ds_sr["y"] = y_sr
ds_sr["z"] = z_sr

# Compute path-integrated reflectivities
ds_sr["zFactorFinal_Ku_cumsum"] = ds_sr["zFactorFinal"].gpm.idecibel.cumsum("range").gpm.decibel
ds_sr["zFactorFinal_Ku_cumsum"] = ds_sr["zFactorFinal_Ku_cumsum"].where(np.isfinite(ds_sr["zFactorFinal_Ku_cumsum"]))
ds_sr["zFactorMeasured_Ku_cumsum"] = ds_sr["zFactorMeasured"].gpm.idecibel.cumsum("range").gpm.decibel
ds_sr["zFactorMeasured_Ku_cumsum"] = ds_sr["zFactorMeasured_Ku_cumsum"].where(
    np.isfinite(ds_sr["zFactorMeasured_Ku_cumsum"]),
)
ds_sr["zFactorFinal_S_cumsum"] = ds_sr["zFactorFinal_S"].gpm.idecibel.cumsum("range").gpm.decibel
ds_sr["zFactorFinal_S_cumsum"] = ds_sr["zFactorFinal_S_cumsum"].where(np.isfinite(ds_sr["zFactorFinal_S_cumsum"]))
ds_sr["zFactorMeasured_S_cumsum"] = ds_sr["zFactorMeasured_S"].gpm.idecibel.cumsum("range").gpm.decibel
ds_sr["zFactorMeasured_S_cumsum"] = ds_sr["zFactorMeasured_S_cumsum"].where(
    np.isfinite(ds_sr["zFactorMeasured_S_cumsum"]),
)

# Add variables to the spaceborne dataset
z_variables = ["zFactorFinal_Ku", "zFactorFinal_S", "zFactorMeasured_Ku", "zFactorMeasured_S"]
sr_variables = [
    *z_variables,
    "zFactorCorrection_Ku",
    "zFactorCorrection_S",
    "precipRate",
    "airTemperature",
    "zFactorFinal_Ku_cumsum",
    "zFactorMeasured_Ku_cumsum",
    "zFactorFinal_S_cumsum",
    "zFactorMeasured_S_cumsum",
    "hres",
    "vres",
    "gate_volume",
    "x",
    "y",
    "z",
]

# Initialize Dataset where to add aggregated SR gates
ds_sr_match_ppi = xr.Dataset()

# Mask SR gates not matching the GR PPI
ds_sr_ppi = ds_sr[sr_variables].where(mask_matched_ppi_3d)

# Compute aggregation statistics
ds_sr_ppi_min = ds_sr_ppi.min("range")
ds_sr_ppi_max = ds_sr_ppi.max("range")
ds_sr_ppi_sum = ds_sr_ppi.sum(dim="range", skipna=True)
ds_sr_ppi_mean = ds_sr_ppi.mean("range")
ds_sr_ppi_std = ds_sr_ppi.std("range")

# Aggregate reflectivities (in mm6/mm3)
for var in z_variables:
    ds_sr_ppi_mean[var] = ds_sr_ppi[var].gpm.idecibel.mean("range").gpm.decibel
    ds_sr_ppi_std[var] = ds_sr_ppi[var].gpm.idecibel.std("range").gpm.decibel
    # If only 1 value, std=0, log transform become -inf --> Set to 0
    is_inf = np.isinf(ds_sr_ppi_std[var])
    ds_sr_ppi_std[var] = ds_sr_ppi_std[var].where(~is_inf, 0)

# Compute counts and fractions above sensitivity thresholds
sr_sensitivity_thresholds = [10, 12, 14, 16, 18]
ds_sr_match_ppi["SR_counts"] = mask_matched_ppi_3d.sum(dim="range")
ds_sr_match_ppi["SR_counts_valid"] = (~np.isnan(ds_sr_ppi["zFactorFinal_Ku"])).sum(dim="range")
for var in z_variables:
    for thr in sr_sensitivity_thresholds:
        fraction = (ds_sr_ppi[var] >= thr).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
        ds_sr_match_ppi[f"SR_{var}_fraction_above_{thr}dBZ"] = fraction

# Compute fraction of hydrometeor types
da_hydro_class = ds_sr["flagHydroClass"].where(mask_matched_ppi_3d)
ds_sr_match_ppi["SR_fraction_no_precip"] = (da_hydro_class == 0).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_rain"] = (da_hydro_class == 1).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_snow"] = (da_hydro_class == 2).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_hail"] = (da_hydro_class == 3).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_melting_layer"] = (da_hydro_class == 4).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_clutter"] = (da_hydro_class == 5).sum(dim="range") / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_below_isotherm"] = (ds_sr_ppi["airTemperature"] >= 273.15).sum(
    dim="range",
) / ds_sr_match_ppi["SR_counts"]
ds_sr_match_ppi["SR_fraction_above_isotherm"] = (ds_sr_ppi["airTemperature"] < 273.15).sum(
    dim="range",
) / ds_sr_match_ppi["SR_counts"]

# Add aggregation statistics
for var in ds_sr_ppi_mean.data_vars:
    ds_sr_match_ppi[f"SR_{var}_mean"] = ds_sr_ppi_mean[var]
for var in ds_sr_ppi_std.data_vars:
    ds_sr_match_ppi[f"SR_{var}_std"] = ds_sr_ppi_std[var]
for var in ds_sr_ppi_sum.data_vars:
    ds_sr_match_ppi[f"SR_{var}_sum"] = ds_sr_ppi_sum[var]
for var in ds_sr_ppi_max.data_vars:
    ds_sr_match_ppi[f"SR_{var}_max"] = ds_sr_ppi_max[var]
for var in ds_sr_ppi_min.data_vars:
    ds_sr_match_ppi[f"SR_{var}_min"] = ds_sr_ppi_min[var]
for var in ds_sr_ppi_min.data_vars:
    ds_sr_match_ppi[f"SR_{var}_range"] = ds_sr_ppi_max[var] - ds_sr_ppi_min[var]

# Compute coefficient of variation
for var in z_variables:
    ds_sr_match_ppi[f"SR_{var}_cov"] = ds_sr_match_ppi[f"SR_{var}_std"] / ds_sr_match_ppi[f"SR_{var}_mean"]

# Add SR L2 variables (useful for final filtering and analysis)
var_l2 = [
    "flagPrecip",
    "flagPrecipitationType",
    "dataQuality",
    "sunLocalTime",
    "SCorientation",
    "qualityFlag",
    "qualityTypePrecip",
    "qualityBB",
    "pathAtten",
    "piaFinal",
    "reliabFlag",
]
for var in var_l2:
    ds_sr_match_ppi[f"SR_{var}"] = ds_sr[var]

# Add SR time
ds_sr_match_ppi["SR_time"] = ds_sr_match_ppi["time"]
ds_sr_match_ppi = ds_sr_match_ppi.drop("time")

# Mask SR beams not matching the GR PPI
ds_sr_match_ppi = ds_sr_match_ppi.where(mask_matched_ppi_2d)

# Remove unecessary coordinates
unecessary_coords = [
    "radar_frequency",
    "sweep_mode",
    "prt_mode",
    "follow_mode",
    "latitude",
    "longitude",
    "altitude",
    "crs_wkt",
    "time",
    "crsWGS84",
    "dataQuality",
    "sunLocalTime",
    "SCorientation",
]
for coord in unecessary_coords:
    if coord in ds_sr_match_ppi:
        ds_sr_match_ppi = ds_sr_match_ppi.drop(coord)
    if coord in mask_matched_ppi_2d.coords:
        mask_matched_ppi_2d = mask_matched_ppi_2d.drop(coord)

# Stack aggregated dataset to beam dimension index
ds_sr_stack = ds_sr_match_ppi.stack(sr_beam_index=("along_track", "cross_track"))
da_mask_matched_ppi_stack = mask_matched_ppi_2d.stack(sr_beam_index=("along_track", "cross_track"))

# Drop beams not matching the GR PPI
ds_sr_match = ds_sr_stack.isel(sr_beam_index=da_mask_matched_ppi_stack)

Display aggregated SR variables#

[ ]:

# Original SR near surface reflectivity
p = ds_sr["zFactorFinalNearSurface"].gpm.plot_map(vmin=10, vmax=40)
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

# Original GR reflectivities
p = ds_gr["DBZH"].where(mask_gr).xradar_dev.plot_map(vmin=10, vmax=40)
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))
ds_sr.gpm.plot_swath_lines(ax=p.axes)

<matplotlib.lines.Line2D at 0x7fe79d54fdd0>

../_images/tutorials_tutorial_03_SR_GR_Matching_110_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_110_2.png

[ ]:

# SR reflectivity statistics matching the GR PPI (using xr.Dataset)
p = ds_sr_ppi_max["zFactorFinal_Ku"].gpm.plot_map(cmap=cmap, vmin=10, vmax=40)
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

p = ds_sr_ppi_mean["zFactorFinal_Ku"].gpm.plot_map(cmap=cmap, vmin=10, vmax=40)
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

p = ds_sr_ppi_mean["zFactorMeasured_Ku"].gpm.plot_map(cmap=cmap, vmin=10, vmax=40)
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

../_images/tutorials_tutorial_03_SR_GR_Matching_111_0.png

../_images/tutorials_tutorial_03_SR_GR_Matching_111_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_111_2.png

[ ]:

# Show difference between using Ku mean and max !
diff = ds_sr_ppi_max["zFactorFinal_Ku"] - ds_sr_ppi_mean["zFactorFinal_Ku"]
diff.name = "Difference between Z max and Z mean"
p = diff.gpm.plot_map(cmap="Spectral_r")
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

../_images/tutorials_tutorial_03_SR_GR_Matching_112_0.png

[ ]:

# Show SR mean precip rate
p = ds_sr_ppi_mean["precipRate"].gpm.plot_map()
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

../_images/tutorials_tutorial_03_SR_GR_Matching_113_0.png

[ ]:

# Display reflectivity in 3D

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(
    ds_sr_match["SR_x_mean"].values,
    ds_sr_match["SR_y_mean"].values,
    ds_sr_match["SR_z_mean"].values,
    c=ds_sr_match["SR_zFactorFinal_Ku_mean"].values,
    cmap="turbo",
    vmin=10,
    vmax=40,
    s=5,
)

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fe7d8c501d0>

../_images/tutorials_tutorial_03_SR_GR_Matching_114_1.png

[ ]:

# Disply SR footprint horizontal resolution
p = ds_sr_ppi_max["hres"].gpm.plot_map(cmap="turbo")
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_115_0.png

[ ]:

# Display SR depth over which data are integrated
p = ds_sr_ppi_sum["vres"].gpm.plot_map(cmap="turbo")
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_116_0.png

[ ]:

# Display SR volume over which data are integrated
p = ds_sr_ppi_sum["gate_volume"].gpm.plot_map(cmap="turbo")
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_117_0.png

[ ]:

# Display GR radar volume
p = vol_gr.xradar_dev.plot_map()
p.axes.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))
plt.show()

../_images/tutorials_tutorial_03_SR_GR_Matching_118_0.png

Retrieve the SR footprints polygons#

Here we retrieve SR footprints polygons and then we construct a geopandas.DataFrame with the aggregated SR reflectivities.

[ ]:

# Retrieve SR footprint polygons (using the footprint radius in AEQD x,y coordinates)
xy_mean_sr = np.stack([ds_sr_match["SR_x_mean"], ds_sr_match["SR_y_mean"]], axis=-1)
footprint_radius = ds_sr_match["SR_hres_max"].to_numpy() / 2
sr_poly = [Point(x, y).buffer(footprint_radius[i]) for i, (x, y) in enumerate(xy_mean_sr)]

[ ]:

# Create geopandas DataFrame
df_sr_match = ds_sr_match.to_dataframe()
gdf_sr = gpd.GeoDataFrame(df_sr_match, crs=crs_gr, geometry=sr_poly)
gdf_sr.columns

Index(['SR_counts', 'SR_counts_valid',
       'SR_zFactorFinal_Ku_fraction_above_10dBZ',
       'SR_zFactorFinal_Ku_fraction_above_12dBZ',
       'SR_zFactorFinal_Ku_fraction_above_14dBZ',
       'SR_zFactorFinal_Ku_fraction_above_16dBZ',
       'SR_zFactorFinal_Ku_fraction_above_18dBZ',
       'SR_zFactorFinal_S_fraction_above_10dBZ',
       'SR_zFactorFinal_S_fraction_above_12dBZ',
       'SR_zFactorFinal_S_fraction_above_14dBZ',
       ...
       'lat', 'gpm_id', 'gpm_granule_id', 'gpm_cross_track_id',
       'gpm_along_track_id', 'sweep_number', 'sweep_fixed_angle',
       'along_track', 'cross_track', 'geometry'],
      dtype='object', length=165)

[ ]:

# Extract radar gate polygon on the range-azimuth axis
sr_poly = np.array(gdf_sr.geometry)

Retrieve the GR gates polygons#

Here we retrieve GR gates polygons and then we construct a geopandas.DataFrame containing the GR measurements.

Currently the GR gates polygons are derived inferring the quadmesh around the GR gate centroids (using the x and y GR CRS coordinates).

More accurate GR gates polygons could be derived using the GR beamwidth. We welcome contributions addressing this point !

[ ]:

# Add variables to GR dataset
ds_gr["gate_volume"] = vol_gr
ds_gr["vres"] = v_res_gr
ds_gr["hres"] = h_res_gr

# Add path_integrated_reflectivities
ds_gr["Z_cumsum"] = ds_gr["DBZH"].gpm.idecibel.cumsum("range").gpm.decibel
ds_gr["Z_cumsum"] = ds_gr["Z_cumsum"].where(np.isfinite(ds_gr["Z_cumsum"]))

# Mask reflectivites above minimum GR Z threshold
mask_gr = ds_gr["DBZH"] > z_min_threshold_gr  # important !

# Subset gates in range distance interval
ds_gr_subset = ds_gr.sel(range=slice(min_gr_range_lb, max_gr_range_ub)).where(mask_gr)

# Retrieve geopandas dataframe
gdf_gr = ds_gr_subset.xradar_dev.to_geopandas()  # here we currently infer the quadmesh using the x,y coordinates
display(gdf_gr)

		DBZH	ZDR	PHIDP	RHOHV	CCORH	gate_volume	vres	hres	Z_cumsum	sweep_mode	...	altitude	crs_wkt	x	y	z	lon	lat	range	azimuth	geometry
azimuth	range
0.252686	2125.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	azimuth_surveillance	...	320	0	9.370927	2124.819145	338.193914	-117.041701	32.938177	2125.0	0.252686	POLYGON ((0.096 1999.907, 17.543 1999.754, 19....
	2375.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	azimuth_surveillance	...	320	0	10.473386	2374.797266	340.369314	-117.041689	32.940431	2375.0	0.252686	POLYGON ((0.108 2249.894, 19.736 2249.722, 21....
	2625.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	azimuth_surveillance	...	320	0	11.575845	2624.775259	342.552070	-117.041678	32.942685	2625.0	0.252686	POLYGON ((0.12 2499.882, 21.929 2499.691, 24.1...
	2875.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	azimuth_surveillance	...	320	0	12.678303	2874.753123	344.742182	-117.041666	32.944939	2875.0	0.252686	POLYGON ((0.132 2749.869, 24.122 2749.659, 26....
	3125.0	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	azimuth_surveillance	...	320	0	13.780761	3124.730858	346.939649	-117.041654	32.947193	3125.0	0.252686	POLYGON ((0.144 2999.857, 26.315 2999.627, 28....
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
359.744568	148875.0	32.0	20.0	263.742454	1.051667	-8.0	1.325715e+09	2598.424805	2598.424805	54.805861	azimuth_surveillance	...	320	0	-663.487562	148825.386039	2879.986614	-117.049005	34.260814	148875.0	359.744568	POLYGON ((-1290.197 148695.029, -35.665 148705...
	149125.0	27.0	20.0	291.597606	0.528333	-8.0	1.330171e+09	2602.788330	2602.788330	54.813053	azimuth_surveillance	...	320	0	-664.601337	149075.214475	2886.478208	-117.049017	34.263066	149125.0	359.744568	POLYGON ((-1292.364 148944.848, -35.725 148955...
	149375.0	29.5	20.0	327.562486	1.018333	-8.0	1.334635e+09	2607.151855	2607.151855	54.825813	azimuth_surveillance	...	320	0	-665.715111	149325.042530	2892.977147	-117.049029	34.265318	149375.0	359.744568	POLYGON ((-1294.532 149194.668, -35.785 149205...
	149625.0	27.5	20.0	346.602717	0.845000	-8.0	1.339106e+09	2611.515137	2611.515137	54.833845	azimuth_surveillance	...	320	0	-666.828883	149574.870202	2899.483431	-117.049042	34.267571	149625.0	359.744568	POLYGON ((-1296.699 149444.486, -35.845 149455...
	149875.0	29.5	20.0	306.759271	1.011667	-6.0	1.343585e+09	2615.878662	2615.878662	54.846543	azimuth_surveillance	...	320	0	-667.942653	149824.697491	2905.997060	-117.049054	34.269823	149875.0	359.744568	POLYGON ((-1298.867 149694.305, -35.905 149705...

426240 rows × 28 columns

[ ]:

# Extract radar gate polygon on the range-azimuth axis
gr_poly = np.array(gdf_gr.geometry)

Aggregate GR data on SR footprints#

To aggregate GR radar gates onto SR beam polygon footprints, we will use the PolyAggregator. The PolyAggregator enable custom aggregation of source (GR) polygons data intesecting the target (SR) polygons.

When aggregating the data, you can use a combination of weighting schemes By default, PolyAggregator applies area_weighting, which assigns weights based on the proportion of source polygon area that intersects the target polygon.

Custom Weighting Schemes

Several additional weighting schemes can be utilized, typically derived as a function of the horizontal distance between the centroids of the source and target polygons. Examples include:

Inverse Distance Weighting: \(\text{weights} = \frac{1}{\text{dist}^{\text{order}}}\)
Barnes Gaussian Weighting: \(\text{weights} = \exp\left(-\frac{\text{dist}^2}{\kappa^2}\right)\)
Cressman Weighting: \(\text{weights} = \frac{\text{max\_distance}^2 - \text{dist}^2}{\text{max\_distance}^2 + \text{dist}^2}\)

For example, the Barnes Gaussian Distance Weighting Scheme weights values proportionally to the distance measured from the center of the SR footprint to the center of the GR radar gate. This approach can partially account for the nonuniform distribution of power within the SR beam.

Custom weights can however also be passed to PolyAggregator, which will automatically handle normalization.

In the context of SR/GR matching, it may be beneficial to adopt a Volume Weighting Scheme, using the volumes of GR radar gates as weights so that larger volumes are weighted more heavily.

In the code here below, we currently apply only the Area Weighting Scheme:

[ ]:

# Define PolyAggregator
aggregator = PolyAggregator(source_polygons=gr_poly, target_polygons=sr_poly)

[ ]:

# Aggregate GR reflecitvities and compute statistics
# - Timestep of acquisition
time_gr = aggregator.first(values=gdf_gr["time"])

# - Total number of gates aggregated
counts = aggregator.counts()
counts_valid = aggregator.apply(lambda x, weights: np.sum(~np.isnan(x)), values=gdf_gr["DBZH"])  # noqa

# - Total gate volume
sum_vol = aggregator.sum(values=gdf_gr["gate_volume"])

# - Fraction of SR area covered
fraction_covered_area = aggregator.fraction_covered_area()

# - Reflectivity statistics
z_mean = wrl.trafo.decibel(aggregator.average(values=wrl.trafo.idecibel(gdf_gr["DBZH"])))
z_std = wrl.trafo.decibel(aggregator.std(values=wrl.trafo.idecibel(gdf_gr["DBZH"])))
z_std[np.isinf(z_std)] = 0  # If only 1 value, std=0, log transform become -inf --> Set to 0
z_max = aggregator.max(values=gdf_gr["DBZH"])
z_min = aggregator.min(values=gdf_gr["DBZH"])
z_range = z_max - z_min
z_cov = z_std / z_mean  # coefficient of variation

Now we create a geopandas.DataFrame with aggregated GR statistics:

[ ]:

# Create DataFrame with GR matched statistics
df = pd.DataFrame(
    {
        "GR_Z_mean": z_mean,
        "GR_Z_std": z_std,
        "GR_Z_max": z_max,
        "GR_Z_min": z_min,
        "GR_Z_range": z_range,
        "GR_Z_cov": z_cov,
        "GR_time": time_gr,
        "GR_gate_volume_sum": sum_vol,
        "GR_fraction_covered_area": fraction_covered_area,
        "GR_counts": counts,
        "GR_counts_valid": counts_valid,
    },
    index=gdf_sr.index,
)
gdf_gr_match = gpd.GeoDataFrame(df, crs=crs_gr, geometry=aggregator.target_polygons)
gdf_gr_match.head()

# Add GR range statistics
gdf_gr_match["GR_range_max"] = aggregator.max(values=gdf_gr["range"])
gdf_gr_match["GR_range_mean"] = aggregator.mean(values=gdf_gr["range"])
gdf_gr_match["GR_range_min"] = aggregator.min(values=gdf_gr["range"])

# Fraction above sensitivity thresholds
gr_sensitivity_thresholds = [8, 10, 12, 14, 16, 18]
for thr in sr_sensitivity_thresholds:
    fraction = aggregator.apply(lambda x, weights: np.sum(x > thr), values=gdf_gr["DBZH"]) / counts  # noqa
    gdf_gr_match[f"GR_Z_fraction_above_{thr}dBZ"] = fraction

# Compute further aggregation statistics
stats_var = ["vres", "hres", "x", "y", "z", "Z_cumsum"]
for var in stats_var:
    gdf_gr_match[f"GR_{var}_mean"] = aggregator.average(values=gdf_gr[var])
    gdf_gr_match[f"GR_{var}_min"] = aggregator.min(values=gdf_gr[var])
    gdf_gr_match[f"GR_{var}_max"] = aggregator.max(values=gdf_gr[var])
    gdf_gr_match[f"GR_{var}_std"] = aggregator.std(values=gdf_gr[var])
    gdf_gr_match[f"GR_{var}_range"] = gdf_gr_match[f"GR_{var}_max"] - gdf_gr_match[f"GR_{var}_min"]

# Compute horizontal distance between centroids
func_dict = {"min": np.min, "mean": np.mean, "max": np.max}
for suffix, func in func_dict.items():
    arr = np.zeros(aggregator.n_target_polygons) * np.nan
    arr[aggregator.target_intersecting_indices] = np.array(
        [func(dist) for i, dist in aggregator.dict_distances.items()],
    )
    gdf_gr_match[f"distance_horizontal_{suffix}"] = arr

To compute custom statistics, you can use the following approaches. Remember that the function must have the values and weights arguments. Add the weights argument also if they are not used in your function !

[ ]:

def weighted_mean(values, weights):
    return np.average(values, weights=weights)


aggregated_values = aggregator.apply(weighted_mean, values=gdf_gr["DBZH"], area_weighting=True)
aggregated_values = aggregator.apply(
    lambda x, weights: np.average(x, weights=weights),
    values=gdf_gr["DBZH"],
    area_weighting=True,
)

Here below we compare the usage of Inverse Distance Weighting, Barnes Gaussian Distance Weighting, Cressman Weighting, and Volume Weighting.

[ ]:

from gpm.utils.zonal_stats import BarnesGaussianWeights, CressmanWeights, InverseDistanceWeights, PolyAggregator

gdf_gr_match1 = gdf_gr_match.copy()

# Weighted by intersection area
gdf_gr_match1["GR_Z_area"] = aggregator.mean(values=gdf_gr["DBZH"], area_weighting=True)

# Weighted by ntersection area and Inverse Distance Weighting
gdf_gr_match1["GR_Z_idw"] = aggregator.mean(
    values=gdf_gr["DBZH"],
    area_weighting=True,
    distance_weighting=InverseDistanceWeights(order=1),
)
# Weighted by intersection area and Cressman Weighting
gdf_gr_match1["GR_Z_cressman"] = aggregator.mean(
    values=gdf_gr["DBZH"],
    area_weighting=True,
    distance_weighting=CressmanWeights(max_distance=2000),
)
# Weighted by intersection area and Barnes Gaussian Weighting (kappa=SR_hres_mean)
gdf_gr_match1["GR_Z_barnes"] = aggregator.mean(
    values=gdf_gr["DBZH"],
    area_weighting=True,
    distance_weighting=BarnesGaussianWeights(kappa=gdf_sr["SR_hres_mean"]),
)
# Weighted Barnes Gaussian Weighting (kappa=SR_radius_mean)
gdf_gr_match1["GR_Z_barnes_without_area"] = aggregator.mean(
    values=gdf_gr["DBZH"],
    area_weighting=False,
    distance_weighting=BarnesGaussianWeights(kappa=gdf_sr["SR_hres_mean"]),
)
# Weighted by intersection area and by GR volume
gdf_gr_match1["GR_Z_volume"] = aggregator.mean(
    values=gdf_gr["DBZH"],
    weights=gdf_gr["gate_volume"],
    area_weighting=True,
)

[ ]:

gdf_gr_match1.plot(column="GR_Z_area", legend=True, cmap="Spectral_r", vmin=0, vmax=40)
gdf_gr_match1.plot(column="GR_Z_idw", legend=True, cmap="Spectral_r", vmin=0, vmax=40)
gdf_gr_match1.plot(column="GR_Z_cressman", legend=True, cmap="Spectral_r", vmin=0, vmax=40)
gdf_gr_match1.plot(column="GR_Z_barnes", legend=True, cmap="Spectral_r", vmin=0, vmax=40)
gdf_gr_match1.plot(column="GR_Z_volume", legend=True, cmap="Spectral_r", vmin=0, vmax=40)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_137_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_137_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_137_3.png

../_images/tutorials_tutorial_03_SR_GR_Matching_137_4.png

../_images/tutorials_tutorial_03_SR_GR_Matching_137_5.png

[ ]:

# Impact of Barnes with and without area-weighting
gdf_gr_match1["DIFF_Barnes_Without_Area"] = gdf_gr_match1["GR_Z_barnes"] - gdf_gr_match1["GR_Z_barnes_without_area"]
gdf_gr_match1.plot(column="DIFF_Barnes_Without_Area", legend=True, cmap="Spectral_r", vmin=-0.25, vmax=0.25)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_138_1.png

[ ]:

# Impact of Barnes with area-weighting and default area weighting
gdf_gr_match1["DIFF_Barnes_Naive"] = gdf_gr_match1["GR_Z_area"] - gdf_gr_match1["GR_Z_barnes"]
gdf_gr_match1.plot(column="DIFF_Barnes_Naive", legend=True, cmap="Spectral_r", vmin=-0.25, vmax=0.25)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_139_1.png

Create the SR/GR Database#

Let’s back to our main task and merge SR/GR aggregated variables into an unique database:

[ ]:

# Create DataFrame with matched variables
gdf_match = gdf_gr_match.merge(gdf_sr, on="geometry")

# Compute ratio SR/GR volume
gdf_match["VolumeRatio"] = gdf_match["SR_gate_volume_sum"] / gdf_match["GR_gate_volume_sum"]

# Compute difference in SR/GR gate volume
gdf_match["VolumeDiff"] = gdf_match["SR_gate_volume_sum"] - gdf_match["GR_gate_volume_sum"]

# Compute time difference [in seconds]
gdf_match["time_difference"] = gdf_match["GR_time"] - gdf_match["SR_time"]
gdf_match["time_difference"] = gdf_match["time_difference"].dt.total_seconds().astype(int)

# Compute lower bound and upper bound height
gdf_match["GR_z_lower_bound"] = gdf_match["GR_z_min"] - gdf_match["GR_vres_mean"] / 2
gdf_match["GR_z_upper_bound"] = gdf_match["GR_z_max"] + gdf_match["GR_vres_mean"] / 2

gdf_match["SR_z_lower_bound"] = gdf_match["SR_z_min"] - gdf_match["SR_vres_mean"] / 2
gdf_match["SR_z_upper_bound"] = gdf_match["SR_z_max"] + gdf_match["SR_vres_mean"] / 2

Analyse the SR/GR Database#

[ ]:

# Show fraction SR footprint horizontal area covered by GR
gdf_match.plot(column="GR_fraction_covered_area", legend=True, cmap="Spectral_r")

# Show ratio SR/GR volume
gdf_match.plot(column="VolumeRatio", legend=True, cmap="Spectral_r")

# Show difference in total gate volume
gdf_match.plot(column="VolumeDiff", legend=True, cmap="Spectral_r")

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_143_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_143_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_143_3.png

[ ]:

gdf_gr_match.plot(column="distance_horizontal_max", legend=True, cmap="Spectral_r")
gdf_gr_match.plot(column="distance_horizontal_mean", legend=True, cmap="Spectral_r")
gdf_gr_match.plot(column="distance_horizontal_min", legend=True, cmap="Spectral_r")

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_144_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_144_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_144_3.png

[ ]:

# Analyse variations in gate resolution
gdf_match.plot(column="GR_vres_mean", legend=True, cmap="Spectral_r")
gdf_match.plot(column="SR_vres_mean", legend=True, cmap="Spectral_r")

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_145_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_145_2.png

[ ]:

# Analyse variations in minimum and maximum gate elevation
gdf_match.plot(column="GR_z_upper_bound", legend=True, cmap="Spectral_r")
gdf_match.plot(column="SR_z_upper_bound", legend=True, cmap="Spectral_r")

gdf_match.plot(column="GR_z_lower_bound", legend=True, cmap="Spectral_r")
gdf_match.plot(column="SR_z_lower_bound", legend=True, cmap="Spectral_r")

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_146_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_146_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_146_3.png

../_images/tutorials_tutorial_03_SR_GR_Matching_146_4.png

[ ]:

# Analyse fraction hydrometeors and sensitivity
gdf_match.plot(column="SR_counts", legend=True, cmap="Spectral_r", vmin=0)
gdf_match.plot(column="SR_counts_valid", legend=True, cmap="Spectral_r", vmin=0)

gdf_match.plot(column="SR_fraction_clutter", legend=True, cmap="Spectral_r", vmin=0, vmax=1)

# gdf_match.plot(column="SR_fraction_melting_layer", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_fraction_hail", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_fraction_snow", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_fraction_rain", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_fraction_no_precip", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_fraction_below_isotherm", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_fraction_above_isotherm", legend=True, cmap="Spectral_r", vmin=0, vmax=1)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_147_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_147_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_147_3.png

[ ]:

# Analyse sensitivity
gdf_match.plot(column="GR_Z_fraction_above_10dBZ", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="GR_Z_fraction_above_12dBZ", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="GR_Z_fraction_above_14dBZ", legend=True, cmap="Spectral_r", vmin=0, vmax=1)

gdf_match.plot(column="SR_zFactorFinal_Ku_fraction_above_12dBZ", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_zFactorFinal_Ku_fraction_above_14dBZ", legend=True, cmap="Spectral_r", vmin=0, vmax=1)
# gdf_match.plot(column="SR_zFactorFinal_Ku_fraction_above_16dBZ", legend=True, cmap="Spectral_r", vmin=0, vmax=1)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_148_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_148_2.png

[ ]:

# Analyse SR attenuation statistics
gdf_match.plot(column="SR_reliabFlag", legend=True, cmap="Spectral_r")
gdf_match.plot(column="SR_pathAtten", legend=True, cmap="Spectral_r", vmin=0, vmax=5)
gdf_match.plot(column="SR_piaFinal", legend=True, cmap="Spectral_r", vmin=0, vmax=5)
gdf_match.plot(column="SR_zFactorCorrection_Ku_max", legend=True, cmap="Spectral_r", vmin=0, vmax=5)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_149_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_149_2.png

../_images/tutorials_tutorial_03_SR_GR_Matching_149_3.png

../_images/tutorials_tutorial_03_SR_GR_Matching_149_4.png

[ ]:

# Analyse SR NUBF
gdf_match.plot(column="SR_zFactorFinal_Ku_std", legend=True, cmap="Spectral_r", vmin=0, vmax=8)
gdf_match.plot(column="GR_Z_std", legend=True, cmap="Spectral_r", vmin=0, vmax=8)

# gdf_match.plot(column="SR_zFactorFinal_Ku_range", legend=True, cmap="Spectral_r", vmin=0, vmax=8)
# gdf_match.plot(column="GR_Z_range", legend=True, cmap="Spectral_r", vmin=0, vmax=20)

# gdf_match.plot(column="SR_zFactorFinal_Ku_cov", legend=True, cmap="Spectral_r")
# gdf_match.plot(column="GR_Z_cov", legend=True, cmap="Spectral_r")

# gdf_match["SR_zFactorFinal_Ku_range_rel"] = gdf_match["SR_zFactorFinal_Ku_range"]/gdf_match["SR_zFactorFinal_Ku_mean"]
# gdf_match.plot(column="SR_zFactorFinal_Ku_range_rel", legend=True, cmap="Spectral_r", vmin=0, vmax=1)

# gdf_match["GR_Z_range_rel"] = gdf_match["GR_Z_range"]/gdf_match["GR_Z_mean"]
# gdf_match.plot(column="GR_Z_range_rel", legend=True, cmap="Spectral_r")

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_150_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_150_2.png

[ ]:

# Compare reflectivities
extent_xy = gdf_sr.total_bounds[[0, 2, 1, 3]]
p = gdf_match.plot(
    column="SR_zFactorFinal_Ku_mean",
    legend=True,
    legend_kwds={"label": "SR Reflectivity (dBz)"},
    cmap="Spectral_r",
    vmin=10,
    vmax=40,
)
p.axes.set_xlim(extent_xy[0:2])
p.axes.set_ylim(extent_xy[2:4])
p.axes.set_title("SR")
plt.xlabel("x (m)", fontsize=12)
plt.ylabel("y (m)", fontsize=12)
plt.grid(lw=0.25, color="grey")

p = gdf_match.plot(
    column="GR_Z_mean",
    legend=True,
    legend_kwds={"label": "GR Reflectivity (dBz)"},
    cmap="Spectral_r",
    vmin=10,
    vmax=40,
)
p.axes.set_xlim(extent_xy[0:2])
p.axes.set_ylim(extent_xy[2:4])
p.axes.set_title("GR")
plt.xlabel("x (m)", fontsize=12)
plt.ylabel("y (m)", fontsize=12)
plt.grid(lw=0.25, color="grey")

../_images/tutorials_tutorial_03_SR_GR_Matching_151_0.png

../_images/tutorials_tutorial_03_SR_GR_Matching_151_1.png

[ ]:

# Compute difference (also with bias removed)
gdf_match["DIFF_Z_SR_GR"] = gdf_match["SR_zFactorFinal_Ku_mean"] - gdf_match["GR_Z_mean"]
gdf_match["DIFF_Z_GR_SR"] = gdf_match["GR_Z_mean"] - gdf_match["SR_zFactorFinal_Ku_mean"]
bias_gr_sr = np.nanmedian(gdf_match["DIFF_Z_GR_SR"]).round(2)
bias_sr_gr = np.nanmedian(gdf_match["DIFF_Z_SR_GR"]).round(2)
p = gdf_match.plot(
    column="DIFF_Z_SR_GR",
    legend_kwds={"label": "SR-GR Reflectivity Difference (dBz)"},
    cmap="RdBu",
    legend=True,
    vmax=10,
)
p.axes.set_xlim(extent_xy[0:2])
p.axes.set_ylim(extent_xy[2:4])
p.axes.set_title("SR - GR")
plt.xlabel("x (m)", fontsize=12)
plt.ylabel("y (m)", fontsize=12)
plt.grid(lw=0.25, color="grey")

gdf_match["DIFF_Z_GR_SR_unbiased"] = gdf_match["DIFF_Z_GR_SR"] - bias_gr_sr
p = gdf_match.plot(
    column="DIFF_Z_GR_SR_unbiased",
    legend_kwds={"label": "SR-GR Reflectivity Difference (wihout bias) (dBz)"},
    cmap="RdBu",
    legend=True,
    vmin=-4,
    vmax=4,
)
p.axes.set_xlim(extent_xy[0:2])
p.axes.set_ylim(extent_xy[2:4])
p.axes.set_title("SR - GR")
plt.xlabel("x (m)", fontsize=12)
plt.ylabel("y (m)", fontsize=12)
plt.grid(lw=0.25, color="grey")

../_images/tutorials_tutorial_03_SR_GR_Matching_152_0.png

../_images/tutorials_tutorial_03_SR_GR_Matching_152_1.png

Here we show how to explore interactively the reflectivity fields using Folium:

[ ]:

gdf_match1 = gdf_match.copy()
gdf_match1.explore(column="GR_Z_mean", legend=True, cmap="Spectral_r", vmin=0, vmax=40)

Make this Notebook Trusted to load map: File -> Trust Notebook

Here we show how to plot cartopy maps in WGS84 and GR CRS AEQD.

[ ]:

# Define Cartopy AEQD CRS
ccrs_aeqd = ccrs.AzimuthalEquidistant(
    central_longitude=ds_gr["longitude"].item(),
    central_latitude=ds_gr["latitude"].item(),
)

[ ]:

# Define Cartopy WGS84 CRS
ccrs_wgs84 = ccrs.PlateCarree()

[ ]:

# Print CRS of geopandas geometries
gdf_match.crs

<Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Wor ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Modified Azimuthal Equidistant
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

[ ]:

# Plot in AEQD CRS (same crs of geopandas geometries)
fig, ax = plt.subplots(subplot_kw={"projection": ccrs_aeqd})
gdf_match.plot(ax=ax, column="GR_Z_mean", cmap="turbo", legend=True)

<GeoAxes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_159_1.png

[ ]:

# Convert geopandas dataframe to another crs
gdf_match_wgs84 = gdf_match.to_crs(crs_sr)
# Plot in WGS84 CRS
fig, ax = plt.subplots(subplot_kw={"projection": ccrs_wgs84})
# plot_cartopy_background(ax)
gdf_match_wgs84.plot(ax=ax, column="GR_Z_mean", cmap="turbo", legend=True)
ax.set_extent(ds_gr.xradar_dev.extent(max_distance=max_gr_range))

../_images/tutorials_tutorial_03_SR_GR_Matching_160_0.png

Filter the SR/GR Database#

When calibrating or comparing SR and GR data, there is a necessary tradeoff between the strictness of filtering criteria filtering and the number of available samples. Here below we provide some general recommendations for effective filtering:

Sensitivity Thresholds: Retain only those radar beams where the aggregated gate reflectivities exceed the instrument sensitivities. The GR_Z_fraction_above_<thr>dBZ and SR_zFactorFinal_Ku_fraction_above_<thr>dBZ variables can be used filter the samples.
Stratiform Precipitation: If the purpose of your analysis is to assess the calibration bias of GR data, it is suggested to focus on stratiform precipitation that ocurrs outside the melting layer (i.e. avoiding the bright band). Only stratiform samples above the melting layer are used in ground validation of the GPM DPR (W. Petersen 2017, personal communication). Convective SR footprints are typically excluded to avoid dealing with:
- the high spatial variability in the precipitation field and issues with non-uniform beam filling (NUBF),
- the potential biases introduced by the SR attenuation correction,
- the potential biases in GR reflectivities, especially at C and X band, due to beam attenuation,
- the multiple scattering signature caused by hail particles.
Reflectivity Range: According to Warren et al. (2018), volume-averaged SR and GR reflectivity values should be selected within the range of 24 to 36 dBZ. This range minimizes the impact of low SR sensitivity, SR beam attenuation and non-rayleigh scattering effects.
Clutter Removal: Exclude SR/GR samples that are contaminated by ground clutter, anomalous propagation, or beam blockage.
Volume matching: Exclude SR/GR samples where there are excessive differences in the total gate volume.

[ ]:

# Define mask
masks = [
    # Select SR scan with "normal" dataQuality (for entire cross-track scan)
    gdf_match["SR_dataQuality"] == 0,
    # Select SR beams with detected precipitation
    gdf_match["SR_flagPrecip"] > 0,
    # Select only 'high quality' SR data
    # - qualityFlag == 1 indicates low quality retrievals
    # - qualityFlag == 2 indicates bad/missing retrievals
    gdf_match["SR_qualityFlag"] == 0,
    # Select SR beams with SR above minimum reflectivity
    gdf_match["SR_zFactorFinal_Ku_fraction_above_12dBZ"] > 0.95,
    # Select SR beams with GR above minimum reflectivity
    gdf_match["GR_Z_fraction_above_10dBZ"] > 0.95,
    # Select only SR beams with detected bright band
    gdf_match["SR_qualityBB"] == 1,
    # Select only beams with confident precipitation type
    gdf_match["SR_qualityTypePrecip"] == 1,
    # Select only stratiform precipitation
    # - SR_flagPrecipitationType == 2 indicates convective
    gdf_match["SR_flagPrecipitationType"] == 1,
    # Select only SR beams with reliable attenuation correction
    gdf_match["SR_reliabFlag"].isin((1, 2)),  # or == 1
    # Select only beams with reduced path attenuation
    # gdf_match["SR_zFactorCorrection_Ku_max"]
    # gdf_match["SR_piaFinal"]
    # gdf_match["SR_pathAtten"]
    # Select only SR beams with matching SR gates with no clutter
    gdf_match["SR_fraction_clutter"] == 0,
    # Select only SR beams with matching SR gates not in the melting layer
    gdf_match["SR_fraction_melting_layer"] == 0,
    # Select only SR beams with matching SR gates with precipitation
    gdf_match["SR_fraction_no_precip"] == 0,
    # Select only SR beams with matching SR gates with no hail
    # gdf_match["SR_fraction_hail"] == 0,
    # Select only SR beams with matching SR gates with rain
    # gdf_match["SR_fraction_hail"] == 1,
    # Select only SR beams with matching SR gates with snow
    # gdf_match["SR_fraction_snow"] == 1,
    # Discard SR beams with high NUBF
    # gdf_match["SR_zFactorFinal_Ku_range"] > 5,
    # gdf_match["SR_zFactorFinal_Ku_cov"] < 0.5,
    # gdf_match["GR_Z_cov"] < 0.5,
    # gdf_match["GR_Z_range"] > 5,
    # Select only interval of reflectivities
    # - Warren et al., 2018:  between 24 and 36 dBZ
    # (gdf_match["SR_zFactorFinal_Ku_mean"] > 24) & (gdf_match["SR_zFactorFinal_Ku_mean"] < 36),
    # (gdf_match["GR_Z_mean"] > 24) & (gdf_match["GR_Z_mean"] < 36),
    # Select SR beams only within given GR radius interval
    # (gdf_match["GR_range_min"] > min_gr_range) & (gdf_match["GR_range_max"] < max_gr_range),
    # Discard SR beams where scanning time difference (provided in seconds)
    gdf_match["time_difference"] < 5 * 60,  # 5 min
    # Discard SR beams where GR gates does not cover 80% of the horizontal area
    # gdf_match["GR_fraction_covered_area"] > 0.8,
    # Filter footprints where volume ratio exceeds 60
    # gdf_match["VolumeRatio"] > 60,
    # Filter footprints where GR affected by beam blockage
    # TODO: Filtering by beam blockage
]

[ ]:

# Define final mask
mask_final = reduce(np.logical_and, masks)

[ ]:

# Dsplay final filtering mask
gdf_match["filtering_mask"] = mask_final
gdf_match.plot(column="filtering_mask", legend=True, cmap="Spectral")
gdf_match.plot(column="SR_zFactorFinal_Ku_mean", legend=True, cmap="Spectral")

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_164_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_164_2.png

[ ]:

# Filter matching dataset
gdf_final = gdf_match[mask_final]

[ ]:

# Plot results
gdf_final.plot(column="SR_zFactorFinal_Ku_mean", legend=True, cmap="Spectral_r", vmin=0, vmax=40)
gdf_final.plot(column="GR_Z_mean", legend=True, cmap="Spectral_r", vmin=0, vmax=40)

<Axes: >

../_images/tutorials_tutorial_03_SR_GR_Matching_166_1.png

../_images/tutorials_tutorial_03_SR_GR_Matching_166_2.png

Determine the GR Calibration Bias#

[ ]:

# Compute statistics
bias = np.nanmean(gdf_match["DIFF_Z_GR_SR"]).round(2)
robBias = np.nanmedian(gdf_match["DIFF_Z_GR_SR"]).round(2)

# R2 (pearson, spearman rank)
print(f"Th ZH Calibration offset is (mean): {bias} dB")
print(f"Th ZH Calibration offset is (median): {robBias} dB")

Th ZH Calibration offset is (mean): -5.14 dB
Th ZH Calibration offset is (median): -4.36 dB

[ ]:

# Comparison of GR vs SR reflectivities
hue_variable = "VolumeDiff"
hue_title = "Difference in volume between SR and GR [m³]"

[ ]:

# - Histograms
fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(121, aspect="equal")
plt.scatter(
    gdf_match["GR_Z_mean"],
    gdf_match["SR_zFactorFinal_Ku_mean"],
    marker="+",
    c=gdf_match[hue_variable],
    cmap="turbo",
)
plt.colorbar(label=hue_title)
plt.plot([0, 60], [0, 60], linestyle="solid", color="black")
plt.xlim(10, 50)
plt.ylim(10, 50)
plt.xlabel("GR reflectivity (dBZ)")
plt.ylabel("SR reflectivity (dBZ)")
plt.title(f"Offset GR-SR: {bias} dBZ")
ax = fig.add_subplot(122)
plt.hist(
    gdf_match["GR_Z_mean"][gdf_match["GR_Z_mean"] > z_min_threshold_gr],
    bins=np.arange(10, 60, 5),
    edgecolor="None",
    label="GR",
)
plt.hist(
    gdf_match["SR_zFactorFinal_Ku_mean"][gdf_match["SR_zFactorFinal_Ku_mean"] > z_min_threshold_sr],
    bins=np.arange(10, 60, 5),
    edgecolor="red",
    facecolor="None",
    label="SR",
)
plt.xlabel("Reflectivity (dBZ)")
plt.legend()
fig.suptitle("GR vs SR")

Text(0.5, 0.98, 'GR vs SR')

../_images/tutorials_tutorial_03_SR_GR_Matching_170_1.png

Next Steps#

Congratulations on reaching the end of this tutorial! If you have suggestions for improving this tutorial or the SR/GR matching process, please feel free to open a Pull Request or create a GitHub issue on the GPM-API GitHub repository.

The code for the SR/GR matching process is encapsulated in the volume_matching function contained in the gpm.gv module . User are free to use or adapt such function to accomplish their tasks. The Spaceborne-Ground Radar Calibration Applied Tutorial shows how to assess the GR calibration bias by abstracting away the complexities discussed in this tutorial.

References#

Citation#

This notebook is part of the GPM-API documentation.

Large portions of this tutorials were adapted and derived from an old \(\omega radlib\) tutorial.

Copyright: \(\omega radlib\) and GPM-API developers. Distributed under the MIT License. See GPM-API license for more info.

Spaceborne-Ground Radar Matching (Methodology)

Contents