Georeferencing single sweep data and converting to numpy array

Python
2

Hi @atmoboran,

this is the simplest approach to create cartesian representation of the polar data. Note that this is just Nearest neigbour interpolation. More sophisticated interpolation might be used for proper handling.

Upper image shows polar source data, lower image cartesian gridded data.

import numpy as np
import wradlib as wrl
import xarray as xr
import matplotlib.pyplot as plt

fname = "/home/kai/Downloads/ras07-pcpng01_sweeph5onem_dbzh_00-2023112617253300-ess-10410-hd5"

# load data
swp = xr.open_dataset(fname, engine="odim")
swp = swp.set_coords("sweep_mode")
# georeference to wanted projection
crs = wrl.georef.epsg_to_osr(32632)
swp = swp.wrl.georef.georeference(crs=crs)

# plot
swp.DBZH.wrl.vis.plot(vmin=0, vmax=60)

# stack source data
st = swp.stack(onedim=("azimuth", "range"))
src = np.stack([st.x, st.y], axis=-1)

# create target grid
xtrg = np.linspace(swp.x.min().values, swp.x.max().values, 1000)
ytrg = np.linspace(swp.y.min().values, swp.y.max().values, 1000)
mx, my = np.meshgrid(xtrg, ytrg)
mx0 = mx.ravel()
my0 = my.ravel()
trg = np.stack([mx0, my0], axis=-1)

# setup interpolator
ip = wrl.ipol.Nearest(src, trg)
out = ip(st.DBZH.values, maxdist=1500)
out = out.reshape(mx.shape)

# plot
plt.figure()
plt.pcolormesh(mx, my, out, cmap="HomeyerRainbow", vmin=0, vmax=60)
plt.colorbar()

source
target

Update: Please note, that we use Azimuthal Equidistant Projection with the radar at the center (provided by wrl.georef.georeference(). You can also rely on swp.xradar.georeference(), which is implemented in xradar-package. The cartesian resolution here is 1000m.

Update2: We now use a given projection (epsg 32632) instead of the default Azimuthal Equidistant Projection.