Example of the rain event detection and rainfall retrieval using the nearby appraoch from Overeem et al. (2016)

[1]:

%matplotlib inline

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

import pycomlink.processing.wet_dry.nearby_wetdry as nearby_wetdry
import pycomlink.processing.nearby_rain_retrival as nearby_rain

Calculate the distance matrix

Calculating distances between all cml endpoints and plotting the neighbors used for wet-dry classification depending on the distance

[ ]:

ds_dist = nearby_wetdry.calc_distance_between_cml_endpoints(
    cml_ids=cmls.cml_id.values,
    site_a_latitude=cmls.site_a_latitude,
    site_a_longitude=cmls.site_a_longitude,
    site_b_latitude=cmls.site_b_latitude,
    site_b_longitude=cmls.site_b_longitude,
)

500it [00:01, 350.47it/s]

[6]:

r = 15  # radius in km
ds_dist["within_r"] = (
    (ds_dist.a_to_all_a < r)
    & (ds_dist.a_to_all_b < r)
    & (ds_dist.b_to_all_a < r)
    & (ds_dist.b_to_all_b < r)
)

[7]:

ds_dist.within_r.sum(dim="cml_id2").plot.hist(
    bins=int(ds_dist.within_r.sum(dim="cml_id2").max())
)
plt.vlines(4, ymin=0, ymax=50, color="red")
plt.annotate(
    text="sufficient dense\nCML network for\nr>3",
    xy=(4, 45),
    xytext=(10, 43),
    arrowprops=dict(arrowstyle="<-"),
)
plt.xlabel("CMLs within radius r")
plt.ylabel("count");

../_images/notebooks_Nearby_link_approach_processing_example_13_0.png

Classification of wet and dry time steps

using default paramters from Overeem et al. (2016)

[8]:

wet, F = nearby_wetdry.nearby_wetdry(
    pmin=pmin,
    ds_dist=ds_dist,
    radius=15,
    thresh_median_P=-1.4,
    thresh_median_PL=-0.7,
    min_links=3,
    interval=15,
    timeperiod=24,
)

100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 500/500 [00:10<00:00, 48.77it/s]

Plotting instantaneous data, min-max data and the rain event detection from the nearby wetdry approach

[9]:

t_start, t_end = "2018-05-13", "2018-05-15"
for cmlid in ["5", "66", "78", "100"]:
    (cmls.sel(cml_id=cmlid).tsl - cmls.sel(cml_id=cmlid).rsl).isel(channel_id=0).sel(
        time=slice(t_start, t_end)
    ).plot(figsize=(8, 3), label="TL: total loss")
    pmin.sel(cml_id=cmlid, time=slice(t_start, t_end)).isel(channel_id=0).plot(
        label="pmin: 15-minute minimum of rsl - tsl"
    )
    pmax.sel(cml_id=cmlid, time=slice(t_start, t_end)).isel(channel_id=0).plot(
        label="pmax: 15-minute maximum of rsl - tsl"
    )
    (
        (wet.isel(channel_id=0).sel(cml_id=cmlid, time=slice(t_start, t_end)) * 200)
        - 100
    ).plot(label="classified wet", alpha=0.5)

    plt.legend()

../_images/notebooks_Nearby_link_approach_processing_example_17_0.png

../_images/notebooks_Nearby_link_approach_processing_example_17_1.png

../_images/notebooks_Nearby_link_approach_processing_example_17_2.png

../_images/notebooks_Nearby_link_approach_processing_example_17_3.png

Baseline estimation (pref)

The baseline is set to the rolling median of the average of pmin and pmax during the last n_average_dry time steps. Default is 24 hours (n_average_dry=96 because 15 minute temporal resolution) and a baseline will be calculated when at least one time step is available (min_periods=1) If pmax is not available e.g. because the min-max data is derived from instanteanous sampled CML data and has the same temporal resolution as the instanteanous CML data, substitute pmax with pmin so pmin and pmax are identical.

[10]:

pref = nearby_rain.nearby_determine_reference_level(pmin, pmax, wet, n_average_dry=96, min_periods=1)

Correction if min-max data

Correcting pmin and pmax so that no rainfall estimation is carried out during dry time steps. All time steps of pmin which are not classified wet and pmin is smaller than pref are set to pref. Similarly, all time steps of pmax where either the corrected pmin (p_c_min) is not smaller than pref or pmax is not smaller than pref are set to pref. This ensures that only wet time steps are used for rainfall estimation an and that pmax is not above pref which would lead to an overestimation of rainfall.

[11]:

p_c_min, p_c_max = nearby_rain.nearby_correct_received_signals(pmin, pmax, wet, pref)

Plot pref (baseline) and corrected pmin and pmax

[12]:

t_start, t_end = "2018-05-13", "2018-05-15"
for cmlid in ["5", "66", "78", "100"]:
    pref.sel(cml_id=cmlid, time=slice(t_start, t_end)).isel(channel_id=0).plot(
        label="pref: baseline level",figsize=(8, 3),color='black'
    )
    p_c_min.sel(cml_id=cmlid, time=slice(t_start, t_end)).isel(channel_id=0).plot(
        label="p_c_min: corrected pmin"
    )
    p_c_max.sel(cml_id=cmlid, time=slice(t_start, t_end)).isel(channel_id=0).plot(
        label="p_c_max: corrected pmax"
    )
    plt.legend()

../_images/notebooks_Nearby_link_approach_processing_example_23_0.png

../_images/notebooks_Nearby_link_approach_processing_example_23_1.png

../_images/notebooks_Nearby_link_approach_processing_example_23_2.png

../_images/notebooks_Nearby_link_approach_processing_example_23_3.png

Calculate rain rates from attenuation data

Also, correcting for wet antenna attenuation and setting the alpha value which defines how close to the minimum attenuation of each intervall the rain rate should be set.

[13]:

R = nearby_rain.nearby_rainfall_retrival(
    pref,
    p_c_min,
    p_c_max,
    F,
    length=pmin.length,
    f_GHz=pmin.frequency / 1e9,
    pol=pmin.polarization,
    waa_max=2.3,
    alpha=0.33,
    F_value_threshold =-32.5,
)

Compare derived rain rates with reference data

As reference, 5-minute path-averaged rainfall sums along the CMLs paths from RADKLIM-YW are provided. This data is resampled to 15 minute rainfall intensities to match the CML derived rainfall intensities. Here CML timeseries are compared individually against reference timeseries. (for simplicity only channel 1 is evaluated here)

[14]:

path_ref = xr.open_dataset(data_path + "/example_path_averaged_reference_data.nc")

[15]:

for i in ["5", "66", "78", "100"]:
    # Plot reference rainfall amount (converted to 15-minute rainfall rate)
    (path_ref.sel(cml_id=i).rainfall_amount.resample(time="15min").sum() * 4).plot(
        label="RADKLIM_YW", color="C3", figsize=(12, 3)
    )
    # Plot 15-minute mean rainfall rates from CMLs
    (R.sel(cml_id=i, channel_id="channel_1")).plot(
        x="time", label="CML_nearby", color="C0"
    )

    plt.xlim(np.datetime64("2018-05-13"), np.datetime64("2018-05-15"))
    plt.ylabel("5-min rainfall rate (mm/h)")
    plt.legend();