Differential Phase Processing#

import matplotlib.pyplot as plt
import xarray as xr

import radarpolarimetry as rp

Differential Phase Overview#

Load Dataset (Overview)#

Showing different problems with differential phase. We use 4.5° scan to show bright band backscattering.

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

Calculate Median#

medwin = 21
phimed = phimask.pipe(rp.filter_data, medwin)
phimed.attrs = phimask.attrs
fig = plt.figure()
ax = fig.add_subplot()
rp.plot_moment(phimed, ax=ax)

<matplotlib.collections.QuadMesh at 0x7f3a9c4df9d0>

../../_images/98bd01062b0c903ebb29394bfc09bc3364fec21f661938b3c2490a7f0683e3f2.png

Plot moments#

swpx = swpm.sel(range=slice(0, 100e3))
rp.plot_2by2(swpx, moments=["DBTH", "PHIDP_OC1", "ZDR", "RHOHV"])

../../_images/e3a81ad38a332d7a99e89a81575cf2b45493400204c0c4b731a2387a79480704.png

Visualize single radial#

azi = 200

swpx = swpm.sel(range=slice(0, 80e3))
phimedx = phimed.sel(range=slice(0, 80e3))
fig = plt.figure(figsize=(10, 8))

ax = fig.add_subplot(232)
rp.plot_moment(swpx.DBZH, ax=ax)
az = swpx.isel(azimuth=azi)
plt.gca().plot(az.x, az.y, "k", lw=0.85)


# create host axis with ylabel
ax0 = rp.create_lineplot(fig, subplot=212, ylabel=r"$PHI_{DP}$ (deg)")
ax0.set_ylim(-25, 50)

# plot line on first axis
line = phimedx[azi].plot(c="b", marker=".", label=r"smoothed $PHI_{DP}$ (deg)", ax=ax0)
line = swp.PHIDP_OC1[azi].plot(c="b", lw=0.5, label=r"raw $PHI_{DP}$ (deg)", ax=ax0)

# add twin axis of host axis with label and position
ax1 = rp.add_axis(ax0, ylabel=r"$RHO_{HV}$", pos=1.0)
ax1.grid()
# add another twin axis of host axis with label and position
ax2 = rp.add_axis(ax0, ylabel=r"$UZ_H$ (dBZ)", pos=1.1)
# plot line on second axis
# line = swpx.RHOHV_NC[azi].plot(c="r", label=r"$RHO_{HV}$ - nc", ax=ax1)
line = swpx.RHOHV[azi].plot(c="r", label=r"$RHO_{HV}$", lw=1.0, ax=ax1)
line = swp.RHOHV[azi].plot(c="r", label=r"$RHO_{HV}$", lw=0.25, ax=ax1)
# plot line on third axis
line = swpx.DBZH[azi].plot(c="g", label=r"$UZ_H$ (dBZ)", ax=ax2)
line = swp.DBZH[azi].plot(c="g", label=r"$UZ_H$ (dBZ)", lw=0.25, ax=ax2)

# set title
ax0.set_title(f"{swp.time.min().values} Azi {azi}")
ax1.set_title("")
ax2.set_title("")

# legend
h1, l1 = ax0.get_legend_handles_labels()
h2, l2 = ax1.get_legend_handles_labels()
h3, l3 = ax2.get_legend_handles_labels()
ax1.legend(h1 + h2 + h3, l1 + l2 + l3, loc="lower right")

<matplotlib.legend.Legend at 0x7f3a9d9f8ad0>

../../_images/8dc1f3b3b75bb33b6b4f9ee5f66c5a26779898bb4361c3d39575b1541d6c1499.png

Exercise 1: Use other azimuths and compare#

Go back to cell above and check other azimuths and compare. Also play with different masking thresholds.

Exercise 2: Play with medwin, kwidth and sigma#

Go back above and play with different values for smoothing.

2d Filtering using xr_rolling#

Can clean out outliers of phase in single rays.

%%time
window = 41
window2 = 3
phi_median = swp_masked.PHIDP_OC1.copy().pipe(
    rp.xr_rolling, window, window2=window2, method="median", skipna=True, min_periods=3
)
# gaussian convolution 1d - smoothing
# play with sigma and kwidth
kwidth = 41
sigma = 10
gkern = rp.gauss_kernel(kwidth, sigma)

