Example for georeferencing a radar dataset#

import warnings

import matplotlib as mpl
import matplotlib.pyplot as plt
import wradlib as wrl
import xradar as xd
from matplotlib.patches import Rectangle

warnings.filterwarnings("ignore")

1st step: Compute centroid coordinates and vertices of all radar bins in WGS84 (longitude and latitude).

We can now generate the polygon vertices of the radar bins - with each vertex in lon/lat coordinates.

… or we can compute the corresponding centroids of all bins - - with each centroid in lon/lat coordinates.

In order to understand how vertices and centroids correspond, we can plot them together.

fig = plt.figure(figsize=(16, 16))
site = (polygons.longitude.values, polygons.latitude.values)

aspect = (centroids[..., 0].max() - centroids[..., 0].min()) / (
    centroids[..., 1].max() - centroids[..., 1].min()
)
ax = fig.add_subplot(121, aspect=aspect)
polycoll = mpl.collections.PolyCollection(
    polygons.isel(xy=slice(0, 2)), closed=True, facecolors="None", linewidth=0.1
)
ax.add_collection(polycoll, autolim=True)
# ax.plot(centroids[..., 0], centroids[..., 1], 'r+')
plt.title("Zoom in\n(only possible for interactive plots).")
ax.add_patch(
    Rectangle(
        (site[0] + 0.25, site[1] + 0.25),
        0.2,
        0.2 / aspect,
        edgecolor="red",
        facecolor="None",
        zorder=3,
    )
)
plt.xlim(centroids[..., 0].min(), centroids[..., 0].max())
plt.ylim(centroids[..., 1].min(), centroids[..., 1].max())

ax = fig.add_subplot(122, aspect=aspect)
polycoll = mpl.collections.PolyCollection(
    polygons.isel(xy=slice(0, 2)), closed=True, facecolors="None"
)
ax.add_collection(polycoll, autolim=True)
ax.plot(centroids[..., 0], centroids[..., 1], "r+")
plt.title("Zoom into red box of left plot")
plt.xlim(site[0] + 0.25, site[0] + 0.25 + 0.2)
plt.ylim(site[1] + 0.25, site[1] + 0.25 + 0.2 / aspect)

(46.422541, 46.56147892646726)

../../_images/bc5ae30f969e4c99bc7f7b0264c9839f5a1716b1c16d07fbaa3dfee95a923a2e.png

2nd step: Reproject the centroid coordinates to Gauss-Krueger Zone 3 (i.e. EPSG-Code 31467).