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
from xarray.backends.api import open_datatree
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)

Load SR and GR data¶

In this tutorial, we exploit a coincident GPM DPR and NEXRAD overpass over San Diego, USA, during the landfall of Tropical Storm Hilary in 2023, to illustrate step-by-step the SR/GR volume matching procedure.

To facilitate your hands-on experience, we’ve preprocessed the native NEXRAD radar data using the xradar package and the coincident GPM DPR overpass using GPM-API. The resulting data has been uploaded in Zarr format to the Pythia Cloud Bucket, allowing you to start the tutorial directly in the Binder environment.

If you wish to adapt this tutorial to your own case study, you can use the xradar package to read your GR radar data into the expected xarray format. Similarly, you can retrieve the relevant GPM data using GPM-API. The GPM-API Spaceborne-Ground Radar Matching Tutorial shows how to read NEXRAD GR data directly from the AWS S3 bucket and how to download the GPM data of interest using GPM-API.

To start this tutorial let’s open the NEXRAD ground radar data and the GPM granule with the coincident overpass.

# Open fsspec connection to the bucket
bucket_url = "https://js2.jetstream-cloud.org:8001"
fs = fsspec.filesystem("s3", anon=True, client_kwargs=dict(endpoint_url=bucket_url))

# Open the GR datatree
file = fs.get_mapper(f"pythia/radar/erad2024/gpm_api/KNKX20230820_221341_V06.zarr")
dt_gr = open_datatree(file, engine="zarr", consolidated=True, chunks={})
display(dt_gr)

# Open the SR dataset
file = fs.get_mapper(
    f"pythia/radar/erad2024/gpm_api/2A.GPM.DPR.V9-20211125.20230820-S213941-E231213.053847.V07B.zarr"
)
ds_sr = xr.open_zarr(file, consolidated=True, chunks={})

During this tutorial we will use only the GPM DPR Ku band reflectivities, therefore we subset the dataset here to facilitate the analysis.

Please note that the GPM DPR reflectivites measurements are provided by the ZFactorFinaland ZFactorMeasured variables:

The ZFactorMeasured variable contains the reflectivity measurements with sidelobe clutter removed. However, it still includes ground clutter, background noise, and is not corrected for atmospheric and hydrometeors attenuation.
The ZFactorFinal represents the reflectivity measurements with background noise, ground clutter, and sidelobe clutter removed. The reflectivites have been corrected atmospheric and hydrometeors attenuation.

ds_sr = ds_sr.sel({"radar_frequency": "Ku"})

display(ds_sr)

Here, we define key parameters and settings necessary for the SR/GR comparison. Minimum reflectivity thresholds are applied to exclude certain radar gates from the analysis, primarily to distinguish precipitating regions.

We strongly recommend specifying an “optimistic” minimum reflectivity sensitivity threshold, as reflectivities may be biased, and this approach allows for the assessment of non-uniform beam filling (NUBF) during the SR/GR matching process.

More restrictive filtering can be applied after the SR/GR gate matching is completed.

# Define half power beam width angle (in degrees)
azimuth_beamwidth_gr = 1.0
elevation_beamwidth_gr = 1.0

# Define Z min threshold
z_min_threshold_gr = 0  # to remove dry case
z_min_threshold_sr = 10

# GR matching ROI
min_gr_range_lb = 0
max_gr_range_ub = 150_000

# Define SR CRS
crs_sr = ds_sr.gpm.pyproj_crs

Preprocess GR data with xradar¶

Now we select the GR sweep, we georeference the dataset, we set the GR CRS and we extract GR radar gates coordinates.

sweep_idx = 0
sweep_group = f"sweep_{sweep_idx}"  # GR sweep (elevation to be used)

# Select sweep dataset
ds_gr = dt_gr[sweep_group].to_dataset().compute()

print("Sweep Mode:", ds_gr["sweep_mode"].to_numpy())
print("Sweep Elevation Angle:", ds_gr["sweep_fixed_angle"].to_numpy())

# Set auxiliary info as coordinates
# - To avoid broadcasting i.e. when masking
possible_coords = ["sweep_mode", "sweep_number", "prt_mode", "follow_mode", "sweep_fixed_angle"]
extra_coords = [coord for coord in possible_coords if coord in ds_gr]
ds_gr = ds_gr.set_coords(extra_coords)

Sweep Mode: azimuth_surveillance
Sweep Elevation Angle: 0.4833984375

# Georeference the data on a azimuthal_equidistant projection centered on the radar
ds_gr = ds_gr.xradar.georeference()
ds_gr["crs_wkt"].attrs

