ZPHI-Method#

Testud, J., Le Bouar, E., Obligis, E., & Ali-Mehenni, M. (2000). The Rain Profiling Algorithm Applied to Polarimetric Weather Radar, Journal of Atmospheric and Oceanic Technology, 17(3), 332-356. Retrieved Nov 24, 2021, from https://doi.org/10.1175/1520-0426(2000)017<0332:TRPAAT>2.0.CO;2

Ryzhkov, A., M. Diederich, P. Zhang, and C. Simmer, 2014: Potential Utilization of Specific Attenuation for Rainfall Estimation, Mitigation of Partial Beam Blockage, and Radar Networking. J. Atmos. Oceanic Technol., 31, 599–619, https://doi.org/10.1175/JTECH-D-13-00038.1.

Diederich, M., Ryzhkov, A., Simmer, C., Zhang, P., & Trömel, S. (2015). Use of Specific Attenuation for Rainfall Measurement at X-Band Radar Wavelengths.: Part I: Radar Calibration and Partial Beam Blockage Estimation. Journal of Hydrometeorology, 16(2), 487–502. https://www.jstor.org/stable/24914953

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import radarpolarimetry as rp

Plot Overview#

fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(121)
im1 = rp.plot_moment(swp.PHIDP_CORR1, ax=ax1)
t = plt.title(r"Corrected $\phi_{DP}$")
t.set_y(1.1)

ax2 = fig.add_subplot(122)
im2 = rp.plot_moment(swp.DBZH, ax=ax2)
t = plt.title(r"Corrected $Z_{H}$")
t.set_y(1.1)
fig.suptitle(swp.time.min().values.astype("<M8[s]"), fontsize=14)
fig.tight_layout()

../../_images/3158f763a886908dbb4b28453f6f77dfdfdf727e3b91a815fd68737e2c5d99c0.png

Formulas ZPHI-Method#

see Testud et.al. (chapter 4. p. 339ff.), Ryzhkov et.al. (chapter 3. p. 601 ff) and Diederich et.al. (chapter 3. p. 492 ff).

There is an equational difference between the first and the other two papers, which can be solved like this:

\(\begin{equation} f\Delta\phi_{DP} = 10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}} - 1 \tag{1} \end{equation}\)

\(\begin{equation} C(b, PIA) = \exp[{0.23 \cdot b \cdot (PIA)}] - 1 \tag{2} \end{equation}\)

with

\(\begin{equation} PIA = \alpha \cdot \Delta\phi_{DP} \tag{3} \end{equation}\)

\(\begin{equation} C(b, PIA) = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{4} \end{equation}\)

Both expressions are used equivalently:

\(\begin{equation} 10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}} - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{5} \end{equation}\)

Using logarithmic identities:

\(\begin{equation} \ln {u^r} = r \cdot \ln {u} \tag{6a} \end{equation}\)

\(\begin{equation} \exp {\ln x} = x \tag{6b} \end{equation}\)

the left hand side can be further expressed as:

\(\begin{equation} \exp [\ln {10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}}] - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{7a} \end{equation}\)

\(\begin{equation} \exp[0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP} \cdot \ln {10}] - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{7b} \end{equation}\)

leading to equality

\(\begin{equation} \exp[0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}] - 1 = \exp[{0.23 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}}] - 1 \tag{7c} \end{equation}\)

We are using the following equations derived from the above papers:

\(\begin{equation} \Delta \phi_{DP} = \phi_{DP}(r2) - \phi_{DP} (r1) \tag{1} \end{equation}\)

\(\begin{equation} f\Delta\phi_{DP} = 10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}} - 1 \tag{2} \end{equation}\)

\(\begin{equation} A_{H}(r) = \frac{\left [Z_a(r) \right ]^b \cdot f(\Delta \phi_{DP})}{0.46b \int_{r1}^{r2} \left [Z_a(s) \right ]^b ds + f(\Delta \phi_{DP}) \cdot 0.46b \int_{r}^{r2} \left [Z_a(s) \right ]^b ds} \tag{3} \end{equation}\)

