## Background#

There are various ways to correct specific errors and artifacts in radar-based quantitative precipitation estimates (radar QPE). Alternatively, you might want to correct your radar QPE regardless of the error source - by using ground truth, or, more specifically, rain gauge observations. Basically, you define the error of your radar QPE at a rain gauge location by the discrepancy between rain gauge observation (considered as “the truth”) and radar QPE at that very location. Whether you consider this “discrepancy” as an additive or multiplicative error is somehow arbitrary - typically, it’s a mix of both. If you quantify this error at various locations (i.e. rain gauges), you can go ahead and construct correction fields for your radar QPE. You might compute a single correction factor for your entire radar domain (which would e.g. make sense in case of hardware miscalibration), or you might want to compute a spatially variable correction field. This typically implies to interpolate the error in space.

$$\omega radlib$$ provides different error models and different spatial interpolation methods to address the adjustment problem. For details, please refer to $$\omega radlib's$$ library reference.

[1]:
import numpy as np
import matplotlib.pyplot as plt

try:
get_ipython().run_line_magic("matplotlib inline")
except:
plt.ion()

## Example for the 1-dimensional case#

Looking at the 1-D (instead of 2-D) case is more illustrative.

### Create synthetic data#

First, we create synthetic data: - true rainfall, - point observations of the truth, - radar observations of the truth.

The latter is disturbed by some kind of error, e.g. a combination between systemtic and random error.

[2]:
obs_coords = np.array([5, 10, 15, 20, 30, 45, 65, 70, 77, 90])

# true rainfall
np.random.seed(1319622840)
truth = np.abs(1.5 + np.sin(0.075 * radar_coords)) + np.random.uniform(
)

errormult = 0.75 + 0.015 * radar_coords

# gage observations are assumed to be perfect
obs = truth[obs_coords]

# add a missing value to observations (just for testing)
obs[1] = np.nan

• multiplicative error, spatially variable (AdjustMultiply)

• mixed error, spatially variable (AdjustMixed)

• multiplicative error, spatially uniform (AdjustMFB)

[3]:
# number of neighbours to be used
nnear_raws = 3

)

)
weights = 1.0 / self.dists**self.p

[4]:
# Enlarge all label fonts
font = {"size": 15}
plt.rc("font", **font)

plt.figure(figsize=(10, 5))
plt.plot(
"k",
linewidth=2.0,
linestyle="dashed",
)
plt.plot(
truth,
"k-",
linewidth=2.0,
label="True rainfall",
)
plt.plot(
obs_coords,
obs,
"o",
markersize=10.0,
markerfacecolor="grey",
label="Gage observation",
)
plt.plot(
)
plt.plot(
)
plt.plot(
"-",
color="blue",
)
plt.xlabel("Distance (km)")
plt.ylabel("Rainfall intensity (mm/h)")
leg = plt.legend(prop={"size": 10})

### Verification#

We use the verify module to compare the errors of different adjustment approaches.

Here, we compare the adjustment to the “truth”. In practice, we would carry out a cross validation.

[5]:
# Verification for this example
[6]:
# Helper function for scatter plot
def scatterplot(x, y, title=""):
"""Quick and dirty helper function to produce scatter plots"""
plt.scatter(x, y)
plt.plot([0, 1.2 * maxval], [0, 1.2 * maxval], "-", color="grey")
plt.xlabel("True rainfall (mm)")
plt.ylabel("Estimated rainfall (mm)")
plt.xlim(0, maxval + 0.1 * maxval)
plt.ylim(0, maxval + 0.1 * maxval)
plt.title(title)
[7]:
# Verification reports
maxval = 4.0
# Enlarge all label fonts
font = {"size": 10}
plt.rc("font", **font)
fig = plt.figure(figsize=(14, 8))
ax.text(0.2, maxval, "Nash=%.1f" % rawerror.nash(), fontsize=12)
ax.text(0.2, maxval, "Nash=%.1f" % adderror.nash(), fontsize=12)
ax.text(0.2, maxval, "Nash=%.1f" % multerror.nash(), fontsize=12)
ax.text(0.2, maxval, "Nash=%.1f" % mixerror.nash(), fontsize=12)
scatterplot(mfberror.obs, mfberror.est, title="Mean Field Bias adjustment")
ax.text(0.2, maxval, "Nash=%.1f" % mfberror.nash(), fontsize=12)
plt.tight_layout()

## Example for the 2-dimensional case#

For the 2-D case, we follow the same approach as before:

• create synthetic data: truth, rain gauge observations, radar-based rainfall estimates

• verification

The way these synthetic data are created is totally arbitrary - it’s just to show how the methods are applied.

### Create 2-D synthetic data#

[8]:
# grid axes
xgrid = np.arange(0, 10)
ygrid = np.arange(20, 30)

# number of observations
num_obs = 10

# create grid
gridshape = len(xgrid), len(ygrid)
grid_coords = wrl.util.gridaspoints(ygrid, xgrid)

# Synthetic true rainfall
truth = np.abs(10.0 * np.sin(0.1 * grid_coords).sum(axis=1))

# Creating radar data by perturbing truth with multiplicative and
# YOU CAN EXPERIMENT WITH THE ERROR STRUCTURE
np.random.seed(1319622840)
radar = 0.6 * truth + 1.0 * np.random.uniform(low=-1.0, high=1, size=len(truth))

# indices for creating obs from raw (random placement of gauges)
obs_ix = np.random.uniform(low=0, high=len(grid_coords), size=num_obs).astype("i4")

# creating obs_coordinates
obs_coords = grid_coords[obs_ix]

# creating gauge observations from truth
obs = truth[obs_ix]

[9]:

# Multiplicative Error Model

[10]:
# Helper functions for grid plots
def gridplot(data, title):
"""Quick and dirty helper function to produce a grid plot"""
xplot = np.append(xgrid, xgrid[-1] + 1.0) - 0.5
yplot = np.append(ygrid, ygrid[-1] + 1.0) - 0.5
grd = ax.pcolormesh(xplot, yplot, data.reshape(gridshape), vmin=0, vmax=maxval)
ax.scatter(
obs_coords[:, 0],
obs_coords[:, 1],
c=obs.ravel(),
marker="s",
s=50,
vmin=0,
vmax=maxval,
)
# plt.colorbar(grd, shrink=0.5)
plt.title(title)
[11]:
# Maximum value (used for normalisation of colorscales)

# open figure
fig = plt.figure(figsize=(10, 6))

# True rainfall
gridplot(truth, "True rainfall")

plt.tight_layout()
[12]:
# Open figure
fig = plt.figure(figsize=(6, 6))

# Scatter plot radar vs. observations