Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Combined usage of packages

import matplotlib.pyplot as plt

import poligrain as plg
import pypwsqc
import mergeplg
import pycomlink as pycml

# Suppress numpy and pypwsqc runtime warnings
import warnings
warnings.filterwarnings('ignore', category=RuntimeWarning, module='numpy')
warnings.filterwarnings('ignore', category=RuntimeWarning, module='pypwsqc')

Load open OpenMRG example data and project coordinates

We load 8 days of the OpenMRG dataset and do some quick visualization of the data using poligrain.

(
    ds_rad,
    ds_cmls,
    ds_gauges_municp,
    ds_gauge_smhi,
) = plg.example_data.load_openmrg(data_dir="example_data", subset="8d")
File already exists at example_data/openmrg_cml_8d.nc
Not downloading!
File already exists at example_data/openmrg_rad_8d.nc
Not downloading!
File already exists at example_data/openmrg_municp_gauge_8d.nc
Not downloading!
File already exists at example_data/openmrg_smhi_gauge_8d.nc
Not downloading!

# Project data to coordinate system
# UTM32N: https://epsg.io/32632
ref_str = "EPSG:32632"

# gauges
ds_gauges_municp.coords["x"], ds_gauges_municp.coords["y"] = plg.spatial.project_point_coordinates(
    ds_gauges_municp.lon, ds_gauges_municp.lat, ref_str
)

# CMLs
ds_cmls.coords["site_0_x"], ds_cmls.coords["site_0_y"] = plg.spatial.project_point_coordinates(
    ds_cmls.site_0_lon, ds_cmls.site_0_lat, ref_str
)
ds_cmls.coords["site_1_x"], ds_cmls.coords["site_1_y"] = plg.spatial.project_point_coordinates(
    ds_cmls.site_1_lon, ds_cmls.site_1_lat, ref_str
)

# Radar
ds_rad.coords["x_grid"], ds_rad.coords["y_grid"] = plg.spatial.project_point_coordinates(
    ds_rad.lon, ds_rad.lat, ref_str
)
plg.plot_map.plot_plg(
    da_grid=ds_rad.R.resample(time='1h').mean().sum(dim="time"),
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim='time'),
    da_cmls=ds_cmls,
    vmin=0,
    vmax=100,
    cmap='YlGnBu',
    kwargs_cmls_plot={"line_color": "k", "line_width": 1},
    kwargs_gauges_plot={"edge_color": "r", "s": 30},
    colorbar_label="rainfall sum (mm)",
)
<Axes: xlabel='lon', ylabel='lat'>
<Figure size 640x480 with 2 Axes>

CML processing

Here we do a standard CML processing with, what we believe, are reasonable methods and parameters. Note that, we use the radar as wet-dry indicator and we use poligrain to extract the path-averaged radar data for each CML which is needed for that.

ds_cmls['tl'] = ds_cmls.tsl - ds_cmls.rsl

# get radar along CMLs
grid_at_lines = plg.spatial.GridAtLines(
    da_gridded_data=ds_rad.R,
    ds_line_data=ds_cmls,
)
da_radar_along_cmls = grid_at_lines(ds_rad.R)

# wet-dry classification based on 
# smoothed radar rain rates along each CML
ds_cmls['wet'] = (
    (da_radar_along_cmls > 0.1)
    .rolling(time=6*3, center=True)
    .max()
)
ds_cmls['wet'] = ds_cmls.wet.ffill(dim='time')

#ds_cmls['wet'] = ds_cmls.tl.rolling(time=10*60).std() > 0.3

# calculate baseline
ds_cmls["baseline"] = pycml.processing.baseline.baseline_constant(
    trsl=ds_cmls.tl,
    wet=ds_cmls.wet,
    n_average_last_dry=5,
)

# calculate wet antenna effect
ds_cmls["A_obs"] = ds_cmls.tl - ds_cmls.baseline
ds_cmls["A_obs"] = ds_cmls.A_obs.where(ds_cmls.A_obs >= 0, 0)
ds_cmls["waa"] = pycml.processing.wet_antenna.waa_pastorek_2021_from_A_obs(
    A_obs=ds_cmls.A_obs,
    f_Hz=ds_cmls.frequency * 1e6,
    pol=ds_cmls.polarization.data,
    L_km=ds_cmls.length / 1000,
    A_max=6,
    zeta=0.7,
    d=0.15,
)

