Attenuation correction module

Attenuation due to precipitation plays a major role in underestimating the precipitation. Therefore, an attenuation correction according to the ZPHI method (Testud et. al, 2001) is performed.

Now we start step by step with the attenuation correction for clutter corrected data of the weather radar FURUNO.

import wrainfo as wrf

At first, we have to import following packages.

[1]:
import datetime as dt
import wradlib as wrl
import warnings
warnings.filterwarnings('ignore')

For the individual steps of the clutter correction, see the Clutter detection module. We continue working with the created dataset at the end of the clutter correction.

Note

The functions work also for a xarray dataset which is linking over the time.

Phase processing

First of all, a phase processing has to be done with the phase_zphi function. Therefore, a binary array (True = rain, False = not-rain) rolling a range-long sum from PHIDP is generated. Afterwards the first occurrence of the maximum (from front and back) of PHIDP binary array is searched for. This gives the indices of the centers of the rolling window. Then min, max, median and mean of these are calculated. At least the true start and end indices of the phase PHIDP are calculated (+- half window length).

[2]:
cphase = wrf.attenuation_corr.phase_zphi(ds_clutter_corr.PHIDP,
                                         rng=1000.0,
                                         start_range=0.0)

display(cphase)
<xarray.Dataset>
Dimensions:       (azimuth: 720, range: 936)
Coordinates:
  * azimuth       (azimuth) float64 0.25 0.75 1.25 1.75 ... 358.8 359.2 359.8
    elevation     (azimuth) float64 0.5 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5 0.5 0.5
    rtime         (azimuth) datetime64[ns] 2022-02-16T13:45:28.941817600 ... ...
  * range         (range) float32 37.5 112.5 187.5 ... 7.009e+04 7.016e+04
    time          datetime64[ns] 2022-02-16T13:45:01
    sweep_mode    <U20 'azimuth_surveillance'
    longitude     float64 13.24
    latitude      float64 53.55
    altitude      float64 38.0
Data variables: (12/15)
    phib          (azimuth, range) float64 nan nan nan nan ... nan nan nan nan
    offset        float64 450.0
    offset_idx    int64 6
    start_range   (azimuth) float32 7.988e+03 7.838e+03 ... 7.988e+03 7.988e+03
    stop_range    (azimuth) float32 2.651e+04 2.621e+04 ... 2.689e+04 2.689e+04
    first_min     (azimuth) float32 32.34 32.34 36.56 ... 36.56 30.94 33.75
    ...            ...
    first_idx     (azimuth) int64 106 104 106 105 105 ... 124 123 123 106 106
    last_min      (azimuth) float32 32.34 30.94 29.53 ... 36.56 35.16 33.75
    last_max      (azimuth) float32 47.81 43.59 43.59 45.0 ... 46.41 45.0 49.22
    last_mean     (azimuth) float32 37.87 35.46 35.76 ... 40.08 39.38 38.97
    last_median   (azimuth) float32 35.86 34.45 36.56 ... 39.38 39.38 38.67
    last_idx      (azimuth) int64 352 348 346 325 325 ... 366 357 357 357 357

ZPHI - Method

In the zphi function the PHIDP is recalcuted and also the KDP. The specific attenuation (AH_ZPHI) is generated and can also be used for the precipitation estimation.

[3]:
ds_zphi=wrf.attenuation_corr.zphi(ds_clutter_corr.DBZH_no_clutter,
                                  cphase=cphase,
                                  alphax=0.28,
                                  bx=0.78)
ds_zphi
[3]:
<xarray.Dataset>
Dimensions:       (azimuth: 720, range: 936)
Coordinates:
  * azimuth       (azimuth) float64 0.25 0.75 1.25 1.75 ... 358.8 359.2 359.8
    elevation     (azimuth) float64 0.5 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5 0.5 0.5
    rtime         (azimuth) datetime64[ns] 2022-02-16T13:45:28.941817600 ... ...
  * range         (range) float32 37.5 112.5 187.5 ... 7.009e+04 7.016e+04
    time          datetime64[ns] 2022-02-16T13:45:01
    sweep_mode    <U20 'azimuth_surveillance'
    longitude     float64 13.24
    latitude      float64 53.55
    altitude      float64 38.0
Data variables:
    PHIDP_RECALC  (azimuth, range) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    AH_ZPHI       (azimuth, range) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    KDP_RECALC    (azimuth, range) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
[4]:
import matplotlib.pyplot as plt

# plot PHIDP Recalc and PHIDP raw
# georefenced dataset for plot
ds_zphi_georef=ds_zphi.pipe(wrl.georef.georeference_dataset)
ds_clutter_corr_georef = ds_clutter_corr.pipe(wrl.georef.georeference_dataset)

