KDP estimation using Lanczos Derivatives#
import matplotlib.pyplot as plt
import numpy as np
import scipy
import xarray as xr
from matplotlib.widgets import Slider
from scipy.ndimage import convolve1d
from wradlib.util import _lanczos_differentiator as kernel
import radarpolarimetry as rp
Theoretical Background#
To estimate \(K_{DP}\) (the derivative of \(Phi_{DP}\)) we use Low-noise Lanczos differentiators.
For further reading please see Differentiation by integration using orthogonal polynomials, a survey and Low-noise Lanczos differentiators.
The Lanczos differentiators can be written in general form for any N as:
\[
f'(x^*) \approx \frac{3}{h} \sum\limits_{k=1}^{m} k \frac{f_k - f_{-k}}{m(m+1)(2m+1)}, m = \frac{N - 1}{2}
\]
For a window width of N=5 this calculates to
\[
f'(x^*) \approx \frac{3}{h} \frac{k_1(f_1-f_{-1}) + k_2(f_2-f_{-2})}{2(2+1)(4+1)}, k_1 = 1, k_2 = 2
\]
which can be written as:
\[
f'(x^*) \approx \frac{-k_2f_{-2} -k_1f_{-1} + k_1f_1 + k_2f_2}{10h}
\]
We add \(\frac{k_0}{10h}f_{0}\) for completeness (\(k_0\) is 0).
\[
f'(x^*) \approx \frac{-k_2}{10h}f_{-2} + \frac{-k_1}{10h}f_{-1} + \frac{k_0}{10h}f_{0} + \frac{k_1}{10h}f_{1} + \frac{k_2}{10h}f_{2}
\]
\(f_{-2}\), \(f_{-1}\), \(f_{1}\) and \(f_{2}\) are the discrete values of \(Phi_{DP}\) around \(f_0\) for which we want to calculate the derivation.
Assuming the factors as some discrete function \(g\), we can write the discrete convolution sum like this:
\[
(f*g)[n] = \sum\limits_{k=-m}^{m} f[n-k]g[k]
\]
This way we can use reasonable fast convolution functions to estimate the derivative of \(Phi_{DP}\).
The implementation to calculate the factors for the Lanczos differentiator can be like this:
The factors for several window-length’s:
for w in range(3, 12, 2):
print(f"Width: {w}")
print(kernel(w)[::-1])
print("------------------------------------------------------------------")
Width: 3
[-0.5 0. 0.5]
------------------------------------------------------------------
Width: 5
[-0.2 -0.1 0. 0.1 0.2]
------------------------------------------------------------------
Width: 7
[-0.10714286 -0.07142857 -0.03571429 0. 0.03571429 0.07142857
0.10714286]
------------------------------------------------------------------
Width: 9
[-0.06666667 -0.05 -0.03333333 -0.01666667 0. 0.01666667
0.03333333 0.05 0.06666667]
------------------------------------------------------------------
Width: 11
[-0.04545455 -0.03636364 -0.02727273 -0.01818182 -0.00909091 0.
0.00909091 0.01818182 0.02727273 0.03636364 0.04545455]
------------------------------------------------------------------
Synthetic truth#
dr = 0.5e3
r = np.arange(0, 100e3, dr)
azi = np.array([0])
kdp_true = np.sin(0.3e-3 * r)
kdp_true[kdp_true < 0] = 0.0
phidp_true = np.cumsum(kdp_true) * 2 * (dr * 1e-3)
# Synthetic observation of PhiDP with a random noise and gaps
phidp_raw = phidp_true + np.random.uniform(-1.5, 1.5, len(phidp_true))
gaps = np.random.uniform(0, len(r), 20).astype("int")
# phidp_raw[gaps] = np.nan
ds = xr.Dataset(
data_vars=dict(
kdp_true=(["azimuth", "range"], kdp_true[None, :].copy()),
phidp_true=(["azimuth", "range"], phidp_true[None, :].copy()),
phidp_raw=(["azimuth", "range"], phidp_raw[None, :].copy()),
),
coords=dict(
range=("range", r),
azimuth=("azimuth", azi),
),
)
Reconstruct PhiDP and KDP#
L = 11
phidp_re, kdp_re = ds.phidp_raw.copy().wrl.dp.phidp_kdp_vulpiani(winlen=L)
kdp_l = convolve1d(
np.ma.masked_invalid(phidp_raw.copy()),
kernel(L),
mode="nearest",
axis=0,
)
kdp_c = ds.phidp_raw.copy().wrl.dp.kdp_from_phidp(winlen=L)
kdp_r = ds.phidp_raw.copy().wrl.dp.kdp_from_phidp(winlen=L, method="lstsq")
phidp_c = kdp_c.wrl.dp.phidp_from_kdp()
phidp_c.attrs = phidp_re.attrs
phidp_r = kdp_r.wrl.dp.phidp_from_kdp()
phidp_r.attrs = phidp_re.attrs
Plot results#
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(211)
ds.kdp_true.plot(c="g", ls="-", label="True KDP", ax=ax)
kdp_re.plot(c="r", ls="-", label="Vulpiani KDP")
ax.plot(r, kdp_l, c="b", ls="--", label="Lanczos KDP")
kdp_c.plot(c="k", ls="-.", label="Lanczos KDP1", ax=ax)
kdp_r.plot(c="c", ls=":", label="Linregress KDP", ax=ax)
plt.grid()
lg = plt.legend()
ax = fig.add_subplot(212)
ds.phidp_true.plot(c="r", ls="--", hue="azimuth", label="True PhiDP", ax=ax)
ds.phidp_raw.plot(c="b", ls="-", hue="azimuth", label="Raw PhiDP", ax=ax)
phidp_re.plot(c="g", ls="-", hue="azimuth", label="Reconstructed PhiDP", ax=ax)
phidp_c.plot(c="c", ls="-", hue="azimuth", label="Reconstructed PhiDP1", ax=ax)
phidp_r.plot(c="y", ls="-", hue="azimuth", label="Reconstructed PhiDP2", ax=ax)
plt.grid()
lg = plt.legend(loc="lower right")
txt = plt.xlabel("Range (km)")
Interactive LANCZOS demo#
Warning
This following works only in an interactive Jupyter session! No interaction possible in the rendered documentation!
%matplotlib ipympl
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 on_change(val):
v = int(val)
dss = dsc.isel(range=slice(v - win2, v + win2 + 1))
dss = dss.assign(
lanczos=(
[
"range",
],
window[::-1],
)
)
cx = dss.phidp_raw * dss.lanczos
d0 = cx.sum(skipna=True, min_count=2) # lanczos sum
d1 = fit_line(np.arange(len(window)), dss.phidp_raw[0].values)
# upper
line3.set_data(dss.range.values[win2 : win2 + 1], np.array([d0]))
line4.set_data(dss.range.values[win2 : win2 + 1], np.array([d1]))
vline1.set_data([cx.range.min().values, cx.range.min().values], [0, 1])
vline2.set_data([cx.range.max().values, cx.range.max().values], [0, 1])
# lower
line7.set_data(dss.range.values, dss.phidp_raw.values)
line9.set_data(dss.range.values, dss.lanczos.values)
line8.set_data(cx.range.values, cx.values)
ax3.set_xlim(cx.range.min().values - 250, cx.range.max().values + 250)
global fill1, fill2
fill1.remove()
fill2.remove()
fill1 = ax4.fill_between(
cx.range.values, cx.values[0], where=cx.values[0] <= 0, color="b", alpha=0.7
)
fill2 = ax4.fill_between(
cx.range.values, cx.values[0], where=cx.values[0] >= 0, color="r", alpha=0.7
)
fig = plt.figure(figsize=(8, 8))
fig.subplots_adjust(
left=None, bottom=0.15, right=0.8, top=0.97, wspace=None, hspace=0.25
)
L = 11
v = 5
win2 = int(L // 2)
window = kernel(L)
dsc = ds.copy().assign(kdp_lanczos=kdp_c)
dss = dsc.isel(range=slice(v - win2, v + win2 + 1))
dss = dss.assign(
lanczos=(
[
"range",
],
window[::-1],
)
)
cx = dss.phidp_raw * dss.lanczos
d0 = cx.sum(skipna=True, min_count=2) # lanczos sum
d1 = fit_line(np.arange(len(window)), dss.phidp_raw)
# top plot
ax1 = fig.add_subplot(211)
ax2 = ax1.twinx()
ax1.set_ylim(-20, 70)
(line1,) = ds.phidp_raw.plot(ax=ax1, label="phidp raw")
(line2,) = ds.phidp_true.plot(ax=ax1, label="phidp true")
(line3,) = ax2.plot(
dss.range.values[win2 : win2 + 1], d0, "yo", zorder=2, label="lanczos conv sum"
)
(line4,) = ax2.plot(
dss.range.values[0 : win2 + len(d1)], d1, "kx", zorder=2, label="linregress"
)
vline1 = ax1.axvline(cx.range.min().values, 0, 1, color="k", lw=0.5)
vline2 = ax1.axvline(cx.range.max().values, 0, 1, color="k", lw=0.5)
(line5,) = ds.kdp_true.plot(ax=ax2, c="g", label="kdp true")
(line5a,) = dsc.kdp_lanczos.plot(ax=ax2, c="c", label="kdp lanczos")
# added these three lines
lns = [line1, line3, line4, line5, line5a]
labs = [ln.get_label() for ln in lns]
ax1.legend(lns, labs)
# bottom plot
ax3 = fig.add_subplot(212)
ax4 = ax3.twinx()
ax5 = rp.add_axis(ax3, pos=1.1)
x1 = np.linspace(-win2, win2, len(window))
(line7,) = dss.phidp_raw.plot(ax=ax3, marker=".", label="phidp raw")
(line8,) = cx.plot(ax=ax4, c="r", marker="o", zorder=2, label="lanczos conv")
(line9,) = dss.lanczos.plot(ax=ax5, c="b", marker="o", label="lanczos coeff")
fill1 = ax4.fill_between(
cx.range.values, cx.values[0], where=cx.values[0] <= 0, color="b", alpha=0.7
)
fill2 = ax4.fill_between(
cx.range.values, cx.values[0], where=cx.values[0] >= 0, color="r", alpha=0.7
)
# added these three lines
lns = [line7, line8, line9]
labs = [ln.get_label() for ln in lns]
ax3.legend(lns, labs, loc="upper left")
ax3.set_xlim(cx.range.min().values - 250, cx.range.max().values + 250)
ax3.set_ylim(-15, 75)
ax4.set_ylim(-4.5, 4.5)
ax1.grid(True)
ax3.grid(True)
slider_ax = plt.axes([0.1, 0.05, 0.8, 0.02])
slider = Slider(slider_ax, "Window", win2, 200 - win2 - 1, valinit=5, valfmt="%1d")
slider.on_changed(on_change)
plt.show()