# calculate attenuation caused by rain and remove negative attenuation
ds_cmls["A"] = ds_cmls.tl - ds_cmls.baseline - ds_cmls.waa
ds_cmls["A"].data[ds_cmls.A < 0] = 0
# derive rain rate via the k-R relation
ds_cmls["R"] = pycml.processing.k_R_relation.calc_R_from_A(
    A=ds_cmls.A,
    L_km=ds_cmls.length.astype(float) / 1000,
    f_GHz=ds_cmls.frequency / 1000,
    pol=ds_cmls.polarization,
)
ds_cmls["R"].data[ds_cmls.R < 0.1] = 0

# Calcualte rainfall sums for plotting
rad_sum = ds_rad.R.resample(time='1h').mean().sum(dim='time')
cml_sum = ds_cmls.R.isel(sublink_id=0).resample(time='1h').mean().sum(dim='time', skipna=False)
rad_along_cml_sum = da_radar_along_cmls.resample(time='1h').mean().sum(dim='time')
# Plot data on maps
fig, axs = plt.subplots(1, 3, figsize=(18, 5), sharex=True, sharey=True)

plg.plot_map.plot_plg(
    da_grid=rad_sum,
    da_cmls=cml_sum, 
    vmin=0, 
    vmax=100,
    ax=axs[0],
    kwargs_cmls_plot={"line_width": 1},
)
plg.plot_map.plot_plg(
    da_grid=rad_sum,
    da_cmls=rad_along_cml_sum,
    kwargs_cmls_plot={"line_width": 1},
    vmin=0, 
    vmax=100,
    ax=axs[1],
)
plg.plot_map.plot_plg(
    #da_grid=rad_sum,
    da_cmls=cml_sum - rad_along_cml_sum,
    kwargs_cmls_plot={"line_width": 3, "vmin": -30, "vmax": 30, "cmap": 'RdBu'},
    vmin=0, 
    vmax=100,
    #cmap="Grays",
    ax=axs[2],
)
<Axes: >
<Figure size 1800x500 with 6 Axes>
# plot for paper
import matplotlib as mpl
R_cmap = mpl.colormaps.get_cmap('YlGnBu')
R_cmap.set_bad(color='gray')
diff_cmap = mpl.colormaps.get_cmap('RdBu')
diff_cmap.set_bad(color='gray')

# Plot data on maps
fig, axs = plt.subplots(2, 1, figsize=(5, 8), sharex=True, sharey=True)

plg.plot_map.plot_plg(
    da_grid=rad_sum,
    da_cmls=cml_sum, 
    vmin=0, 
    vmax=100,
    ax=axs[0],
    cmap=R_cmap,
    kwargs_cmls_plot={"line_width": 1},
    colorbar_label="rainfall sum (mm)",

)

plg.plot_map.plot_plg(
    #da_grid=rad_sum,
    da_cmls=cml_sum - rad_along_cml_sum,
    kwargs_cmls_plot={"line_width": 1, "vmin": -30, "vmax": 30, "cmap": diff_cmap},
    vmin=0, 
    vmax=100,
    ax=axs[1],
    colorbar_label="CML-Radar Difference (mm)",
);
<Figure size 500x800 with 4 Axes>

Apply the indicator correlation filter to detect outliers

Note that the indicator correlation filter was developed for quality control of PWS data, but we will use it here for the processed CML rainfall data to examplify the usage of tools.

import pypwsqc.indicator_correlation as ic
# Calculate CML midpoint which is then used to treat CMLs as virtual points with coordinate x and y
ds_cmls.coords["x"] = (ds_cmls.site_0_x + ds_cmls.site_1_x) / 2
ds_cmls.coords["y"] = (ds_cmls.site_0_y + ds_cmls.site_1_y) / 2
# Aggregate to 15 minutes
da_gauges = 4 * ds_gauges_municp.rainfall_amount.resample(time='15min').sum()
da_cmls = (
    ds_cmls.R.isel(sublink_id=0, drop=True)
    .resample(time='15min').mean()
    .rename({'cml_id': 'id'})
) 