{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Azimuthal Equidistant",ID["EPSG",1125]],PARAMETER["Latitude of natural origin",32.919017791748,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-117.041801452637,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.314245179,
 'inverse_flattening': 298.257223563,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'unknown',
 'horizontal_datum_name': 'World Geodetic System 1984',
 'projected_crs_name': 'unknown',
 'grid_mapping_name': 'azimuthal_equidistant',
 'latitude_of_projection_origin': 32.91901779174805,
 'longitude_of_projection_origin': -117.04180145263672,
 'false_easting': 0.0,
 'false_northing': 0.0}

# Get the GR CRS
crs_gr = ds_gr.xradar_dev.pyproj_crs

# Add lon/lat coordinates to GR
# - Useful for plotting ;)
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)
ds_gr["lon"] = lon_gr
ds_gr["lat"] = lat_gr
ds_gr = ds_gr.set_coords(["lon", "lat"])

# Set GR gates with Z < 0 to NaN
# - Following Morris and Schwaller 2011 recommendation
ds_gr = ds_gr.where(ds_gr["DBZH"] >= 0)

# Show some GR fields
ds_gr["DBZH"].xradar_dev.plot_map()
ds_gr["RHOHV"].xradar_dev.plot_map()

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

# Retrieve GR extent
extent_gr = ds_gr.xradar_dev.extent(max_distance=None, crs=None)
extent_gr

[-121.95467110437517,
 -112.12893180089884,
 28.771031109975677,
 37.06425438465268]

Preprocess SR data¶

Let’s start by cropping the SR granule on the GR extent and put data into memory:

ds_sr = ds_sr.gpm.crop(extent_gr).compute()

We can further subset the GPM overpass selecting the scans occured between the following timesteps:

start_time = datetime.datetime.strptime("2023-08-20 22:12:00", "%Y-%m-%d %H:%M:%S")
end_time = datetime.datetime.strptime("2023-08-20 22:13:45", "%Y-%m-%d %H:%M:%S")
ds_sr = ds_sr.gpm.sel(time=slice(start_time, end_time))

Now let’s display some variables illustrating the GPM overpass

ds_sr["precipRateNearSurface"].gpm.plot_map()
ds_sr["zFactorMeasured"].gpm.slice_range_at_bin(ds_sr["binClutterFreeBottom"]).gpm.plot_map(
    vmax=40,
)  # zFactorMeasuredNearSurface
ds_sr["zFactorFinalNearSurface"].gpm.plot_map(vmax=40)
ds_sr["flagPrecip"].gpm.plot_map()

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

In the following map and cross-sections, we can see that the zFactorMeasured variable is contaminated by ground clutter:

# Display the maximum observed reflectivity
ds_sr["zFactorMeasured"].max(dim="range").gpm.plot_map()
ds_sr["zFactorFinal"].max(dim="range").gpm.plot_map(vmax=40)

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

# Compute SR attenuation correction
ds_sr["zFactorCorrection"] = ds_sr["zFactorFinal"] - ds_sr["zFactorMeasured"]

# Display cross-section
ds_sr["zFactorMeasured"].isel(cross_track=24).gpm.plot_cross_section(zoom=False)
ds_sr["zFactorFinal"].isel(cross_track=24).gpm.plot_cross_section(zoom=False)
ds_sr["zFactorCorrection"].isel(cross_track=24).gpm.plot_cross_section(vmin=0, vmax=3, cmap="Spectral_r", zoom=False)

So now we mask out gates impacted by ground clutter, adn we verify with maps and cross-sections that the ground clutter is effectively removed:

# Mask SR gates with ground clutter
ds_sr["zFactorMeasured"] = ds_sr["zFactorMeasured"].gpm.mask_below_bin(
    bins=ds_sr["binClutterFreeBottom"],
    strict=False,
)  # do not mask the specified bin
ds_sr["zFactorCorrection"] = ds_sr["zFactorCorrection"].gpm.mask_below_bin(
    bins=ds_sr["binClutterFreeBottom"],
    strict=False,
)  # do not mask the specified bin

# Display the maximum observed reflectivity
ds_sr["zFactorMeasured"].max(dim="range").gpm.plot_map()
ds_sr["zFactorFinal"].max(dim="range").gpm.plot_map(vmax=40)

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

# Display cross-section with ground clutter masked out
ds_sr["zFactorMeasured"].isel(cross_track=24).gpm.plot_cross_section(zoom=False)
ds_sr["zFactorFinal"].isel(cross_track=24).gpm.plot_cross_section(zoom=False)
ds_sr["zFactorCorrection"].isel(cross_track=24).gpm.plot_cross_section(vmin=0, vmax=3, cmap="Spectral_r", zoom=False)

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 0x7f5f43953110>

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

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()

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

