System Differential Phase \(\phi_{DP}^{sys}\)#

Correct retrieval of \(\phi_{DP}^{sys}\) (Offset) is crucial for correct processing of further derivations. Normally \(\phi_{DP}^{sys}\) is more or less constant. It can have azimuthal/elevational deviations depending on the antenna near field.

Retrieval of \(\phi_{DP}^{sys}\) is sometimes tedious, due to the contamination of the signal with clutter and other artifacts.

Import Section#

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

import radarpolarimetry as rp

Open Dataset#

Use data from prior lectures if available.

stage = "01_pre"
dtime = "20170810T1100"
fname = rp.data_path() / "lecture" / f"{stage}_{dtime}.zarr"
swp = xr.open_dataset(fname)

# start_time = dt.datetime(2017, 8, 10, 11)
# scan = "n_ppi_010deg"
# swp = rp.load_single_sweep(start_time, scan)

Inspect Dataset#

display(swp)

<xarray.Dataset> Size: 47MB
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/22)
    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 ...
    ...                 ...
    sweep_number       int64 8B ...
    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.73

System Differential Phase \(\phi_{DP}^{sys}\) via first precipitating bins#

This is a most common algorithm. But there are several approaches.

N consecutive radar bins with rhoHV > threshold
maximum number of valid bins in a N-size window
first N valid bins (not necessarily consecutive)

Mask source data#

Use CMAP moment, if available

Value	Bit	Name	Description
0x0001	0	NOISE-H	Is set when the horizontal SNR is below NOISEThreshold
0x0002	1	NOISE-V	Is set when the vertical SNR is below NOISEThreshold
0x0004	2	CCOR-H	Is set when horizontal CCOR is lower than CCORThreshold
0x0008	3	CCOR-V	Is set when vertical CCOR is lower than CCORThreshold
0x0010	4	SQI1-H	Is set when the horizontal SQI is lower than SQI1Threshold
0x0020	5	SQI1-V	Is set when the vertical SQI is lower than SQI1Threshold
0x0040	6	SQI2-H	Is set when the horizontal SQI is lower than SQI2Threshold
0x0080	7	SQI2-V	Is set when the vertical SQI is lower than SQI2Threshold
0x0100	8	SIGPOW-H	Is set when the horizontal SIGPOW is lower than SIGPOWThreshold
0x0200	9	SIGPOW-V	Is set when the vertical SIGPOW is lower than SIGPOWThreshold
0x0400	10	URHOHV	Is set when UrhoHV is lower than RHOHVThreshold
0x0800	11	CMS	Censored by Clutter Micro Suppression

If CMAP is not available we’ll only rely on CMAP2 from Notebook Radar Clutter and Noise. Otherwise we combine CMAP/CMAP2.

CMAP = swp.CMAP.astype("uint16")
bit = 0
cmp = (CMAP >> bit) & 0x1
cmp.wrl.vis.plot()
plt.gca().set_title(f"CMAP bit {bit}")

Text(0.5, 1.0, 'CMAP bit 0')

../../_images/76fa898d7278280704b705c655c74a61d76b9874fba3143c9a4c060f37939580.png

Exercise 1: Check whether other CMAP bits are set and plot#

Use plt.subplots and plot all bits in in one figure.

Exercise 2: Check CMAP2 and plot#

Use plt.subplots and plot everything in one figure.

mask phase#

no_clutter_bb = (cmp == 0) & (swp.CMAP2 == 0)
phimask = swp.PHIDP.where(no_clutter_bb & (swp.RHOHV >= 0.9))
# create binary array
phib = xr.where(np.isnan(phimask), 0, 1)

Exercise 3: Change range and process again#

decrease and/or increase range to see how this affects the estimated system differential phase

Exercise 4: Change range and process again#

decrease and/or increase range to see how this affects the estimated system differential phase

Exercise 5: Change number of valid bins and process again#

decrease and/or increase number of valid bins to see how this affects the estimated system differential phase.

System Differential Phase \(\phi_{DP}^{sys}\) via phase histogram#

This is a quite new technique. The idea behind is:

\(\phi_{DP}\) is constantly increasing
\(\phi_{DP}\) is inherently noisy

That means, the majority of phase measurements (precipitating bins) should lie around \(\phi^{sys}_{DP}\). It’s extremely robust since it does not rely on finding precipitating bins in each ray or the like.

from xhistogram.xarray import histogram as xhist

hlist = []
phase_res = 0.1
bins = np.arange(-180, 180, phase_res)

phi_hist = xhist(swp.PHIDP, dim=("range",), bins=[bins])
display(phi_hist)