# Calculate indicator correlation and distance matrix
dist_mtx_ref, indcorr_mtx_ref = ic.indicator_distance_matrix(
    da_gauges,
    da_gauges,
    max_distance=30e3,
    prob=0.95,
    min_valid_overlap=6 * 24,
)
  0%|          | 0/10 [00:00<?, ?it/s]
 70%|███████   | 7/10 [00:00<00:00, 62.27it/s]
100%|██████████| 10/10 [00:00<00:00, 59.73it/s]


plt.scatter(dist_mtx_ref, indcorr_mtx_ref, color="red", s=10, label="Ref")
plt.ylim(0, 1)
plt.xlim(0, 30e3)
plt.ylabel("Indicator Correlation [-]")
plt.xlabel("Distance [m]")
plt.title("Indicator correlation for reference data vs reference data")
plt.legend();
<Figure size 640x480 with 1 Axes>
# Distance and indicator correlations of reference data
dist_mtx, indcorr_mtx = ic.indicator_distance_matrix(
    da_gauges,
    da_cmls,
    max_distance=30e3,
    prob=0.95,
    min_valid_overlap=6 * 24,
)
  0%|          | 0/10 [00:00<?, ?it/s]
 10%|█         | 1/10 [00:00<00:06,  1.49it/s]
 20%|██        | 2/10 [00:01<00:05,  1.44it/s]
 30%|███       | 3/10 [00:02<00:04,  1.47it/s]
 40%|████      | 4/10 [00:02<00:03,  1.63it/s]
 50%|█████     | 5/10 [00:03<00:02,  1.70it/s]
 60%|██████    | 6/10 [00:03<00:02,  1.84it/s]
 70%|███████   | 7/10 [00:03<00:01,  2.00it/s]
 80%|████████  | 8/10 [00:04<00:01,  1.85it/s]
 90%|█████████ | 9/10 [00:05<00:00,  1.71it/s]
100%|██████████| 10/10 [00:05<00:00,  1.68it/s]
100%|██████████| 10/10 [00:05<00:00,  1.70it/s]


plt.scatter(dist_mtx, indcorr_mtx, color="b", alpha=0.2, s=10, label="CML")
plt.scatter(dist_mtx_ref, indcorr_mtx_ref, color="red", s=10, label="Ref")
plt.ylim(0, 1.05)
plt.xlim(0, 30e3)
plt.ylabel("Indicator Correlation [-]")
plt.xlabel("Distance [m]")
plt.title("Indicator correlation for reference data vs reference data and CML vs reference data")
plt.legend();
<Figure size 640x480 with 1 Axes>
# Filter out CMLs with low indicator correlation
indcorr_results = ic.ic_filter(
    indicator_correlation_matrix_ref=indcorr_mtx_ref,
    distance_correlation_matrix_ref=dist_mtx_ref,
    indicator_correlation_matrix=indcorr_mtx,
    distance_matrix=dist_mtx,
    max_distance=40000,
    bin_size=1000,
    quantile_bin_ref=0.5,
    quantile_bin_pws=0.5,
    threshold=0.0,
);

indcorr_results
Loading...
print(
    str(indcorr_results.indcorr_good.data.sum())
    + " of "
    + str(len(indcorr_results.indcorr_good))
    + " CMLs were accepted"
)
256 of 364 CMLs were accepted

cml_ids_bad_indcorr = indcorr_results.id_neighbor.where(
    ~indcorr_results.indcorr_good & (indcorr_results.dist.max(dim='id') < 16000), drop=True,
)
cml_ids_bad_indcorr
Loading...
# Plot data of one CML with low indicator correlation
i = 5
cml_id = cml_ids_bad_indcorr[i].data

cml_id = 10058
gauge_id = 'Lbom'

fig, axs = plt.subplots(3, 1, figsize=(6, 10))
plt.sca(axs[0])
indcorr_results.plot.scatter(x="dist", y="indcorr", color="grey", alpha=0.5, s=10, label='all CML vs gauges')
plt.scatter(dist_mtx_ref, indcorr_mtx_ref, color="r", alpha=0.5, s=10, label="gauges vs gauges")
indcorr_results.sel(id_neighbor=cml_id).plot.scatter(
    x="dist", y="indcorr", color="lime", s=15, label=f'CML ID {cml_id:.0f} vs. gauges'
)
plt.legend()
plt.title('')
plt.sca(axs[1])
da_cmls.sel(id=cml_id).plot(label=f'CML {cml_id}')
da_gauges.sel(id=gauge_id).plot(label=f'Gauge {gauge_id}')
plt.legend()