\(\begin{equation} \phi_{DP}^{cal}(r_i, \alpha) = 2 \cdot \int_{r1}^{r2} \frac{A_H(s; \alpha)}{\alpha}ds \tag{4} \end{equation}\)

\(\begin{equation} \Delta = \sum_{i=1}^{N} \left|{\phi_{DP}^{cal}(r_i, \alpha, b)- \phi_{DP}^{m}(r_i)}\right| \tag{5} \end{equation}\)

\(\begin{equation} A_{H} = \beta Z_{H}^{b} \tag{6} \end{equation}\)

\(\begin{equation} A_{H} = \alpha K_{DP} \tag{7} \end{equation}\)

Retrieving \(\Delta \phi_{DP}\)#

See Delta PhiDP Notebook.

Todo: Use cphase from delta phi notebook, so we do not recalculate here.

CMAP = swp.CMAP.astype("uint16")
noise = ((CMAP >> 0) & 0x1).astype(bool)
clutter = ((CMAP >> 2) & 0x1).astype(bool)
masked = noise | clutter
valid = ~masked
phi = swp.PHIDP.where((swp.CMAP2 == 0) & valid & (swp.RHOHV > 0.9))
rng = 2000.0

cphase = rp.phase_zphi(phi, rng=rng)
dphi = cphase.dphi.where(cphase.dphi >= 0, 1e-4)

\(\Delta \phi_{DP}\) - Linear Plots#

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8))
cphase.first.plot(label="first", ax=ax1)
cphase.last.plot(label="last", ax=ax1)
ax1.grid()
ax1.legend()
dphi.plot(ls="-", marker=".", label="delta", ax=ax2)
ax2.grid()
ax2.legend()
plt.tight_layout()

../../_images/26afa2180dec93a938a0f2cc3617d0521928410a2bf0e3a93ecdb23ba59ad240.png

\(\Delta \phi_{DP}\) - Polar Plots#

This visualizes first and last indizes including \(\Delta \phi_{DP}\).

fig = plt.figure(figsize=(20, 14))
ax1 = plt.subplot(231, projection="polar")
ax2 = plt.subplot(232, projection="polar")
ax3 = plt.subplot(233, projection="polar")
# set the lable go clockwise and start from the top
ax1.set_theta_zero_location("N")
ax2.set_theta_zero_location("N")
ax3.set_theta_zero_location("N")
# clockwise
ax1.set_theta_direction(-1)
ax2.set_theta_direction(-1)
ax3.set_theta_direction(-1)

theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
ax1.plot(theta, cphase.first_idx, color="b", linewidth=3)
_ = ax1.set_title("Delta PHIDP - First Index")
ax2.plot(theta, cphase.last_idx, color="r", linewidth=3)
_ = ax2.set_title("Delta PHIDP - Last Index")
ax3.plot(theta, dphi, color="g", linewidth=3)
_ = ax3.set_title("Delta PHIDP")

../../_images/9dfe7e157702db968ec3cf0d5180e09b412304b237ac877c71a7ed3c7d6ddc08.png

Calculating \(f\Delta\phi_{DP}\)#

\[f\Delta\phi_{DP} = 10^{0.1 \cdot b \cdot \alpha \cdot \Delta\phi_{DP}} - 1\]

alphax = np.arange(0.1, 0.35, 0.01)
alphax = xr.DataArray(alphax, coords={"alpha": alphax}, dims="alpha")
betax = 0.05
bx = 0.78
# need to expand alphax to dphi-shape
fdphi = 10 ** (0.1 * bx * alphax * dphi) - 1
display(fdphi)

