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

We will use the simple method of finding the first and the last non NAN values per ray from. This is the most simple and probably not very robust method. See also PHIDP Processing Notebook.

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"Uncorrected $Z_{H}$")
t.set_y(1.1)
fig.suptitle(swp.time.min().values.astype("<M8[s]"), fontsize=14)
fig.tight_layout()

../../_images/730ca1e5c1ff78a5d303b0dbe8ff1552ebea5ce3c3defa59c0e80028477bb6bb.png

Choose Differential Phase Moment#

We also mask using CMAP/CMAP2 and RHOHV.

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

Set Range Region#

rng = 2000.0

Get Range parameters#

range_step = np.diff(phi.range)[0]
nprec = int(rng / range_step)
if nprec % 2 == 0:
    nprec += 1

offset = nprec // 2 * np.diff(phi.range)[0]
offset_idx = nprec // 2
print(nprec, offset, offset_idx)

11 1000.0 5

Create binary array#

phib = xr.where(np.isnan(phi), 0, 1)
phib.wrl.vis.plot()

<matplotlib.collections.QuadMesh at 0x7f9d5831d590>

../../_images/08fb3c0feb92840bb979c6b1184855108935603e70627f52b568bff2ae93cef3.png

take nprec range bins and calculate rolling sum#

phib_sum = phib.rolling(range=nprec, center=True, min_periods=nprec).sum(skipna=True)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
phib_sum.plot(x="azimuth", ax=ax1)
secs = "<M8[s]"
t0 = phi.time.min().values.astype(secs)
ax1.set_title(f"Valid Count - {t0}")
phi.where(phib_sum == nprec).plot(x="azimuth", ax=ax2)
ax2.set_title(f"{phi.name} - {t0}")
plt.tight_layout()

../../_images/ec86bf73830a06b3d8163cbb191f31bafef6fa1e18cf840143652386af0a9868.png

Get start range per each ray and retrieve phase values#

First consecutive valid N bins.

start_range = phib_sum.idxmax(dim="range") - offset
start_range_idx = phib_sum.argmax(dim="range") - offset_idx
first = phi.where(
    (phi.range >= start_range.load()) & (phi.range <= start_range.load() + rng),
)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8))
phi.where(phib_sum >= offset_idx).plot(x="azimuth", vmin=80, vmax=90, ax=ax1)
start_range.plot(lw=0.8, c="k", ax=ax1)
ax1.set_ylim(0, 60000)
secs = "<M8[s]"
t0 = phi.time.min().values.astype(secs)
ax1.set_title(f"{phi.name} - {t0}")
first.plot(x="azimuth", ax=ax2)
start_range.plot(lw=0.8, c="k", ax=ax2)
ax2.set_title(f"{phi.name} first N - {t0}")
ax2.set_ylim(0, 60000)
plt.tight_layout()

../../_images/cfee0384f9ce29e3174ea17ae541b6cb3b06b8f0329cb8aa9ab09b1cfa799f50.png

Get stop range per each ray and retrieve phase values#

We just do this in reversed range order. Last consecutive valid N bins.

stop_range = phib_sum[..., ::-1].idxmax(dim="range") - offset
stop_range_idx = (
    len(phib_sum.range) - (phib_sum[..., ::-1].argmax(dim="range") - offset_idx) - 2
)
last = phi.where(
    (phi.range >= stop_range.load() - rng) & (phi.range <= stop_range.load()),
)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 8))
phi.where(phib_sum >= offset_idx).plot(x="azimuth", ax=ax1)
stop_range.plot(lw=0.8, c="k", ax=ax1)
secs = "<M8[s]"
t0 = phi.time.min().values.astype(secs)
ax1.set_title(f"{phi.name} - {t0}")
last.plot(x="azimuth", ax=ax2)
stop_range.plot(lw=0.8, c="k", ax=ax2)
ax2.set_title(f"{phi.name} first N - {t0}")
plt.tight_layout()

../../_images/3f93545a069fa41a3be13bf2d1576266605ffbc7a1872d4e373eeb4bcefc32fd.png

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

first = first.median("range", skipna=True)
last = last.median("range", skipna=True)
dphi = last - first

Merge sweep Data#

cphase = xr.Dataset(
    dict(
        phib=phib_sum,
        offset=offset,
        offset_idx=offset_idx,
        start_range=start_range,
        stop_range=stop_range,
        first=first,
        first_idx=start_range_idx,
        last=last,
        last_idx=stop_range_idx,
        dphi=dphi,
    )
)
swp = xr.merge([swp, cphase])
display(swp)

<xarray.Dataset> Size: 71MB
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 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>
    ...                 ...
    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'
Data variables: (12/46)
    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>
    ...                 ...
    stop_range         (azimuth) float32 1kB 5.51e+04 5.53e+04 ... 6.19e+04
    first              (azimuth) float64 3kB dask.array<chunksize=(180,), meta=np.ndarray>
    first_idx          (azimuth) int64 3kB dask.array<chunksize=(180,), meta=np.ndarray>
    last               (azimuth) float64 3kB dask.array<chunksize=(180,), meta=np.ndarray>
    last_idx           (azimuth) int64 3kB dask.array<chunksize=(180,), meta=np.ndarray>
    dphi               (azimuth) float64 3kB dask.array<chunksize=(180,), 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

\(\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()
cphase.dphi.plot(ls="-", marker=".", label="delta", ax=ax2)
ax2.grid()
ax2.legend()
plt.tight_layout()

../../_images/50610fc8f99bda12e78f895489d0f942ddd2b2e555f35de03133af815c58a287.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")
ax4 = plt.subplot(234, projection="polar")
ax5 = plt.subplot(235, projection="polar")
ax6 = plt.subplot(236, 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")
ax4.set_theta_zero_location("N")
ax5.set_theta_zero_location("N")
ax6.set_theta_zero_location("N")
# clockwise
ax1.set_theta_direction(-1)
ax2.set_theta_direction(-1)
ax3.set_theta_direction(-1)
ax4.set_theta_direction(-1)
ax5.set_theta_direction(-1)
ax6.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")

ax4.plot(theta, cphase.start_range, color="b", linewidth=3)
_ = ax4.set_title("Delta PHIDP - First Range")
ax5.plot(theta, cphase.stop_range, color="r", linewidth=3)
_ = ax5.set_title("Delta PHIDP - Last Range")
ax6.plot(theta, cphase.dphi, color="g", linewidth=3)
_ = ax6.set_title("Delta PHIDP")

../../_images/06e519d09411cb91af2d09d1fa2e502943c12da666a8f2786138f93f6d60e612.png

Save data#

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

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

\Delta \phi_{DP}

Contents