How to use wradlib’s ipol module for interpolation tasks?

import wradlib as wrl
import numpy as np
import matplotlib.pyplot as plt
import datetime as dt
import warnings

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

1-dimensional example

Includes Nearest Neighbours, Inverse Distance Weighting, and Ordinary Kriging.

# Synthetic observations
xsrc = np.arange(10)[:, None]
vals = np.sin(xsrc).ravel()

# Define target coordinates
xtrg = np.linspace(0, 20, 100)[:, None]

# Set up interpolation objects
#   IDW
idw = wrl.ipol.Idw(xsrc, xtrg)
#   Nearest Neighbours
nn = wrl.ipol.Nearest(xsrc, xtrg)
#   Linear
ok = wrl.ipol.OrdinaryKriging(xsrc, xtrg)

# Plot results
plt.figure(figsize=(10, 5))
plt.plot(xsrc.ravel(), vals, "bo", label="Observation")
plt.plot(xtrg.ravel(), idw(vals), "r-", label="IDW interpolation")
plt.plot(xtrg.ravel(), nn(vals), "k-", label="Nearest Neighbour interpolation")
plt.plot(xtrg.ravel(), ok(vals), "g-", label="Ordinary Kriging")
plt.xlabel("Distance", fontsize="large")
plt.ylabel("Value", fontsize="large")
plt.legend(loc="lower right")

<matplotlib.legend.Legend at 0x7f10fa76ffe0>

2-dimensional example

Includes Nearest Neighbours, Inverse Distance Weighting, Linear Interpolation, and Ordinary Kriging.

# Synthetic observations and source coordinates
src = np.vstack((np.array([4, 7, 3, 15]), np.array([8, 18, 17, 3]))).transpose()
np.random.seed(1319622840)
vals = np.random.uniform(size=len(src))

# Target coordinates
xtrg = np.linspace(0, 20, 40)
ytrg = np.linspace(0, 20, 40)
trg = np.meshgrid(xtrg, ytrg)
trg = np.vstack((trg[0].ravel(), trg[1].ravel())).T

# Interpolation objects
idw = wrl.ipol.Idw(src, trg)
nn = wrl.ipol.Nearest(src, trg)
linear = wrl.ipol.Linear(src, trg)
ok = wrl.ipol.OrdinaryKriging(src, trg)


# Subplot layout
def gridplot(interpolated, title=""):
    pm = ax.pcolormesh(xtrg, ytrg, interpolated.reshape((len(xtrg), len(ytrg))))
    plt.axis("tight")
    ax.scatter(src[:, 0], src[:, 1], facecolor="None", s=50, marker="s")
    plt.title(title)
    plt.xlabel("x coordinate")
    plt.ylabel("y coordinate")


# Plot results
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(221, aspect="equal")
gridplot(idw(vals), "IDW")
ax = fig.add_subplot(222, aspect="equal")
gridplot(nn(vals), "Nearest Neighbours")
ax = fig.add_subplot(223, aspect="equal")
gridplot(np.ma.masked_invalid(linear(vals)), "Linear interpolation")
ax = fig.add_subplot(224, aspect="equal")
gridplot(ok(vals), "Ordinary Kriging")
plt.tight_layout()

Using the convenience function ipol.interpolation in order to deal with missing values

(1) Exemplified for one dimension in space and two dimensions of the source value array (could e.g. be two time steps).

# Synthetic observations (e.g. two time steps)
src = np.arange(10)[:, None]
vals = np.hstack((1.0 + np.sin(src), 5.0 + 2.0 * np.sin(src)))
# Target coordinates
trg = np.linspace(0, 20, 100)[:, None]
# Here we introduce missing values in the second dimension of the source value array
vals[3:5, 1] = np.nan
# interpolation using the convenience function "interpolate"
idw_result = wrl.ipol.interpolate(src, trg, vals, wrl.ipol.Idw, nnearest=4)
nn_result = wrl.ipol.interpolate(src, trg, vals, wrl.ipol.Nearest)
# Plot results
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111)
pl1 = ax.plot(trg, idw_result, "b-", label="IDW")
pl2 = ax.plot(trg, nn_result, "k-", label="Nearest Neighbour")
pl3 = ax.plot(src, vals, "ro", label="Observations")

(2) Exemplified for two dimensions in space and two dimensions of the source value array (e.g. time steps), containing also NaN values (here we only use IDW interpolation)

