Georeferencing single sweep data and converting to numpy array

Python
1

Hello,

Is there any recipe on how to (1.) georeference single sweep (!) data,(2.) transform it to cartesian coordinates and then (3.) to a numpy array?

Basically what happened here ( Recipe #2: Reading and visualizing an ODIM_H5 polar volume — wradlib) under [8], but for files that include only one sweep, i.e. 2D files.
To be more precise: the data I want to use this on is the precipitation scan from DWD ( Index of /weather/radar/sites/sweep_pcp_z/ess/hdf5/filter_polarimetric/ (dwd.de)).

Thank you and best regards,
Boran

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.

3

Hey Kai, thank you very much!
This worked well.

But now I ran into a subsequent issue: I would like to have the cartesian data on some coordinate system/ projection (e.g. EPSG 32632), so it is comparable with other data (like the CAPPI from the volume scans or satellite data).

I’m sorry if this is a stupid question, but I really lack knowledge in this whole projecting/ georeferencing/ coordinate system topic.

Thanks,
Boran

4

Hi Boran @atmoboran,

I’ve updated the above example.

HTH,
Kai

5

Hi Kai, thanks so much!

I’ve also done this with the default WG84 projection in lon lat coordinates:

# georeference 
swp = swp.wrl.georef.georeference(crs = wrl.georef.get_default_projection()) # default = WGS 84 in lon lat

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

# create target grid
ytrg = np.arange(swp.y.min().values, swp.y.max().values, 0.01) # 0.01 is the resolution in °
xtrg = np.arange(swp.x.min().values, swp.x.max().values, 0.01)
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=0.02)
out = out.reshape(mx.shape)

Is this done correctly?

Best,

Boran

6

@atmoboran

I’ve added stacked source data and made the code coyp&paste-able. All in all this looks OK, but you do not have an equal-area projection anymore. Is this just for display or is this used for further processing? In that case I’d suggest to stick with an equal area projection.

HTH,
Kai

7

3 posts were split to a new topic: Selecting projection

8

@kmuehlbauer
Hello again,
just want to make sure:
after georeferencing the coordinate system is on ground, i.e. the given distances from radar refer to the distance on the ground and not along the beam?!
Best regards,
Boran

9

Hi @atmoboran,

Yes, after georeferencing the created x,y,z coordinates are in Azimuthal Equidistant (default) or any other given projection with the radar at it’s center.

The slant range (along ray) is still available in the range coordinate.

HTH,
Kai