<xarray.DataArray 'histogram_PHIDP' (azimuth: 360, PHIDP_bin: 3599)> Size: 10MB
array([[1, 0, 0, ..., 2, 2, 0],
       [0, 0, 0, ..., 1, 0, 0],
       [1, 1, 1, ..., 0, 1, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [1, 0, 0, ..., 1, 1, 1],
       [0, 0, 1, ..., 1, 0, 1]], shape=(360, 3599))
Coordinates:
  * azimuth    (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5
  * PHIDP_bin  (PHIDP_bin) float64 29kB -179.9 -179.9 -179.8 ... 179.7 179.8

Summing Plot - Sweep Histogram#

fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 5))
phi_hist.sum("azimuth").plot(ax=ax1)
ax1.grid()
phi_hist.sum("azimuth").plot(ax=ax2)
ax2.set_xlim(80, 100)
ax2.grid()

../../_images/2dfccbafaf0f6a6b5448062b559472c0fb87a1b3afd2b11e20ffe7074610f437.png

Exercise 6: Calculate Histogram with masked data#

what do you expect and what does actually happen

Smoothing, Max and Turning Point#

smooth with window size 5 (0.5 deg) and get phase at maximum
cut by +/- 5deg
normalize with maximum
differentiate, smooth with window size 11 (1.1 deg) and get phase at turning point

phi_hist_smooth = phi_hist.rolling(PHIDP_bin=5, center=True).mean()
pmax = phi_hist_smooth.idxmax("PHIDP_bin")
phi_hist_smooth_cut = phi_hist_smooth.where(
    (phi_hist.PHIDP_bin > (pmax - 5)) & (phi_hist.PHIDP_bin < (pmax + 5)), drop=True
)
phi_hist_smooth_cut = phi_hist_smooth_cut / phi_hist_smooth_cut.max()
sysphi_4 = (
    phi_hist_smooth_cut.differentiate("PHIDP_bin")
    .rolling(PHIDP_bin=11, center=True)
    .mean()
    .idxmax("PHIDP_bin")
)

pmax.plot()
sysphi_4.plot()

[<matplotlib.lines.Line2D at 0x7f9a18e09450>]

../../_images/48f5d64392699dfdfc4392d640cab6e43eab2192633ac20d3d1a58b4f7426431.png

Overview and Diagnostic plots#

add theta in radians to dataset

theta = np.linspace(0, 2 * np.pi, num=360, endpoint=False)
phi_hist = phi_hist.assign_coords(
    theta=(["azimuth"], theta, {"standard_name": "azimuth angle"})
)
phi_1m = phi_1m.assign_coords(
    theta=(["azimuth"], theta, {"standard_name": "azimuth angle"})
)
phi_2m = phi_2m.assign_coords(
    theta=(["azimuth"], theta, {"standard_name": "azimuth angle"})
)
phi_3m = phi_3m.assign_coords(
    theta=(["azimuth"], theta, {"standard_name": "azimuth angle"})
)

fig = plt.figure(figsize=(12, 12))
vmin = 80
vmax = 90
startaz = swp.sortby("time").azimuth[0]

ax1 = rp.get_polar_ax(fig=fig, ax=111)

ax1.axvline(np.radians(startaz.values), c="black", lw=1.0)
phi_1m.plot(x="theta", c="k", lw=0.5, ax=ax1)
phi_2m.plot(x="theta", c="m", lw=0.5, ax=ax1)
phi_3m.plot(x="theta", c="r", lw=0.5, ax=ax1)
ax1.axhline(sysphi_1, c="k").get_path()._interpolation_steps = 180
ax1.axhline(sysphi_2, c="m").get_path()._interpolation_steps = 180
ax1.axhline(sysphi_3, c="r").get_path()._interpolation_steps = 180
ax1.plot(theta, sysphi_4, c="b")

phi_hist.where(phi_hist > 3).plot(
    ax=ax1, x="theta", cmap="HomeyerRainbow", levels=np.arange(0, 30, 2), zorder=1
)
ax1.set_ylim(vmin, vmax)
ax1.grid()
plt.tight_layout()

../../_images/4938bb30c226cded06ba507ce48583c2a0608b9b4e9d989246b2026e2e8e893a.png

Add and save data#

swp = swp.assign({"SYSPHI1": sysphi_1})
swp = swp.assign({"PHIDP_OC1": swp.PHIDP - sysphi_1})
swp.PHIDP_OC1.attrs = swp.PHIDP.attrs
swp = swp.assign({"SYSPHI2": sysphi_2})
swp = swp.assign({"PHIDP_OC2": swp.PHIDP - sysphi_2})
swp.PHIDP_OC2.attrs = swp.PHIDP.attrs
swp = swp.assign({"SYSPHI3": sysphi_3})
swp = swp.assign({"PHIDP_OC3": swp.PHIDP - sysphi_3})
swp.PHIDP_OC3.attrs = swp.PHIDP.attrs
swp = swp.assign({"SYSPHI4": sysphi_4})
swp = swp.assign({"PHIDP_OC4": swp.PHIDP - sysphi_4})
swp.PHIDP_OC4.attrs = swp.PHIDP.attrs

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

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