plt.sca(axs[2])
ds_cmls.sel(cml_id=cml_id).isel(sublink_id=0).tl.plot(label=f'CML {cml_id}')
plt.legend()

axs[1].set_title('')
axs[2].set_title('')
axs[1].set_xlabel('')
axs[1].set_xticklabels([]);
<Figure size 600x1000 with 3 Axes>
# Plot CMLs with low indicator correlation
for cml_id in cml_ids_bad_indcorr.values:
    fig, axs = plt.subplots(2, 1)
    ds_cmls.sel(cml_id=cml_id).isel(sublink_id=0).tl.plot(ax=axs[0])
    ds_cmls.sel(cml_id=cml_id).isel(sublink_id=0).R.plot(ax=axs[1])
    plt.show()
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>
<Figure size 640x480 with 2 Axes>

Adjust radar with CMLs using mergeplg

import mergeplg as mrg

# Set up object to get radar values at rain gauges for validation
grid_at_points = plg.spatial.GridAtPoints(
    da_gridded_data=ds_rad.R,
    da_point_data=ds_gauges_municp.rainfall_amount,
    nnear=1,
)

# Set up the merging object which internally calculates 
# and stores the CML-radar intersection weights
merge_diff_idw = mrg.merge.MergeDifferenceIDW(
    ds_rad=ds_rad.R, 
    ds_cmls=ds_cmls, 
    method='additive',
    nnear=12,
)

# Aggregate to 1h
R_radar = ds_rad.R.resample(time='1h').mean()
R_cmls = ds_cmls.isel(sublink_id=0).R.resample(time='1h').mean()
R_gauges = ds_gauges_municp.rainfall_amount.resample(time='1h').sum()

# Do adjustment for one selected timestamp
t = '2015-07-26 03:00:00'
R_grid_idw = merge_diff_idw(
    da_rad=R_radar.sel(time=t),
    da_cmls=R_cmls.sel(time=t),
)

# Get diff of radar and gauges at gauge location
R_rad_at_gauge = grid_at_points(
    da_gridded_data=R_radar.sel(time=t),
    da_point_data=R_gauges.sel(time=t),
)
R_adjusted_rad_at_gauge = grid_at_points(
    da_gridded_data=R_grid_idw,
    da_point_data=R_gauges.sel(time=t),
)

# Plotting the result
fig, axs = plt.subplots(3, 1, figsize=(5, 13))
vmin, vmax = 0, 15
vmin_diff, vmax_diff = -7, 7

plg.plot_map.plot_plg(
    da_grid=R_radar.sel(time=t),
    da_gauges=R_gauges.sel(time=t),
    vmin=vmin, 
    vmax=vmax,
    ax=axs[0],
    kwargs_cmls_plot={"line_width": 1},
    colorbar_label="1h rainfall sum (mm)",
)
plg.plot_map.plot_plg(

    da_grid=R_grid_idw,
    da_cmls=R_cmls.sel(time=t),
    da_gauges=R_gauges.sel(time=t),
    kwargs_cmls_plot={"line_width": 1},
    vmin=vmin, 
    vmax=vmax,
    ax=axs[1],
    colorbar_label="1h rainfall sum (mm)",

)
axs[2].scatter(
    R_gauges.sel(time=t).data,
    R_rad_at_gauge.data,
    c='k',
    label='unadjusted radar'
)
axs[2].scatter(
    R_gauges.sel(time=t).data,
    R_adjusted_rad_at_gauge.data,
    c='r',
    label='CML-adjusted radar'
)
axs[2].legend()
axs[2].plot([-2, 22], [-2, 22], 'k:', linewidth=0.5)

axs[0].set_title('Unadjusted radar')
axs[1].set_title('CML-adjusted radar')