# plot
fig = plt.figure(figsize=(14, 5))

# first subplot
ax1 = fig.add_subplot(121)
ds_clutter_corr_georef.PHIDP.plot(x="x", y="y", ax=ax1, vmin=0, vmax=40)
t = plt.title(r'Raw PHIDP')
t.set_y(1.1)

# second subplot
ax2 = fig.add_subplot(122)
ds_zphi_georef.PHIDP_RECALC.plot(x="x", y="y", ax=ax2)
t = plt.title(r'Recalculated PHIDP')
t.set_y(1.1)
../_images/jupyter_notebooks_attenuation_correction_12_0.png

Attenuation correction

Now, the attenuation correction can be done.

[5]:
zhcorr = ds_clutter_corr.DBZH_no_clutter + 0.28 * ds_zphi.PHIDP_RECALC
[6]:
ds_clutter_corr["DBZH_CORR"] = xr.DataArray(zhcorr,
                                            dims=["azimuth", "range"],
                                            attrs=ds_clutter_corr.DBZH.attrs)
[7]:
# smoothing of the dataset
ds_smooth = ds_clutter_corr.pad(azimuth=(3, 3),
                                mode="wrap").rolling(azimuth=7,
                                                     min_periods=3,
                                                     center=True).mean().isel(azimuth=slice(3, -3))

ds_smooth
[7]:
<xarray.Dataset>
Dimensions:          (azimuth: 720, range: 936)
Coordinates: (12/15)
  * azimuth          (azimuth) float64 0.25 0.75 1.25 1.75 ... 358.8 359.2 359.8
    elevation        (azimuth) float64 0.5 0.5 0.5 0.5 0.5 ... 0.5 0.5 0.5 0.5
    rtime            (azimuth) datetime64[ns] 2022-02-16T13:45:28.941817600 ....
  * range            (range) float32 37.5 112.5 187.5 ... 7.009e+04 7.016e+04
    time             datetime64[ns] 2022-02-16T13:45:01
    sweep_mode       <U20 'azimuth_surveillance'
    ...               ...
    x                (azimuth, range) float64 0.1636 0.4909 ... -305.8 -306.1
    y                (azimuth, range) float64 37.5 112.5 ... 7.008e+04 7.015e+04
    z                (azimuth, range) float64 38.3 38.95 39.61 ... 939.2 940.3
    gr               (azimuth, range) float64 37.5 112.5 ... 7.008e+04 7.015e+04
    rays             (azimuth, range) float64 0.25 0.25 0.25 ... 359.8 359.8
    bins             (azimuth, range) float32 37.5 112.5 ... 7.009e+04 7.016e+04
Data variables:
    RATE             (azimuth, range) float64 nan nan nan nan ... nan nan nan
    DBZH             (azimuth, range) float64 nan nan nan nan ... nan nan nan
    VRAD             (azimuth, range) float64 nan nan nan nan ... nan nan nan
    ZDR              (azimuth, range) float64 nan nan nan nan ... nan nan nan
    KDP              (azimuth, range) float64 nan nan nan nan ... nan nan nan
    PHIDP            (azimuth, range) float64 nan nan nan nan ... nan nan nan
    RHOHV            (azimuth, range) float64 nan nan nan nan ... nan nan nan
    WRAD             (azimuth, range) float64 nan nan nan nan ... nan nan nan
    QC_INFO          (azimuth, range) float64 nan nan nan nan ... nan nan nan
    FUZZ             (azimuth, range) float64 nan nan nan nan ... nan nan nan
    DBZH_no_clutter  (azimuth, range) float64 nan nan nan nan ... nan nan nan
    DBZH_CORR        (azimuth, range) float64 nan nan nan nan ... nan nan nan
Attributes:
    fixed_angle:  0.5
[8]:
# plot PHIDP Recalc and PHIDP raw

# georefenced dataset for plot
ds_smooth_georef=ds_smooth.pipe(wrl.georef.georeference_dataset)

# plot
fig = plt.figure(figsize=(14, 5))

# first subplot
ax1 = fig.add_subplot(121)
ds_smooth_georef.DBZH.plot(x="x", y="y", ax=ax1, cmap = "jet", vmin=0, vmax=40)
t = plt.title(r'Raw reflectivity')
t.set_y(1.1)

# second subplot
ax2 = fig.add_subplot(122)
ds_smooth_georef.DBZH_CORR.plot(x="x", y="y", ax=ax2, cmap="jet", vmin=0, vmax=40)
t = plt.title(r'Attenuation corrected reflecivity')
t.set_y(1.1)
../_images/jupyter_notebooks_attenuation_correction_18_0.png

Library Reference

Seealso

Get more information about the attenuation correction module in the library reference section.