Solution 1: Check whether other CMAP bits are set and plot#

Use plt.subplots and plot everything in one figure.

fig, axes = plt.subplots(3, 4, sharex=True, sharey=True, figsize=(16, 10))
for bit in range(12):
    cmp = (CMAP >> bit) & 0x1
    ax = axes.flat[bit]
    cmp.wrl.vis.plot(ax=ax)
    ax.set_title(f"CMAP bit {bit}")
plt.tight_layout()

../../_images/523a2699d22c342d311f4b3efa3fb9a0199ea38c4b3a55a9bc016c3f973eb78a.png

Solution 2: Check CMAP2 and plot#

Use plt.subplots and plot everything in one figure.

fig, axes = plt.subplots(1, 2, sharey=True, figsize=(12, 10))
for bit in range(2):
    cmp = (swp.CMAP2 >> bit) & 0x1
    ax = axes.flat[bit]
    cmp.wrl.vis.plot(ax=ax)
    ax.set_title(f"CMAP2 bit {bit}")
plt.tight_layout()

../../_images/224eb5f48292fefb60f7e9e2e895b6d717c0565725c62743231e555d9d6a0a8a.png

Solution 3: Change range and process again#

decrease and/or increase range to see how this affects the estimated system differential phase

rng = 1000.0
N = swp.pipe(rp.get_bins_from_range, rng)
print(N)

# take N range bins and calculate sum, for finding N consecutive bins
phib_sum = phib.rolling(range=N, center=True).sum(skipna=True)

# get maximum along range
smax = phib_sum.max(dim="range", skipna=True)

# if smax == N, we've identified the first consecutive range of valid phase values
# derive start_range and stop_range
rstep = rp.get_range_step(phib_sum)
start_range = phib_sum.idxmax(dim="range").where(smax == N) - smax // 2 * rstep
stop_range = start_range + N * rstep

Solution 4: Change range and process again#

decrease and/or increase range to see how this affects the estimated system differential phase

rng = 3000.0
N = swp.pipe(rp.get_bins_from_range, rng)
print(N)

# take N range bins and calculate sum, for finding N consecutive bins
phib_sum = phib.rolling(range=N, center=True).sum(skipna=True)

# get the maximum along the ray
# derive start_range and stop_range
rstep = rp.get_range_step(phib_sum)
start_range = phib_sum.idxmax(dim="range") - N // 2 * rstep
stop_range = start_range + N * rstep

Solution 5: Change number of valid bins and process again#

decrease and/or increase number of valid bins to see how this affects the estimated system differential phase.

N = 5

Solution 6: Calculate Histogram with masked data#

what do you expect and what does actually happen

phi_histm = xhist(phimask, dim=("range",), bins=[bins])
display(phi_histm)

<xarray.DataArray 'histogram_PHIDP' (azimuth: 360, PHIDP_bin: 3599)> Size: 10MB
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], shape=(360, 3599))
Coordinates:
  * azimuth    (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5
  * PHIDP_bin  (PHIDP_bin) float64 29kB -179.9 -179.9 -179.8 ... 179.7 179.8

fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 5))
phi_hist.sum("azimuth").plot(ax=ax1)
phi_histm.sum("azimuth").plot(ax=ax1)
ax1.grid()
phi_hist.sum("azimuth").plot(ax=ax2)
phi_histm.sum("azimuth").plot(ax=ax2)
ax2.set_xlim(80, 100)
ax2.grid()

../../_images/d17252a3fdfe029ff077fa514b00e2e097097827cbe65b643857a7aaf40baf36.png

Footer#

We’ve had a look into retrieval of System Differential Phase. Four methods have been shown and pros and cons have been discussed.

System Differential Phase \phi_{DP}^{sys}

Contents

System Differential Phase \(\phi_{DP}^{sys}\)#

Import Section#

Open Dataset#

Inspect Dataset#

System Differential Phase \(\phi_{DP}^{sys}\) via first precipitating bins#

Mask source data#

Exercise 1: Check whether other CMAP bits are set and plot#

Exercise 2: Check CMAP2 and plot#

mask phase#

N consecutive valid radar bins#

Exercise 3: Change range and process again#

maximum number of valid radar bins in N-sized window#

Exercise 4: Change range and process again#

first N valid bins#

Exercise 5: Change number of valid bins and process again#

System Differential Phase \(\phi_{DP}^{sys}\) via phase histogram#

Summing Plot - Sweep Histogram#

Exercise 6: Calculate Histogram with masked data#

Smoothing, Max and Turning Point#

Overview and Diagnostic plots#

Add and save data#

Solution 1: Check whether other CMAP bits are set and plot#

Solution 2: Check CMAP2 and plot#

Solution 3: Change range and process again#

Solution 4: Change range and process again#

Solution 5: Change number of valid bins and process again#

Solution 6: Calculate Histogram with masked data#

Footer#