Match spaceborn SR (GPM/TRMM) with ground radars GR#

The idea is to match ground radar (GR) and space-born radar (SR) measurements in order to create spatially and temporally coicident samples without interpolation. The procedure had been suggested by Schwaller and Morris (2011) and is based on the adaption by Warren, et. al. (2017).

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

Quote Warren, et.al: “[…] In this approach, intersections between indi vidual 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. […]”.

This becomes clearer in Fig. 3:

Schwaller, MR, and Morris, KR. 2011. A ground validation network for the Global Precipitation Measurement mission. J. Atmos. Oceanic Technol., 28, 301-319.

Warren, R.A., A. Protat, S.T. Siems, H.A. Ramsay, V. Louf, M.J. Manton, and T.A. Kane, 0: Calibrating ground-based radars against TRMM and GPM. J. Atmos. Oceanic Technol., 0,

Conventions#

This code is based on the following conventions:

gr indicates ground radar
sr indicates space-born precipitation radar (TRMM or GPM)

The base routines are designed to process one GR sweep at a time. If a full GR volume with nelev of sweeps is available, you can iterate over each sweep. In this code, ee is an index that points to one of the nelev sweeps/elevation angles. Accordingly, a GR data set will be organised as an array of shape (nelev_gr, nray_gr, ngate_gr).

A SR data set is typically organised as arrays with dimensions (nscan_sr, nray_sr, ngate_sr).

[1]:

import warnings

#warnings.filterwarnings("ignore")
import wradlib as wrl
import matplotlib.pyplot as plt
import matplotlib as mpl

try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()

import numpy as np
import datetime as dt
from osgeo import osr
import xarray as xr
import xradar as xd

Acquire datafiles#

[2]:

# define GPM data set
gpm_file = wrl.util.get_wradlib_data_file(
    "gpm/2A-CS-151E24S154E30S.GPM.Ku.V7-20170308.20141206-S095002-E095137.004383.V05A.HDF5"
)
# define matching ground radar file
gr2gpm_file = wrl.util.get_wradlib_data_file("hdf5/IDR66_20141206_094829.vol.h5")

# define TRMM data sets
trmm_2a23_file = wrl.util.get_wradlib_data_file(
    "trmm/2A-CS-151E24S154E30S.TRMM.PR.2A23.20100206-S111425-E111526.069662.7.HDF"
)
trmm_2a25_file = wrl.util.get_wradlib_data_file(
    "trmm/2A-CS-151E24S154E30S.TRMM.PR.2A25.20100206-S111425-E111526.069662.7.HDF"
)

# define matching ground radar file
gr2trmm_file = wrl.util.get_wradlib_data_file("hdf5/IDR66_20100206_111233.vol.h5")

Downloading file 'gpm/2A-CS-151E24S154E30S.GPM.Ku.V7-20170308.20141206-S095002-E095137.004383.V05A.HDF5' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/gpm/2A-CS-151E24S154E30S.GPM.Ku.V7-20170308.20141206-S095002-E095137.004383.V05A.HDF5' to '/home/runner/work/wradlib/wradlib/wradlib-data'.
Downloading file 'hdf5/IDR66_20141206_094829.vol.h5' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/IDR66_20141206_094829.vol.h5' to '/home/runner/work/wradlib/wradlib/wradlib-data'.
Downloading file 'trmm/2A-CS-151E24S154E30S.TRMM.PR.2A23.20100206-S111425-E111526.069662.7.HDF' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/trmm/2A-CS-151E24S154E30S.TRMM.PR.2A23.20100206-S111425-E111526.069662.7.HDF' to '/home/runner/work/wradlib/wradlib/wradlib-data'.
Downloading file 'trmm/2A-CS-151E24S154E30S.TRMM.PR.2A25.20100206-S111425-E111526.069662.7.HDF' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/trmm/2A-CS-151E24S154E30S.TRMM.PR.2A25.20100206-S111425-E111526.069662.7.HDF' to '/home/runner/work/wradlib/wradlib/wradlib-data'.
Downloading file 'hdf5/IDR66_20100206_111233.vol.h5' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/hdf5/IDR66_20100206_111233.vol.h5' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

Set SR and GR parameters#

[3]:

# Space-born precipitation radar parameters
sr_pars = {
    "trmm": {
        "zt": 402500.0,  # orbital height of TRMM (post boost)   APPROXIMATION!
        "dr": 250.0,  # gate spacing of TRMM
        "gr_file": gr2trmm_file,
    },
    "gpm": {
        "zt": 407000.0,  # orbital height of GPM                 APPROXIMATION!
        "dr": 125.0,  # gate spacing of GPM
        "gr_file": gr2gpm_file,
    },
}

[4]:

# Set parameters for this procedure
bw_sr = 0.71  # SR beam width
platf = "gpm"  # SR platform/product: one out of ["gpm", "trmm"]
zt = sr_pars[platf]["zt"]  # SR orbit height (meters)
dr_sr = sr_pars[platf]["dr"]  # SR gate length (meters)
gr_file = sr_pars[platf]["gr_file"]
ee = 2  # Index that points to the GR elevation angle to be used

Data Input#

Ground Radar Reflectivity matched#

[63]:

mean = obj3.mean(wrl.trafo.idecibel(swp.DBZH).values.ravel())
obj3.zdata.trg.geo["mean"] = wrl.trafo.decibel(mean)
ref_gr_m = obj3.zdata.trg.geo["mean"]

[64]:

mean = obj3.mean(wrl.trafo.idecibel(ref_gr2sr).values.ravel())
obj3.zdata.trg.geo["mean"] = wrl.trafo.decibel(mean)
ref_gr2sr_m = obj3.zdata.trg.geo["mean"]

[65]:

# todo: title fig.suptitle("GR vs SR (to_gr)") etc

[66]:

plot_matched_refl(ref_sr2gr_m, ref_gr_m, xyp_sr)
plt.suptitle("GR vs SR (to_gr)")

[66]:

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

../../_images/notebooks_workflow_recipe3_103_1.png

[67]:

plot_matched_refl(ref_sr_m, ref_gr2sr_m, xyp_sr)
plt.suptitle("GR (to_sr) vs SR")

[67]:

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

../../_images/notebooks_workflow_recipe3_104_1.png

Diagnostic Plots#

[68]:

volcol = sr_match2_sum.vol_sr

[69]:

fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(121, aspect="equal")
plt.scatter(
    ref_gr_m.values, ref_sr2gr_m.values.ravel(), marker="+", c=volcol, cmap="turbo"
)
plt.colorbar(label="Matching Volume [m³]")
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: {np.nanmedian(ref_gr_m.values - ref_sr2gr_m.values.ravel())}")
ax = fig.add_subplot(122)
plt.hist(
    ref_gr_m[ref_gr_m > 10], bins=np.arange(-10, 50, 5), edgecolor="None", label="GR"
)
plt.hist(
    ref_sr2gr_m[ref_sr2gr_m > -10],
    bins=np.arange(-10, 50, 5),
    edgecolor="red",
    facecolor="None",
    label="SR",
)
plt.xlabel("Reflectivity (dBZ)")
plt.legend()
fig.suptitle("GR vs SR (to_gr)")

[69]:

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

../../_images/notebooks_workflow_recipe3_107_1.png

[70]:

fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(121, aspect="equal")
plt.scatter(ref_gr2sr_m, ref_sr_m, marker="+", c=volcol, cmap="turbo")
plt.colorbar(label="Matching Volume [m³]")
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: {np.nanmedian(ref_gr2sr_m.values - ref_sr_m.values.ravel())}")
ax = fig.add_subplot(122)
plt.hist(
    ref_gr2sr_m[ref_gr2sr_m > 10],
    bins=np.arange(-10, 50, 5),
    edgecolor="None",
    label="GR",
)
plt.hist(
    ref_sr_m[ref_sr_m > -10],
    bins=np.arange(-10, 50, 5),
    edgecolor="red",
    facecolor="None",
    label="SR",
)
plt.xlabel("Reflectivity (dBZ)")
plt.legend()
fig.suptitle("GR (to_sr) vs SR")

[70]:

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

../../_images/notebooks_workflow_recipe3_108_1.png

Match spaceborn SR (GPM/TRMM) with ground radars GR

Contents

Match spaceborn SR (GPM/TRMM) with ground radars GR#

Conventions#

Acquire datafiles#

Set SR and GR parameters#

Data Input#

Ground Radar GR#

Read GR data into datatree using xradar#

Extract relevant GR data and meta-data#

Georeference GR (AEQD)#

Assign lonlat coords to GR#

Get GR Grid Polygons#

Calculate GR Bounding Box#

Satellite data SR#

Read GPM data#

Georeference SR data#

Set fundamental georeferencing parameters#

Subset relevant SR data#

Subset relevant SR data with precip flag#

SR Parallax Correction#

Compute spherical coordinates of SR bins with regard to GR#

Compute SR and GR pulse volumes#

Calculate distance from orbit \(r_s\)#

SR pulse volume#

GR pulse volume#

Calculate horizontal and vertical dimensions \(R_s\) and \(D_s\) of SR bins#

BrightBand Handling#

Median Brightband Width/Height#

Convert SR Ku reflectivities to S-band#

Convert S-band GR reflectivities to Ku-band#

Matching SR/GR#

Identify which SR rays actually intersect with the GR sweep#

New version matching#

Identify which GR bins actually intersect with the SR footprint#

SR Reflectivity matched#

Ground Radar Reflectivity matched#

Diagnostic Plots#