axs[0].set_xlim(11.655, 12.255)
axs[0].set_ylim(57.555, 57.845)
axs[1].set_xlim(11.655, 12.255)
axs[1].set_ylim(57.555, 57.845)
axs[0].set_xlabel('')
axs[0].set_ylabel('')
axs[1].set_xlabel('')
axs[1].set_ylabel('')

axs[2].set_xlim(-1, 21)
axs[2].set_ylim(-1, 21)
axs[2].set_xlabel('Gauge rainfall [mm]')
axs[2].set_ylabel('Radar rainfall [mm]')
axs[2].set_xticks(range(0, 22, 5))
axs[2].set_yticks(range(0, 22, 5));


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[20], line 36
     32 R_rad_at_gauge = grid_at_points(
     33     da_gridded_data=R_radar.sel(time=t),
     34     da_point_data=R_gauges.sel(time=t),
     35 )
---> 36 R_adjusted_rad_at_gauge = grid_at_points(
     37     da_gridded_data=R_grid_idw,
     38     da_point_data=R_gauges.sel(time=t),
     39 )

File ~/work/OpenSenseAction.github.io/OpenSenseAction.github.io/.venv/lib/python3.11/site-packages/poligrain/spatial.py:452, in GridAtPoints.__call__(self, da_gridded_data, da_point_data)
    449 da_point_data_t = da_point_data.sel(time=t)
    451 # copy-paste code from wradlib.adjust.RawAtObs.__call__
--> 452 raw = da_gridded_data_t.to_numpy().flatten()
    453 obs = da_point_data_t.to_numpy()
    454 # get the values of the raw neighbours of obs

File ~/work/OpenSenseAction.github.io/OpenSenseAction.github.io/.venv/lib/python3.11/site-packages/xarray/core/common.py:286, in AttrAccessMixin.__getattr__(self, name)
    284         with suppress(KeyError):
    285             return source[name]
--> 286 raise AttributeError(
    287     f"{type(self).__name__!r} object has no attribute {name!r}"
    288 )

AttributeError: 'Dataset' object has no attribute 'to_numpy'
fig, axs = plt.subplots(2, 2, figsize=(12, 10), sharex=True, sharey=True)

plg.plot_map.plot_plg(
    da_grid=R_radar.sel(time=t),
    da_gauges=R_gauges.sel(time=t),
    vmin=vmin, 
    vmax=vmax,
    ax=axs[0, 0],
    kwargs_cmls_plot={"line_width": 1},
    colorbar_label="1h rainfall sum (mm)",

)
plg.plot_map.plot_plg(

    da_grid=R_grid_idw,
    da_cmls=R_cmls.sel(time=t),
    da_gauges=R_gauges.sel(time=t),
    kwargs_cmls_plot={"line_width": 1},
    vmin=vmin, 
    vmax=vmax,
    ax=axs[1, 0],
    colorbar_label="1h rainfall sum (mm)",
)

plg.plot_map.plot_plg(
    #da_grid=R_radar.sel(time=t),
    da_gauges=R_rad_at_gauge - R_gauges.sel(time=t),
    vmin=vmin,
    vmax=vmax,
    cmap='Grays',
    ax=axs[0, 1],
    colorbar_label="Radar-Gauge Difference (mm)",
    kwargs_gauges_plot={'cmap': 'RdBu', 'vmin': vmin_diff, 'vmax': vmax_diff},
)

plg.plot_map.plot_plg(
    #da_grid=R_grid_idw,
    da_gauges=R_adjusted_rad_at_gauge - R_gauges.sel(time=t),
    vmin=vmin, 
    vmax=vmax,
    cmap='Grays',
    ax=axs[1, 1],
    colorbar_label="Adjusted Radar-Gauge Difference (mm)",
    kwargs_gauges_plot={'cmap': 'RdBu', 'vmin': vmin_diff, 'vmax': vmax_diff},
);

axs[0, 0].set_xlim(11.655, 12.255)
axs[0, 0].set_ylim(57.555, 57.845)

axs[0, 0].set_title('Unadjusted radar')
axs[1, 0].set_title('CML-adjusted radar')
axs[0, 1].set_title('Radar-Gauge Difference')
axs[1, 1].set_title('Adjusted Radar-Gauge Difference')
OpenSense Software
mergeplg quick tour
OpenSense Software
SML processing example