# Just a helper function for repeated subplots
def plotall(ax, trgx, trgy, src, interp, pts, title, vmin, vmax):
    ix = np.where(np.isfinite(pts))
    ax.pcolormesh(
        trgx, trgy, interp.reshape((len(trgx), len(trgy))), vmin=vmin, vmax=vmax
    )
    ax.scatter(
        src[ix, 0].ravel(),
        src[ix, 1].ravel(),
        c=pts.ravel()[ix],
        s=20,
        marker="s",
        vmin=vmin,
        vmax=vmax,
    )
    ax.set_title(title)
    plt.axis("tight")

# Synthetic observations
src = np.vstack((np.array([4, 7, 3, 15]), np.array([8, 18, 17, 3]))).T
np.random.seed(1319622840 + 1)
vals = np.round(np.random.uniform(size=(len(src), 2)), 1)

# Target coordinates
trgx = np.linspace(0, 20, 100)
trgy = np.linspace(0, 20, 100)
trg = np.meshgrid(trgx, trgy)
trg = np.vstack((trg[0].ravel(), trg[1].ravel())).transpose()

result = wrl.ipol.interpolate(src, trg, vals, wrl.ipol.Idw, nnearest=4)

# Now introduce NaNs in the observations
vals_with_nan = vals.copy()
vals_with_nan[1, 0] = np.nan
vals_with_nan[1:3, 1] = np.nan
result_with_nan = wrl.ipol.interpolate(
    src, trg, vals_with_nan, wrl.ipol.Idw, nnearest=4
)
vmin = np.concatenate((vals.ravel(), result.ravel())).min()
vmax = np.concatenate((vals.ravel(), result.ravel())).max()

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(221)
plotall(ax, trgx, trgy, src, result[:, 0], vals[:, 0], "1st dim: no NaNs", vmin, vmax)
ax = fig.add_subplot(222)
plotall(ax, trgx, trgy, src, result[:, 1], vals[:, 1], "2nd dim: no NaNs", vmin, vmax)
ax = fig.add_subplot(223)
plotall(
    ax,
    trgx,
    trgy,
    src,
    result_with_nan[:, 0],
    vals_with_nan[:, 0],
    "1st dim: one NaN",
    vmin,
    vmax,
)
ax = fig.add_subplot(224)
plotall(
    ax,
    trgx,
    trgy,
    src,
    result_with_nan[:, 1],
    vals_with_nan[:, 1],
    "2nd dim: two NaN",
    vmin,
    vmax,
)
plt.tight_layout()

How to use interpolation for gridding data in polar coordinates?

Read polar coordinates and corresponding rainfall intensity from file

filename = wrl.util.get_wradlib_data_file("misc/bin_coords_tur.gz")
src = np.loadtxt(filename)

filename = wrl.util.get_wradlib_data_file("misc/polar_R_tur.gz")
vals = np.loadtxt(filename)

Downloading file 'misc/bin_coords_tur.gz' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/misc/bin_coords_tur.gz' to '/home/runner/work/wradlib/wradlib/wradlib-data'.
Downloading file 'misc/polar_R_tur.gz' from 'https://github.com/wradlib/wradlib-data/raw/pooch/data/misc/polar_R_tur.gz' to '/home/runner/work/wradlib/wradlib/wradlib-data'.

src.shape

(46080, 2)

Define target grid coordinates

xtrg = np.linspace(src[:, 0].min(), src[:, 0].max(), 200)
ytrg = np.linspace(src[:, 1].min(), src[:, 1].max(), 200)
trg = np.meshgrid(xtrg, ytrg)
trg = np.vstack((trg[0].ravel(), trg[1].ravel())).T

Linear Interpolation

ip_lin = wrl.ipol.Linear(src, trg)
result_lin = ip_lin(vals.ravel(), fill_value=np.nan)

IDW interpolation

ip_near = wrl.ipol.Nearest(src, trg)
maxdist = trg[1, 0] - trg[0, 0]
result_near = ip_near(vals.ravel(), maxdist=maxdist)

Plot results

fig = plt.figure(figsize=(15, 6))
fig.subplots_adjust(wspace=0.4)
ax = fig.add_subplot(131, aspect="equal")
vals = wrl.georef.create_xarray_dataarray(vals).wrl.georef.georeference()
vals.wrl.vis.plot(ax=ax)
ax = fig.add_subplot(132, aspect="equal")
plt.pcolormesh(xtrg, ytrg, result_lin.reshape((len(xtrg), len(ytrg))))
ax = fig.add_subplot(133, aspect="equal")
plt.pcolormesh(xtrg, ytrg, result_near.reshape((len(xtrg), len(ytrg))))

<matplotlib.collections.QuadMesh at 0x7f10f3542630>