Spaceborne-Ground Radar Calibration (Application)

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 and assess the GR calibration bias.

As example, we analyze a coincident GPM DPR and NEXRAD overpass over San Diego, USA, during the landfall of Tropical Storm Hilary in 2023. The NEXRAD KNKX radar data appears to be significantly miscalibrated and we will assess the calibration bias.

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.

Additionally, note that the GPM-API’s volume_matching routine can automatically download and open the coincident GPM overpass data if the SR dataset is not provided to the function.

Please read the Spaceborne-Ground Radar Matching Tutorial for a step-by-step guide through the process of obtaining spatially and temporally coincident radar samples.

Now let’s start the tutorial by importing the required packages:

import warnings

warnings.filterwarnings("ignore")

from functools import reduce

import cartopy.crs as ccrs
import fsspec
import numpy as np
import xradar as xd
import xarray as xr
from IPython.display import display
from xarray.backends.api import open_datatree
from gpm.gv import (
    calibration_summary,
    compare_maps,
    reflectivity_scatterplots,
    volume_matching,
)

np.set_printoptions(suppress=True)

Load SR and GR data¶

Now let’s open the NEXRAD GR 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={})

Now, let’s select a GR sweep to match with GPM DPR:

print("Available sweeps:", list(dt_gr.groups))

Available sweeps: ['/', '/sweep_0', '/sweep_1', '/sweep_10', '/sweep_11', '/sweep_12', '/sweep_13', '/sweep_14', '/sweep_15', '/sweep_16', '/sweep_17', '/sweep_18', '/sweep_19', '/sweep_2', '/sweep_3', '/sweep_4', '/sweep_5', '/sweep_6', '/sweep_7', '/sweep_8', '/sweep_9']

sweep_idx = 0
sweep_group = f"sweep_{sweep_idx}"  # GR sweep (elevation to be used)
ds_gr = dt_gr[sweep_group].to_dataset().compute()

and display the horizontally and vertically polarized reflectivity fields:

ds_gr["DBZH"].where(ds_gr["DBZH"] > -10).xradar_dev.plot_map()
ds_gr["DBZV"] = (ds_gr["DBZH"].gpm.idecibel / ds_gr["ZDR"].gpm.idecibel).gpm.decibel
ds_gr["DBZV"].where(ds_gr["DBZV"] > -10).xradar_dev.plot_map()
# ds_gr["ZDR"].xradar_dev.plot_map()

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

Run SR/GR volume matching¶

Here we start defining the required settings for the SR/GR volume matching procedure:

The GR radar_band controls to which frequency SR Ku-band reflectivities will be converted. Valid values are X, C, S, Ku.
The beamwidth_gr refers to the angular beam width of the GR.
The minimum reflectivity thresholds z_min_threshold_sr and z_min_threshold_gr are used to mask out SR and GR gates belows such thresholds.

radar_band = "S"
beamwidth_gr = 1
z_min_threshold_gr = 0
z_min_threshold_sr = 10

The SR/GR volume matching routine typically takes around 40-50 seconds to complete. It returns a geopandas.DataFrame with the matched aggregated reflectivities and the associated statistics.

gdf_match = volume_matching(
    ds_gr=ds_gr,
    ds_sr=ds_sr,
    z_variable_gr="DBZH",
    radar_band=radar_band,
    beamwidth_gr=beamwidth_gr,
    z_min_threshold_gr=z_min_threshold_gr,
    z_min_threshold_sr=z_min_threshold_sr,
    min_gr_range=0,
    max_gr_range=150_000,
    # gr_sensitivity_thresholds=None,
    # sr_sensitivity_thresholds=None,
    download_sr=False,  # require internet connection !
    display_quicklook=True,
)
# Note: The circle radius in the quicklook represents the min_gr_range, max_gr_range and 250_000 m range distances.

SR/GR Matching Elapsed time: 0:00:18.080791 .

display(gdf_match)

Analyse SR/GR database¶

Now let’s first define the variables names of the SR and GR reflectivities:

sr_z_column = f"SR_zFactorFinal_{radar_band}_mean"
gr_z_column = "GR_Z_mean"

Compare spatial reflectivity fields¶

Here we start comparing the spatial reflectivity fields without applying restrictive filtering criteria:

fig = compare_maps(
    gdf_match,
    sr_column=sr_z_column,
    gr_column=gr_z_column,
    sr_label="SR Reflectivity (dBz)",
    gr_label="GR Reflectivity (dBz)",
    cmap="Spectral_r",
    unified_color_scale=True,
    vmin=15,
    # vmax=40
)
fig.tight_layout()

