KDP Derivation#
Specific differential Phase KDP.
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import radarpolarimetry as rp
Load Dataset#
stage = "03_phidp"
dtime = "20170810T1100"
fname = rp.data_path() / "lecture" / f"{stage}_{dtime}.zarr"
swp = xr.open_dataset(fname)
display(swp)
<xarray.Dataset> Size: 64MB
Dimensions: (azimuth: 360, range: 750)
Coordinates: (12/16)
* azimuth (azimuth) float64 3kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
elevation (azimuth) float64 3kB ...
time (azimuth) datetime64[ns] 3kB ...
* range (range) float32 3kB 100.0 300.0 ... 1.497e+05 1.499e+05
bins (azimuth, range) float32 1MB ...
gr (azimuth, range) float64 2MB ...
... ...
altitude float64 8B ...
crs_wkt int64 8B ...
latitude float64 8B ...
longitude float64 8B ...
NOISE_LEVEL float64 8B ...
sweep_mode <U20 80B ...
Data variables: (12/34)
CMAP (azimuth, range) float64 2MB ...
CMAP2 (azimuth, range) uint8 270kB ...
DBTH (azimuth, range) float64 2MB ...
DBTV (azimuth, range) float64 2MB ...
DBZH (azimuth, range) float64 2MB ...
DBZV (azimuth, range) float64 2MB ...
... ...
SYSPHI4 (azimuth) float64 3kB ...
VRADH (azimuth, range) float64 2MB ...
VRADV (azimuth, range) float64 2MB ...
WRADH (azimuth, range) float64 2MB ...
WRADV (azimuth, range) float64 2MB ...
ZDR (azimuth, range) float64 2MB ...
Attributes:
source: gamic
ant_gain_h: 44.0
ant_gain_v: 44.0
noise_power_h: -1.39495
noise_power_v: -1.09706
rx_loss_h: 0.73
rx_loss_v: 0.73L = 11
phi = "PHIDP_CORR1"
KDP from Linear Regression#
Slow, terribly slow! Not usable for fast processing. Do not use this for whole PPI!
import scipy
def fit_line(x, y):
"""slope calculation
Parameters: x and y array"""
slope, intercept, r_value, p_value, err = scipy.stats.linregress(x, y)
return slope
def kdp_calc(p, wl=5, dr=1.0):
"""K_DP calculation for all time steps
Parameter: p: radar matrix for all time steps"""
# print(p.shape)
# get x-array over last dimension of source-array (here [2])
x = np.arange(0, p.shape[1])
# kdp = np.ones_like(p)*(-10)
kdp = np.ones_like(p) * (0)
# iterate over azimuths
for j in range(np.shape(p)[0]):
# print(j)
# iterate along ray
for i in range(np.shape(p)[1] - wl):
# fit
slope = fit_line(x[i : i + wl - 1], p[j, i : i + wl - 1])
# Durch 2 wg Def von Kdp, / dr (Einheit deg/km)
kdp[j, i + int((wl - 1) / 2)] = slope / 2.0 / dr
return kdp
We use only one radial to show the slowness.
%%time
dr = swp.range.diff("range").median() / 1000.0
kdp_linregress = kdp_calc(swp[phi][65].values[None, :], wl=L, dr=dr.values)
CPU times: user 210 ms, sys: 1 μs, total: 210 ms
Wall time: 208 ms
KDP using Least Squares#
swp[phi].wrl.vis.plot(levels=np.arange(0, 30, 1))
<matplotlib.collections.QuadMesh at 0x7f31d6813b60>
dr = swp[phi].range.diff("range").median("range").values / 1000.0
swp[phi]
<xarray.DataArray 'PHIDP_CORR1' (azimuth: 360, range: 750)> Size: 2MB
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]], shape=(360, 750))
Coordinates: (12/16)
* azimuth (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
elevation (azimuth) float64 3kB ...
time (azimuth) datetime64[ns] 3kB ...
* range (range) float32 3kB 100.0 300.0 500.0 ... 1.497e+05 1.499e+05
bins (azimuth, range) float32 1MB ...
gr (azimuth, range) float64 2MB ...
... ...
altitude float64 8B 99.5
crs_wkt int64 8B 0
latitude float64 8B 50.73
longitude float64 8B 7.072
NOISE_LEVEL float64 8B -24.7
sweep_mode <U20 80B 'azimuth_surveillance'
Attributes:
format: UV16
is_dft: 0
unit: deg
units: degrees
long_name: Differential phase HV
standard_name: radar_differential_phase_hv
_Undetect: 0.0%%time
kdp_lstsq = swp[phi].wrl.dp.kdp_from_phidp(winlen=L, min_periods=2, method="lstsq")
CPU times: user 5.38 s, sys: 273 ms, total: 5.65 s
Wall time: 1.36 s
KDP using Lanczos-Derivatives#
By using convolution kernel the Lanczos-Derivative can be calculated:
The processing is some orders of magnitude faster.
%%time
kdp_lanczos = swp[phi].wrl.dp.kdp_from_phidp(winlen=L, min_periods=2)
CPU times: user 70.9 ms, sys: 7.94 ms, total: 78.8 ms
Wall time: 78.3 ms
KDP Vulpiani-Method#
Iterative approach using first guess derivative.
%%time
phidp_v, kdp_vulpiani = swp[phi].wrl.dp.phidp_kdp_vulpiani(winlen=L)
CPU times: user 72.2 ms, sys: 4.01 ms, total: 76.3 ms
Wall time: 75.3 ms
Comparison#
domain wise for KDP
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 10))
kdp_lstsq.wrl.vis.plot(ax=axes.flat[0], vmin=-2, vmax=3)
axes.flat[0].set_title("KDP least squares")
kdp_vulpiani.where(~np.isnan(swp[phi])).wrl.vis.plot(ax=axes.flat[1], vmin=-2, vmax=3)
axes.flat[1].set_title("KDP vulpiani")
kdp_lanczos.wrl.vis.plot(ax=axes.flat[2], vmin=-2, vmax=3)
axes.flat[2].set_title("KDP lanczos")
swp.KDP.wrl.vis.plot(ax=axes.flat[3], vmin=-2, vmax=3)
axes.flat[3].set_title("KDP sigproc")
fig.tight_layout()
Exercise 1: Check other window sizes and corrected phases#
Solution 1: Check other window sizes and corrected phases#
L = 21
phi = "PHIDP_CORR2"
kdp_lstsq = swp[phi].wrl.dp.kdp_from_phidp(winlen=L, min_periods=2, method="lstsq")
kdp_lanczos = swp[phi].wrl.dp.kdp_from_phidp(winlen=L, min_periods=2)
phidp_v, kdp_vulpiani = swp[phi].wrl.dp.phidp_kdp_vulpiani(winlen=L)
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 10))
kdp_lstsq.wrl.vis.plot(ax=axes.flat[0], vmin=-2, vmax=3)
axes.flat[0].set_title("KDP least squares")
kdp_vulpiani.where(~np.isnan(swp[phi])).wrl.vis.plot(ax=axes.flat[1], vmin=-2, vmax=3)
axes.flat[1].set_title("KDP vulpiani")
kdp_lanczos.wrl.vis.plot(ax=axes.flat[2], vmin=-2, vmax=3)
axes.flat[2].set_title("KDP lanczos")
swp.KDP.wrl.vis.plot(ax=axes.flat[3], vmin=-2, vmax=3)
axes.flat[3].set_title("KDP sigproc")
fig.tight_layout()
Add and save data#
swp = swp.assign({"KDP_LANC": kdp_lanczos})
swp = swp.assign({"KDP_VULP": kdp_vulpiani})
dim0 = "azimuth"
dtime = swp.time.min(dim0).dt.floor("5min").dt.strftime("%Y%m%dT%H%M").values
stage = "04_kdp"
fname = rp.data_path() / "lecture" / f"{stage}_{dtime}.zarr"
print(fname)
swp.to_zarr(fname, mode="w", zarr_version=2)
/builds/_NRRzFk9l/0/grp_meteo/radar/radarpolarimetry/data/lecture/04_kdp_20170810T1100.zarr
<xarray.backends.zarr.ZarrStore at 0x7f31e3caf100>