<xarray.DataArray (alpha: 25, azimuth: 360)> Size: 72kB
dask.array<sub, shape=(25, 360), dtype=float64, chunksize=(25, 180), chunktype=numpy.ndarray>
Coordinates:
  * alpha        (alpha) float64 200B 0.1 0.11 0.12 0.13 ... 0.31 0.32 0.33 0.34
  * azimuth      (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation    (azimuth) float64 3kB 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
    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

Calculating Reflectivity Integrals/Sums#

\[za(r) = \left[Z_a(r) \right ]^b\]

\[iza(r,r2) = 0.46 \cdot b \cdot \int_{r}^{r2} \left [Z_a(s) \right ]^b ds\]

We do not restrict (mask) the reflectivities for now, but switch between DBTH and DBZH to see the difference.

zhraw = swp.DBTH.where((swp.CMAP2 == 0) & valid & (swp.RHOHV > 0.9))
zhraw = zhraw.where((swp.range > cphase.start_range) & (swp.range < cphase.stop_range))
zhraw.plot(x="azimuth", vmin=-10, vmax=100)

<matplotlib.collections.QuadMesh at 0x7f3a921f1810>

../../_images/7b7a008685b8eef8a69591b4dd968d53a8f4f5b3dade7154a3449f202420a566.png

# calculate linear reflectivity and ^b
zax = zhraw.wrl.trafo.idecibel().fillna(1e-4)
zh_cal = zax.wrl.trafo.decibel()
zh_cal.wrl.vis.plot()

<matplotlib.collections.QuadMesh at 0x7f3a90172d50>

../../_images/5d3cb7f02a15abc1a0a414bfedc2036c18dfb69042496568f71b9adc0d5ebb54.png

za = zax**bx
za.wrl.vis.plot(vmax=1000)

<matplotlib.collections.QuadMesh at 0x7f3a6876a0d0>

../../_images/45cd4ea7cfa5c55a5b4efeef4c8c864828663585012a8f00a225cf070b4798d8.png

# set masked to zero for integration
za_zero = za.fillna(0)
za_zero.wrl.vis.plot(vmax=1000)

<matplotlib.collections.QuadMesh at 0x7f3a68628e10>

# integrate
iza_x = 0.46 * bx * za_zero.pipe(rp.x_cumtrapz)
iza = iza_x.max("range") - iza_x
iza.wrl.vis.plot()
display(iza)

range res [km]: 0.2

<xarray.DataArray 'DBTH' (azimuth: 360, range: 750)> Size: 2MB
dask.array<sub, shape=(360, 750), dtype=float64, chunksize=(180, 750), chunktype=numpy.ndarray>
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 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
  * range        (range) float32 3kB 100.0 300.0 500.0 ... 1.497e+05 1.499e+05
    bins         (azimuth, range) float32 1MB dask.array<chunksize=(180, 375), meta=np.ndarray>
    gr           (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    ...           ...
    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:         UV8
    is_dft:         0
    unit:           dBZ
    units:          dBZ
    long_name:      Total power H (uncorrected reflectivity)
    standard_name:  radar_equivalent_reflectivity_factor_h
    _Undetect:      0.0

../../_images/adaf0be3ec1c2a81789311f22031c24404703bedd5da88cdba9c1e41a2afac92.png

Calculating Attenuation \(A_{H}\) for whole domain#

\[A_{H}(r) = \frac{\left [Z_a(r) \right ]^b \cdot f(\Delta \phi_{DP})}{0.46b \int_{r1}^{r2} \left [Z_a(s) \right ]^b ds + f(\Delta \phi_{DP}) \cdot 0.46b \int_{r}^{r2} \left [Z_a(s) \right ]^b ds}\]

We can reduce the number of operations by rearranging the equation like this:

\[A_{H}(r) = \frac{\left [Z_a(r) \right ]^b}{\frac{0.46b \int_{r1}^{r2} \left [Z_a(s) \right ]^b ds}{f(\Delta \phi_{DP})} + 0.46b \int_{r}^{r2} \left [Z_a(s) \right ]^b ds}\]

iza_fdphi = iza / fdphi
display(iza_fdphi)

<xarray.DataArray (azimuth: 360, range: 750, alpha: 25)> Size: 54MB
dask.array<truediv, shape=(360, 750, 25), dtype=float64, chunksize=(180, 750, 25), chunktype=numpy.ndarray>
Coordinates: (12/17)
  * azimuth      (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation    (azimuth) float64 3kB 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
  * range        (range) float32 3kB 100.0 300.0 500.0 ... 1.497e+05 1.499e+05
    bins         (azimuth, range) float32 1MB dask.array<chunksize=(180, 375), meta=np.ndarray>
    gr           (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    ...           ...
    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:
    is_dft:     0
    _Undetect:  0.0

# integral from r1 to r2
iza_first = iza_fdphi[:, cphase.first_idx.compute()]
display(iza_first)

<xarray.DataArray (azimuth: 360, alpha: 25)> Size: 72kB
dask.array<transpose, shape=(360, 25), dtype=float64, chunksize=(360, 25), chunktype=numpy.ndarray>
Coordinates: (12/17)
  * azimuth      (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation    (azimuth) float64 3kB 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
    bins         (azimuth) float32 1kB dask.array<chunksize=(360,), meta=np.ndarray>
    gr           (azimuth) float64 3kB dask.array<chunksize=(360,), meta=np.ndarray>
    range        (azimuth) float32 1kB 5.49e+04 5.49e+04 ... 1.21e+04 1.21e+04
    ...           ...
    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:
    is_dft:     0
    _Undetect:  0.0

# Specific attenutation attributes
spec_att = dict(
    units="dB/km",
    standard_name="specific_attenuation_h",
    long_name="Specific attenuation H",
    short_name="AH",
)

ah = za / (iza_first + iza)
ah.attrs = spec_att

display(ah)

<xarray.DataArray (azimuth: 360, range: 750, alpha: 25)> Size: 54MB
dask.array<truediv, shape=(360, 750, 25), dtype=float64, chunksize=(180, 750, 25), chunktype=numpy.ndarray>
Coordinates: (12/17)
  * azimuth      (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation    (azimuth) float64 3kB 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
  * range        (range) float32 3kB 100.0 300.0 500.0 ... 1.497e+05 1.499e+05
    bins         (azimuth, range) float32 1MB dask.array<chunksize=(180, 375), meta=np.ndarray>
    gr           (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    ...           ...
    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:
    units:          dB/km
    standard_name:  specific_attenuation_h
    long_name:      Specific attenuation H
    short_name:     AH

fig = plt.figure(figsize=(8, 8))
ticks_ah = np.arange(0, 1, 0.05)
cmap = rp.get_discrete_cmap(ticks_ah, "HomeyerRainbow")
norm = rp.get_discrete_norm(ticks_ah)
im = ah.sel(alpha=0.28, method="nearest").wrl.vis.plot(
    ax=111, fig=fig, norm=norm, cmap=cmap
)

../../_images/3a98ad8736379999aa4724a96e4b94be33f29f36e7f02be7f45c68a285507027.png

\(K_{DP}\) from \(A_H\) vs. \(K_{DP}\) from \(\phi_{DP}\)#

\(K_{DP} = \frac{A_H}{\alpha}\)
\(K_{DP} = \frac{1}{2}\frac{\mathrm{d}\phi_{DP}}{\mathrm{d}r}\)

What are the benefits of \(K_{DP}(A_H)\)?

no noise artefacts
no \(\delta\)
no negative \(K_{DP}\)
no spatial degradation

filling nan with zero

kdp_ah = ah.fillna(0) / alphax
kdp_ah.attrs = swp.KDP.attrs

Calculate \(\phi_{DP}^{cal}(r, \alpha)\) for whole domain#

\[\phi_{DP}^{cal}(r_i, \alpha) = 2 \cdot \int_{r1}^{r2} \frac{A_H(s; \alpha)}{\alpha}ds\]

phical = kdp_ah.wrl.dp.phidp_from_kdp()
phical = phical.rename("PHICAL")
phical.attrs = swp.PHIDP.attrs
display(phical)

<xarray.DataArray 'PHICAL' (azimuth: 360, alpha: 25, range: 750)> Size: 54MB
dask.array<mul, shape=(360, 25, 750), dtype=float64, chunksize=(180, 25, 750), chunktype=numpy.ndarray>
Coordinates: (12/17)
  * azimuth      (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation    (azimuth) float64 3kB 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
  * alpha        (alpha) float64 200B 0.1 0.11 0.12 0.13 ... 0.31 0.32 0.33 0.34
  * range        (range) float32 3kB 100.0 300.0 500.0 ... 1.497e+05 1.499e+05
    bins         (azimuth, range) float32 1MB dask.array<chunksize=(180, 375), meta=np.ndarray>
    ...           ...
    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

Add Moments#

swp = swp.assign(
    AH=ah.sel(alpha=0.28, method="nearest"),
    KDP_AH=kdp_ah.sel(alpha=0.28, method="nearest"),
    PHIDP_CORR5=phical.sel(alpha=0.28, method="nearest"),
)
display(swp)

<xarray.Dataset> Size: 77MB
Dimensions:            (azimuth: 360, range: 750)
Coordinates: (12/17)
  * azimuth            (azimuth) float64 3kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
    elevation          (azimuth) float64 3kB 0.9998 0.9998 ... 0.9998 0.9998
    time               (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.69900...
  * range              (range) float32 3kB 100.0 300.0 ... 1.497e+05 1.499e+05
    bins               (azimuth, range) float32 1MB dask.array<chunksize=(180, 375), meta=np.ndarray>
    gr                 (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    ...                 ...
    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'
    alpha              float64 8B 0.28
Data variables: (12/49)
    CMAP               (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    CMAP2              (azimuth, range) uint8 270kB dask.array<chunksize=(180, 750), meta=np.ndarray>
    DBTH               (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    DBTV               (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    DBZH               (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    DBZV               (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    ...                 ...
    WRADH              (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    WRADV              (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    ZDR                (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    AH                 (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    KDP_AH             (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    PHIDP_CORR5        (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
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.73

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
im = rp.plot_moment(phical.sel(alpha=0.28, method="nearest"), ax=ax)

../../_images/b69ec49861f5210a75ad007b72f4f1fc152eb2d1709431d5b9ea2f4191de0982.png

Todo: check how to handle this!!!

Check the plot for different azimuth angles and start over from the beginning to iteratively find an adequate value for \(\alpha\).

An automated algorithm would try to minimize the \(\Delta\)-Sum (below).

Plotting different \(\phi_{DP}\) for Radar Domain#

fig = plt.figure(figsize=(10, 8))
fig.suptitle(swp.time.min().values.astype("<M8[s]"), fontsize=14)

ticks1 = np.arange(0, 50, 1)


norm1 = rp.get_discrete_norm(ticks1)
cmap = rp.get_discrete_cmap(ticks1, plt.cm.get_cmap("turbo"))
im1 = swp.PHIDP_OC1.wrl.vis.plot(
    ax=221, fig=fig, norm=norm1, cmap=cmap, add_colorbar=False
)
t = plt.title(r"Uncorrected $\phi_{DP}$")
t.set_y(1.1)

im2 = swp.PHIDP_CORR1.ffill("range").wrl.vis.plot(
    ax=222, fig=fig, cmap=cmap, norm=norm1, add_colorbar=False
)
t = plt.title(r"Thresholded $\phi_{DP}$")
t.set_y(1.1)

fig.subplots_adjust(wspace=0.35, hspace=0.3)

ticks2 = np.arange(-5, 11, 1)
im3 = swp.PHIDP_CORR2.ffill("range").wrl.vis.plot(
    ax=223, fig=fig, norm=norm1, cmap=cmap, add_colorbar=False
)
t = plt.title(r"Corrected $\phi_{DP}$")
t.set_y(1.1)

im4 = phical.sel(alpha=0.28, method="nearest").wrl.vis.plot(
    ax=224, fig=fig, cmap=cmap, norm=norm1, add_colorbar=False
)
t = plt.title(r"Calculated $\phi_{DP}$")
t.set_y(1.1)

# add colorbar
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.825, 0.1, 0.025, 0.8])
cbar = fig.colorbar(im1, cax=cbar_ax, norm=norm1, label=r"$\phi_{DP}$ (deg)")
cbar.locator = plt.FixedLocator(ticks1)
cbar.update_ticks()
# make some space between subplots
fig.subplots_adjust(wspace=0.3, hspace=0.3)
plt.show()

../../_images/3bb5b977a8ddceed12ef6103664c6617b5b84332ee700b90e25ee3be57f0e4a0.png

Calculation of \(\Delta\)-Sum#

\[\Delta = \sum_{i=1}^{N} \left|{\phi_{DP}^{cal}(r_i, \alpha, b)- \phi_{DP}^{m}(r_i)}\right|\]

# calculate delta sum
diff = np.abs(phical - swp.PHIDP_CORR2)
diffsel = diff.where(
    (diff.range > cphase.start_range) & (diff.range < cphase.stop_range)
)
delta = diffsel.sum("range", skipna=True)
display(delta)

<xarray.DataArray (azimuth: 360, alpha: 25)> Size: 72kB
dask.array<sum-aggregate, shape=(360, 25), dtype=float64, chunksize=(90, 25), chunktype=numpy.ndarray>
Coordinates:
  * azimuth      (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    elevation    (azimuth) float64 3kB 0.9998 0.9998 0.9998 ... 0.9998 0.9998
    time         (azimuth) datetime64[ns] 3kB 2017-08-10T11:03:42.699000 ... ...
  * alpha        (alpha) float64 200B 0.1 0.11 0.12 0.13 ... 0.31 0.32 0.33 0.34
    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

fig = plt.figure(figsize=(10, 8))
ax1 = plt.subplot(111, projection="polar")
# set the lable go clockwise and start from the top
ax1.set_theta_zero_location("N")
# clockwise
ax1.set_theta_direction(-1)
theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
ax1.plot(theta, delta, linewidth=1)  # , color="b")
t = ax1.set_title(r"$\Delta$-Sum")
t.set_y(1.1)
ax1.set_rmax(1000)
xticks = ax1.set_xticks(np.pi / 180.0 * np.linspace(0, 360, 36, endpoint=False))

../../_images/aa49b6556825368c6aaf3c23ce5b869583d3009c2bb9e45881551a7ee2d944dd.png

Apply attenuation correction#

Finally, after retrieval of the fitting \(\alpha\) we correct \(Z_H\) and \(Z_{DR}\) for attenuation.

Plot Uncorrected vs Corrected#

corr = 5
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6), (ax7, ax8, ax9)) = plt.subplots(
    3, 3, sharey=True, sharex=True, figsize=(12, 10)
)
rp.plot_moment(swp.DBTH, ax=ax1)
rp.plot_moment(swp[f"DBTH_CORR{corr}"], ax=ax2)
rp.plot_moment(swp[f"DBTH_CORR{corr}"] - swp.DBTH, ax=ax3)
rp.plot_moment(swp.DBZH, ax=ax4)
rp.plot_moment(swp[f"DBZH_CORR{corr}"], ax=ax5)
rp.plot_moment(swp[f"DBZH_CORR{corr}"] - swp.DBZH, ax=ax6)
rp.plot_moment(swp.ZDR, ax=ax7)
rp.plot_moment(swp[f"ZDR_CORR{corr}"], ax=ax8)
rp.plot_moment(swp[f"ZDR_CORR{corr}"] - swp.ZDR, ax=ax9)
plt.tight_layout()

../../_images/a16c2ad214ccb41d4777bca6ab79bb077c4327763ad13c969520c5295fdb53a3.png

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    2, 2, sharex=True, sharey=True, figsize=(12, 9)
)
rp.plot_moment(swp.KDP_AH, ax=ax1)
rp.plot_moment(swp.KDP_VULP, ax=ax2)
rp.plot_moment(swp.KDP_LANC.fillna(0), ax=ax3)
rp.plot_moment(swp.KDP, ax=ax4)
plt.tight_layout()

../../_images/59e4ba943677070df821c4dc50d5da15c0159db075d436eb51f726cb49703653.png

Save data#

dim0 = "azimuth"
dtime = swp.time.min(dim0).dt.floor("5min").dt.strftime("%Y%m%dT%H%M").values
stage = "06_zphi"
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/06_zphi_20170810T1100.zarr

<xarray.backends.zarr.ZarrStore at 0x7f3a9228b740>

ZPHI-Method

Contents