Retrieve SR range distance¶

The range distance from the satellite for each SR radar gate has been precomputed using information in the GPM L1B product and already included in the preprocessed GPM Zarr Store.

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

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=ds_sr["crossTrackBeamWidth"], 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 0x7f5f5c8f5ad0>

# 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 0x7f5f5c6e7410>

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=ds_sr["crossTrackBeamWidth"], 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()

# 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 0x7f5f5cc45a10>

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 0x7f5f40a7e3d0>

# 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

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

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'}

# 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']}

# 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()

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()

# 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()

# 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()

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

Here we convert reflectivities from Ku-band to X and C-band.

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()

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()

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)

# 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")

# 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()

# 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 0x7f5f5c518850>

# 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()

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

{'description': 'Flag for precipitation detection',
 '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'],
 'flag_values': [0, 1, 2, 10, 11, 12, 20, 21, 22],
 'gpm_api_decoded': 'yes',
 'gpm_api_product': '2A-DPR',
 'grid_mapping': 'crsWGS84'}

# 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 0x7f5f50302e90>

# 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 0x7f5f505dfcd0>

# 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 0x7f5f5071e990>

# 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 0x7f5f5cc16950>

# 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 0x7f5f5caead10>

# 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 0x7f5f5c460890>

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)

1305

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)

# 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))

# 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))

# 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))

# 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 0x7f5f43d3fbd0>

# 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()

# 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()

# 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()

# 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()

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',
       ...
       'gpm_cross_track_id', 'gpm_granule_id', 'gpm_id', 'lat', 'lon',
       'sweep_fixed_angle', 'sweep_number', '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)

# 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: >

# 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: >

# 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: >

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: >

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: >

# 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: >

# 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: >

# 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: >

# 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: >

# 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: >

# 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: >

# 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")

# 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")

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)

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: 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: >

# 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))

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: >

# 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: >

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.53 dB
Th ZH Calibration offset is (median): -4.58 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")

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.

References¶

Schwaller, M. R., & Morris, K. R. (2011). A Ground Validation Network for the Global Precipitation Measurement Mission. Journal of Atmospheric and Oceanic Technology, 28(3), 301–319. 10.1175/2010jtecha1403.1
Warren, R. A., Protat, A., Siems, S. T., Ramsay, H. A., Louf, V., Manton, M. J., & Kane, T. A. (2018). Calibrating Ground-Based Radars against TRMM and GPM. Journal of Atmospheric and Oceanic Technology, 35(2), 323–346. 10.1175/jtech-d-17-0128.1
Cao, Q., Hong, Y., Qi, Y., Wen, Y., Zhang, J., Gourley, J. J., & Liao, L. (2013). Empirical conversion of the vertical profile of reflectivity from Ku‐band to S‐band frequency. Journal of Geophysical Research: Atmospheres, 118(4), 1814–1825. 10.1002/jgrd.50138

Spaceborne - Ground Radar Matching and Calibration

Spaceborne-Ground Radar Calibration (Application)

Spaceborne - Ground Radar Matching and Calibration

Pyrad QPE Exercise

Spaceborne-Ground Radar Matching (Methodology)

Relevant software¶

Load SR and GR data¶

Define settings for SR/GR matching¶

Preprocess GR data with xradar¶

Preprocess SR data¶

Retrieve GR gate resolution¶

Define GR range distance interval¶

Display SR slice and GR sweep¶

Retrieve GR (lon, lat) coordinates¶

Retrieve SR (x,y,z) coordinates¶

Retrieve SR (range, azimuth, elevation) coordinates¶

Retrieve SR range distance¶

Retrieve radar gates volumes¶

Retrieve SR gate resolution¶

Retrieve Bright Band (BB) Ratio¶

Retrieve Precipitation and Hydrometeors Types¶

Convert SR Ku-band to S-band¶

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

Convert GR S-band to Ku-band¶

Identify SR gates intersecting GR sweep¶

Define SR Masks¶

Aggregate SR radar gates¶

Display aggregated SR variables¶

Retrieve the SR footprints polygons¶

Retrieve the GR gates polygons¶

Aggregate GR data on SR footprints¶

Create the SR/GR Database¶

Analyse the SR/GR Database¶

Filter the SR/GR Database¶

Determine the GR Calibration Bias¶

Next Steps¶

References¶

Citation¶