phi_clean = phi_median.pipe(rp.smooth_data, gkern)
phi_clean.attrs = phimed.attrs
phi_clean.attrs = swp_masked.PHIDP_OC1.attrs

CPU times: user 1.51 s, sys: 402 ms, total: 1.92 s
Wall time: 1.92 s

time = swp.time.min().values
# create figure
fig = plt.figure(figsize=(10, 8))
# create host axis with ylabel
ax0 = rp.create_lineplot(fig, subplot=111, ylabel=r"$PHI_{DP}$ (deg)")
# set title
ax0.set_title(f"{time} Azi {azi}")
# plot line on first axis
line = swp.PHIDP_OC1[azi].plot(c="lightblue", label=r"raw $PHI_{DP}$ (deg)", ax=ax0)
line = swp_masked.PHIDP_OC1[azi].plot(
    c="b", marker="*", label=r"masked $PHI_{DP}$ (deg)", ax=ax0
)


# add twin axis of host axis with label and position
ax1 = rp.add_axis(ax0, ylabel=r"$RHO_{HV}$", pos=1.0)
# add another twin axis of host axis with label and position
ax2 = rp.add_axis(ax0, ylabel=r"$Z_H$ (dBZ)", pos=1.1)
# plot line on second axis
line = swp_masked.RHOHV_SMOOTH[azi].plot(c="r", label=r"$RHO_{HV}$", zorder=0, ax=ax1)
# plot line on third axis
line = swp_masked.DBTH[azi].plot(c="g", label=r"$Z_H$ (dBZ)", zorder=0, ax=ax2)
ax0.set_ylim(-5, 40)
line = phiclean[azi].plot(
    c="k", marker=".", label=r"clean $PHI_{DP}$ (deg)", zorder=10, ax=ax0
)
line = phi_clean[azi].plot(
    c="y", marker=".", label=r"clean 2D $PHI_{DP}$ (deg)", zorder=10, ax=ax0
)
h1, l1 = ax0.get_legend_handles_labels()
h2, l2 = ax1.get_legend_handles_labels()
h3, l3 = ax2.get_legend_handles_labels()

ax1.legend(h1 + h2 + h3, l1 + l2 + l3)

plt.show()

../../_images/e96948d6ca87bc188e512d40dfd378b5238444b1df9450196afe1f4ab63b18fc.png

Exercise 3: Play with window sizes#

Go back above and play with different window sizes and smoothing parameters.

Overview Plot#

fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(12, 10))
axes = axes.flat
phimed.wrl.vis.plot(vmin=-5, vmax=40, ax=axes[0])
phi_median.wrl.vis.plot(vmin=-5, vmax=40, ax=axes[1])
phiclean.wrl.vis.plot(vmin=-5, vmax=40, ax=axes[2])
phi_clean.wrl.vis.plot(vmin=-5, vmax=40, ax=axes[3])
plt.tight_layout()

../../_images/af8a0e32e23a02022cb73513b7769621bc2d8f95a4ea7e9fa9cd0774def92f7a.png

Add and save data#

swp = swp.assign({"PHIDP_CORR1": phimed})
swp = swp.assign({"PHIDP_CORR2": phi_median})
swp = swp.assign({"PHIDP_CORR3": phiclean})
swp = swp.assign({"PHIDP_CORR4": phi_clean})
dim0 = "azimuth"
dtime = swp.time.min(dim0).dt.floor("5min").dt.strftime("%Y%m%dT%H%M").values
stage = "03_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/03_phidp_20170810T1100.zarr

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

Differential Phase Processing

Contents

Differential Phase Processing#

Differential Phase Overview#

Load Dataset (Overview)#

Inspect Dataset#

Calculate Median#

Plot moments#

Visualize single radial#

Exercise 1: Use other azimuths and compare#

Phase Processing Workflow#

Load dataset (workflow)#

Mask data using DBZH, RHOHV, CMAP/CMAP2#

“Blend” RHOHV#

Filtering and Smoothing#

Median Filtering#

Smoothing#

Exercise 2: Play with medwin, kwidth and sigma#

2d Filtering using xr_rolling#

Exercise 3: Play with window sizes#

Overview Plot#

Add and save data#

Footer#