Attenuation Correction#

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

import radarpolarimetry as rp

Load Data and Plot Initial Situation#

stage = "06_zphi"
dtime = "20170810T1100"
fname = rp.data_path() / "lecture" / f"{stage}_{dtime}.zarr"
swp = xr.open_dataset(fname, chunks={})

display(swp)

<xarray.Dataset> Size: 83MB
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 dask.array<chunksize=(360,), meta=np.ndarray>
    time               (azimuth) datetime64[ns] 3kB dask.array<chunksize=(360,), meta=np.ndarray>
  * 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 ...
    crs_wkt            int64 8B ...
    latitude           float64 8B ...
    longitude          float64 8B ...
    NOISE_LEVEL        float64 8B ...
    sweep_mode         <U20 80B ...
Data variables: (12/52)
    AH                 (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    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>
    DBTH_CORR5         (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>
    ...                 ...
    VRADH              (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    VRADV              (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>
    ZDR_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, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
    nrows=2,
    ncols=2,
    figsize=(12, 10),
    sharex=True,
    sharey=True,
    squeeze=True,
    constrained_layout=True,
)

rp.plot_moment(swp.DBTH, ax=ax1)
rp.plot_moment(swp.DBZH, ax=ax2)
rp.plot_moment(swp.ZDR, ax=ax3)
rp.plot_moment(swp.PHIDP_OC1, ax=ax4)

t = fig.suptitle(
    f"sweep at {swp.sweep_fixed_angle.values:.1f} deg elevation at {swp.time.min().values.astype('<M8[s]')}",
    fontsize=18,
)

../../_images/a5b2b4693ef2014ad70ff3b4308103a41068cd55f41c54dbe73ee38702a806f1.png

Phase Processing#

Phase processing was already done in the prior lectures. So we just use the corrected phase moments from the dataset.

Exercise 1: Check different pre-calculated PHIDP values for specific radials#

use the below pre-created code
add more plots of other pre-calculated PHIDP values
also plot reflectivity values on first twin axis
add a second twin axis of host axis with label and position for RHOHV
select suiting axis limits
move legend to suitable location

azi = 105
azi = np.arange(azi, azi + 1)
fig = plt.figure(figsize=(10, 6))

ax0 = rp.create_lineplot(fig, subplot=111, ylabel=r"$PHI_{DP}$ (deg)")
t = ax0.set_title(f"{swp.time.min().values.astype('<M8[s]')}")
t.set_y(1.05)

line = ax0.plot(swp.PHIDP[azi, :].T, "darkblue", label="$PHI_{DP}^{raw}$")
line = ax0.plot(swp.PHIDP_OC1[azi, :].T, "b", lw=0.5, label="$PHI_{DP}OC^{corr1}$")
line = ax0.plot(swp.PHIDP_CORR1[azi, :].T, "b", label="$PHI_{DP}^{corr1}$")

ax0.set_ylim(-10, 120)

plt.legend(ncol=4, loc=8)

# add twin axis of host axis with label and position
ax1 = rp.add_axis(ax0, ylabel=r"$UZ_H$ (dBZ)", pos=1.0)

# add another twin axis of host axis with label and position

plt.tight_layout()

../../_images/a76b410cd2d995f4d82d1018e35c9728e8c2a15eaa0a4b44f1b2a727d854ebbc.png

Exercise 2: Display PHIDP moments in 2by2 plot#

select different moments
change resolution

Standard Correction Method#

\[ Z^{corr}_{H}(r) = Z^{att}_{H}(r) + \alpha\phi_{DP}(r) \]

\[ Z^{corr}_{DR}(r) = Z^{att}_{DR}(r) + \beta\phi_{DP}(r) \]

\[ \alpha_X = 0.28 \]

\[ \beta_X = 0.05 \]

Apply ZH and ZDR standard correction#

This also adds the corrected moments to the sweep.

uh = "DBTH"
zh = "DBZH"
zdr = "ZDR"
phi_list = [f"PHIDP_CORR{i}" for i in range(1, 5)]
alphax = 0.28
betax = 0.05

with xr.set_options(keep_attrs=True):
    for phi in phi_list:
        uhcorr = swp[uh] + alphax * swp[phi]
        uhcorr.name = f"{swp[uh].name}_{swp[phi].name.split('_')[-1]}"
        zhcorr = swp[zh] + alphax * swp[phi]
        zhcorr.name = f"{swp[zh].name}_{swp[phi].name.split('_')[-1]}"
        zdrcorr = swp[zdr] + betax * swp[phi]
        zdrcorr.name = f"{swp[zdr].name}_{swp[phi].name.split('_')[-1]}"
        swp = swp.assign(
            {uhcorr.name: uhcorr, zhcorr.name: zhcorr, zdrcorr.name: zdrcorr}
        )

ZPHI Method#

See ZPHI Method - Theory Notebook.

Todo: use precalculated data from zphi notebook! 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

zh = swp.DBZH.where((swp.CMAP2 == 0) & valid & (swp.RHOHV > 0.9))

Calculate Differential Attenuation#

ah = rp.zphi(phi, zh, rng=2000)
display(ah)

range res [km]: 0.2

<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/618938f224deef9728b6d1b53620e818fbb23b8379be9d93b0a5ef5c9a8232bd.png

KDP from AH#

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

<matplotlib.collections.QuadMesh at 0x7f65cb1a0550>

../../_images/e47924b68edd4fb66fcac497d8e58fcdb6779e3baa12010f03704207ebbce27b.png

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: 109MB
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>
    ...                 ...
    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/64)
    AH                 (azimuth, range) float64 2MB dask.array<chunksize=(180, 750), meta=np.ndarray>
    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>
    DBTH_CORR5         (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>
    ...                 ...
    DBTH_CORR3         (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    DBZH_CORR3         (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    ZDR_CORR3          (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    DBTH_CORR4         (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    DBZH_CORR4         (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), meta=np.ndarray>
    ZDR_CORR4          (azimuth, range) float64 2MB dask.array<chunksize=(90, 375), 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

Apply ZH and ZDR ZPHI-based correction#

Plot KDP moments overview#

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

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

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/58091b63fc06d0ec60d6e22414578cc2dbccf8e386aa1a7761d3601508a6d19a.png

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/6fcb44a9886eb53c285f4156615fc40f610d79e8949b526a971f3018a2565426.png

Difference Plots#

zdiff0 = "DBTH_CORR5"
zdiff1 = "DBTH_CORR4"
zdrdiff0 = "ZDR_CORR5"
zdrdiff1 = "ZDR_CORR4"

zdiff = swp[zdiff0] - swp[zdiff1]
zdrdiff = swp[zdrdiff0] - swp[zdrdiff1]

fig = plt.figure(figsize=(10, 5))
fig.suptitle(zdiff.time.min().values.astype("<M8[s]"), fontsize=14)
cmap = "HomeyerRainbow"
im = zdiff.wrl.vis.plot(
    ax=121, fig=fig, cmap=cmap, vmin=-2, vmax=2, levels=np.arange(-1, 4, 0.1)
)
t = plt.title("Difference Reflectivity")
t.set_y(1.1)

im = zdrdiff.wrl.vis.plot(
    ax=122, fig=fig, cmap=cmap, vmin=0, vmax=2, levels=np.arange(-1, 2, 0.1)
)
t = plt.title("Difference ZDR")
t.set_y(1.1)

fig.tight_layout()

../../_images/12437c0bf207890fc256b749642420a59838eebac4e63c8077cb391b15fdd497.png

Save data#

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

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

Solution 1: Check different pre-calculated PHIDP values for specific radials#

azi = 104
azi = np.arange(azi, azi + 1)
fig = plt.figure(figsize=(10, 6))
ax0 = rp.create_lineplot(fig, subplot=111, ylabel=r"$PHI_{DP}$ (deg)")
t = ax0.set_title(f"{swp.time.min().values.astype('<M8[s]')}")
t.set_y(1.05)
line = ax0.plot(swp.PHIDP[azi, :].T, "b.", label="$PHI_{DP}^{raw}$")
line = ax0.plot(swp.PHIDP_OC1[azi, :].T, "b", lw=0.5, label="$PHI_{DP}OC^{corr1}$")
line = ax0.plot(swp.PHIDP_OC2[azi, :].T, "r", lw=0.5, label="$PHI_{DP}OC^{corr2}$")
line = ax0.plot(swp.PHIDP_OC3[azi, :].T, "k", lw=0.5, label="$PHI_{DP}OC^{corr3}$")
line = ax0.plot(swp.PHIDP_OC4[azi, :].T, "g", lw=0.5, label="$PHI_{DP}OC^{corr4}$")
line = ax0.plot(swp.PHIDP_CORR1[azi, :].T, "b", label="$PHI_{DP}^{corr1}$")
line = ax0.plot(swp.PHIDP_CORR2[azi, :].T, "r", label="$PHI_{DP}^{corr2}$")
line = ax0.plot(swp.PHIDP_CORR3[azi, :].T, "k", label="$PHI_{DP}^{corr3}$")
line = ax0.plot(swp.PHIDP_CORR4[azi, :].T, "g", label="$PHI_{DP}^{corr4}$")

ax0.set_ylim(-10, 40)

plt.legend(ncol=1, loc=8, bbox_to_anchor=(1.3, 0.0))

# add twin axis of host axis with label and position
ax1 = rp.add_axis(ax0, ylabel=r"$UZ_H$ (dBZ)", pos=1.0)
line = ax1.plot(swp.DBTH[azi, :].T, "g", zorder=0)

# add another twin axis of host axis with label and position
ax2 = rp.add_axis(ax0, ylabel=r"$RHO_{HV}$", pos=1.1)
line = ax2.plot(swp.RHOHV[azi, :].T, "c")
# plt.tight_layout()
plt.show()

../../_images/c3c9753f6c6dffd1dceaa46f3c83441ef62276e1f889a047e83aab9435e3a8f9.png

Solution 2: Display PHIDP moments in 2by2 plot#

select different moments
change resolution

ticks = np.arange(0, 52, 2)
moments = ["PHIDP_OC1", "PHIDP_CORR1", "PHIDP_OC4", "PHIDP_CORR4"]

rp.plot_2by2(swp, moments, ticks=ticks, add_colorbar=False)

# add colorbar
fig = plt.gcf()
im = fig.axes[0].collections[0]
fig.subplots_adjust(right=0.7)
cbar_ax = fig.add_axes([0.75, 0.1, 0.025, 0.8])
cbar = fig.colorbar(im, cax=cbar_ax, label=r"$\phi_{DP}$ (deg)")
cbar.locator = plt.FixedLocator(ticks)
cbar.update_ticks()

# make some space between subplots
fig.subplots_adjust(wspace=0.25, hspace=0.25)

../../_images/e8f8b54f399864476d1179f4a323bfe275afe9dbef60df2e9692e42f1d3032d8.png

Attenuation Correction

Contents