If you wish to create a cartopy map, specify the 'projection' in the subplot_kwargs argument:

ccrs_gr_aeqd = ccrs.AzimuthalEquidistant(
    central_longitude=ds_gr["longitude"].item(),
    central_latitude=ds_gr["latitude"].item(),
)
subplot_kwargs = {}
subplot_kwargs["projection"] = ccrs_gr_aeqd
fig = compare_maps(
    gdf_match,
    sr_column=sr_z_column,
    gr_column=gr_z_column,
    sr_label="SR Reflectivity (dBz)",
    gr_label="GR Reflectivity (dBz)",
    cmap="Spectral_r",
    unified_color_scale=True,
    vmin=15,
    # vmax=40
    subplot_kwargs=subplot_kwargs,
)
fig.tight_layout()

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

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

Generate calibration summary¶

We now create a figure comparing SR/GR aggregated reflectivities volume-by-volume and displaying the overall distributions.

gr_z_column = "GR_Z_mean"
sr_z_column = f"SR_zFactorFinal_{radar_band}_mean"

fig = calibration_summary(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    # Histogram options
    bin_width=2,
    # Scatterplot options
    hue_column="SR_fraction_clutter",
    # gr_range=[15, 50]
    # sr_range=[15, 50]
    marker="+",
    cmap="Spectral",
)
fig.tight_layout()

Investigate filtering criteria¶

When comparing SR and GR data or trying to determine an accurate GR calibration bias, it’s necessary to define a set of filtering criteria. In the figures below, we perform exploratory data analysis to investigate the relationships between the SR/GR reflectivity deviations and sets of variables characterizing radar measurements and SR/GR volume properties. The patterns and deviations observed in the scatterplots will be used to define a set of filtering critera in the next section of the tutorial.

# Matching characteristics
hue_columns = [
    "VolumeRatio",
    "VolumeDiff",
    "GR_vres_max",
    "GR_range_max",
    "GR_gate_volume_sum",
    "SR_gate_volume_sum",
    "SR_vres_sum",
    "GR_fraction_covered_area",
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

# Value validity & sensitivity
hue_columns = [
    "GR_counts",
    "GR_counts_valid",
    "SR_counts",
    "SR_counts_valid",
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

hue_columns = [
    "SR_zFactorFinal_Ku_fraction_above_12dBZ",
    "SR_zFactorFinal_Ku_fraction_above_14dBZ",
    "GR_Z_fraction_above_10dBZ",
    "GR_Z_fraction_above_12dBZ",
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

# NUBF & Variability
hue_columns = [
    "SR_zFactorFinal_Ku_range",
    "GR_Z_range",
    "SR_zFactorFinal_Ku_std",
    "GR_Z_std",
    "SR_zFactorFinal_Ku_cov",
    "GR_Z_cov",
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

# Attenuation
hue_columns = [
    "SR_pathAtten",
    "SR_piaFinal",
    "SR_zFactorCorrection_Ku_max",
    "SR_zFactorCorrection_Ku_mean",
    "SR_reliabFlag",
    "GR_Z_cumsum_max",
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

# Hydrometeors
hue_columns = [
    "SR_fraction_no_precip",
    "SR_fraction_clutter",
    "SR_fraction_rain",
    "SR_fraction_snow",
    "SR_fraction_hail",
    "SR_fraction_below_isotherm",
    "SR_fraction_above_isotherm",
    "SR_fraction_melting_layer",
    # 'SR_airTemperature_min',
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

# SR Quality Flags
hue_columns = [
    "SR_flagPrecipitationType",
    "SR_qualityTypePrecip",
    "SR_qualityBB",
    "SR_flagPrecip",
    "SR_qualityFlag",
    "SR_dataQuality",
]
fig = reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

Define SR/GR filtering criteria¶

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 masks
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_rain"] == 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"] < 15,
    # 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 > 5 minutes
    # - time_difference is in seconds !
    gdf_match["time_difference"] < 60 * 5,
    # 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
reflectivity_scatterplots(
    df=gdf_match,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns="filtering_mask",
    marker="o",
    s=1,
    cmap="viridis_r",
)

Analyze filtered SR/GR database¶

We now filter the SR/GR Database, compare the reflectivity fields and display the new calibration summary:

gdf_final = gdf_match[mask_final]

fig = compare_maps(
    gdf_final,
    sr_column=sr_z_column,
    gr_column=gr_z_column,
    sr_label="SR Reflectivity (dBz)",
    gr_label="GR Reflectivity (dBz)",
    cmap="Spectral_r",
    unified_color_scale=True,
    vmin=15,
    # vmax=40
)
fig.tight_layout()

fig = calibration_summary(
    df=gdf_final,
    sr_z_column=sr_z_column,
    gr_z_column=gr_z_column,
    # Histogram options
    bin_width=2,
    # Scatterplot options
    # hue_column="GR_gate_volume_sum",
    hue_column="GR_range_mean",
    # gr_range=[15, 50]
    # sr_range=[15, 50]
    marker="+",
    cmap="Spectral",
)
fig.tight_layout()

Exploratory data analysis can again be performed to check if the filtering criteria were effective to discard the undesired SR/GR matched volumes and if the filtering criteria should be further expanded.

# Attenuation
hue_columns = [
    "SR_pathAtten",
    "SR_piaFinal",
    "SR_reliabFlag",
    "SR_zFactorCorrection_Ku_max",
    "SR_zFactorCorrection_Ku_mean",
    "GR_Z_cumsum_max",
]
fig = reflectivity_scatterplots(
    df=gdf_final,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    hue_columns=hue_columns,
    ncols=2,
)
fig.tight_layout()

Determine GR calibration bias¶

The GR calibration bias can be obtained by averaging the difference between the matched SR/GR reflectivity measurements:

# Compute average offset
z_offset = np.nanmean(gdf_final[sr_z_column] - gdf_final[gr_z_column]).round(2)
z_offset_robust = np.nanmedian(gdf_final[sr_z_column] - gdf_final[gr_z_column]).round(2)

print(f"Th ZH Calibration offset is (mean): {z_offset} dBZ")
print(f"Th ZH Calibration offset is (median): {z_offset_robust} dBZ")

Th ZH Calibration offset is (mean): 3.94 dBZ
Th ZH Calibration offset is (median): 3.78 dBZ

Volume matching with calibrated GR¶

Now, we remove the estimated GR calibration bias, we rerun the SR/GR volume matching procedure, we compare the matched reflectivities and we assess if we have effectively removed the bias by generating the calibration summary. Please note that here below we don’t filter out undesired SR/GR volumes, but this is operation is recommended once you defined appropriate robust filtering criteria !

# Apply offset
ds_gr["DBZH_c"] = ds_gr["DBZH"] + z_offset

# Apply matching
gdf_match_corrected = volume_matching(
    ds_gr=ds_gr,
    ds_sr=ds_sr,
    z_variable_gr="DBZH_c",
    radar_band=radar_band,
    beamwidth_gr=beamwidth_gr,
    z_min_threshold_gr=z_min_threshold_gr,
    z_min_threshold_sr=z_min_threshold_sr,
    min_gr_range=0,
    max_gr_range=150_000,
    # gr_sensitivity_thresholds=None,
    # sr_sensitivity_thresholds=None,
    download_sr=False,  # require internet connection !
)

SR/GR Matching Elapsed time: 0:00:18.247898 .

# Compare reflectivities
fig = compare_maps(
    gdf_match_corrected,
    sr_column=sr_z_column,
    gr_column=gr_z_column,
    sr_label="SR Reflectivity (dBz)",
    gr_label="GR Reflectivity (dBz)",
    cmap="Spectral_r",
    unified_color_scale=True,
    vmin=15,
    # vmax=40
)
fig.tight_layout()

# Display difference in reflectivity betweeen SR and GR
gdf_match_corrected["z_offset"] = gdf_match_corrected[sr_z_column] - gdf_match_corrected[gr_z_column]
gdf_match_corrected.plot(column="z_offset", legend=True, cmap="RdBu", vmin=-1.5, vmax=1.5)

<Axes: >

# Display calibration summary
fig = calibration_summary(
    df=gdf_match_corrected,
    gr_z_column=gr_z_column,
    sr_z_column=sr_z_column,
    # Histogram options
    bin_width=2,
    # Scatterplot options
    hue_column="SR_fraction_clutter",
    marker="+",
    cmap="Spectral",
)
fig.tight_layout()

Next steps¶

We encourage you to explore and adapt the code provided in this tutorial to infer calibration bias also for the vertically polarized channel, iterate over GR sweeps, analyze multiple GPM overpasses, and assess the long-term calibration bias of your GR network.

Please share any insights or suggestions for improving the matching procedure or filtering criteria with the GPM-API community so we can collaboratively enhance the routines.

We hope you enjoyed the 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¶

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

ARCO Radar Datasets

QPE & QVPs

Spaceborne - Ground Radar Matching and Calibration

Spaceborne-Ground Radar Matching (Methodology)