# Georeferencing single sweep data and converting to numpy array

**URL:** <https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353>\
**Category:** Python\
**Tags:** wradlib\
**Created:** [November 27, 2023, 2:07pm UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353 "2023-11-27T14:07:44Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![atmoboran](https://avatars.discourse-cdn.com/v4/letter/a/ecccb3/32.png) [@atmoboran](https://openradar.discourse.group/u/atmoboran)\
**Post date:** [November 27, 2023, 2:07pm UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/1 "2023-11-27T14:07:44Z")

</div>

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](https://docs.wradlib.org/en/latest/notebooks/workflow/recipe2.html)) 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)](https://opendata.dwd.de/weather/radar/sites/sweep_pcp_z/ess/hdf5/filter_polarimetric/)).

Thank you and best regards,  
Boran

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [November 28, 2023, 7:59pm UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/2 "2023-11-28T19:59:19Z")

</div>

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.

```python
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](https://global.discourse-cdn.com/free1/uploads/openradar/original/1X/ced6ec3e126ea53a01c2d71bffb87215b851aac1.png)  
 ![target](https://global.discourse-cdn.com/free1/uploads/openradar/original/1X/8b24c513cd63b5bddb6528a96c2ef3deda741b39.png)

Update: Please note, that we use Azimuthal Equidistant Projection with the radar at the center (provided by [wrl.georef.georeference()](https://docs.wradlib.org/en/latest/generated/wradlib.georef.polar.georeference.html). 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.

---

<div class="post-metadata">

**Author:** ![atmoboran](https://avatars.discourse-cdn.com/v4/letter/a/ecccb3/32.png) [@atmoboran](https://openradar.discourse.group/u/atmoboran)\
**Post date:** [December 7, 2023, 12:48pm UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/3 "2023-12-07T12:48:13Z")

</div>

> [@kmuehlbauer](#):
>
> Azimuthal Equidistant Projection

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

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [December 8, 2023, 11:46am UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/4 "2023-12-08T11:46:38Z")

</div>

Hi Boran @atmoboran,

I’ve updated the above example.

HTH,  
Kai

---

<div class="post-metadata">

**Author:** ![atmoboran](https://avatars.discourse-cdn.com/v4/letter/a/ecccb3/32.png) [@atmoboran](https://openradar.discourse.group/u/atmoboran)\
**Post date:** [December 11, 2023, 6:46pm UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/5 "2023-12-11T18:46:21Z")

</div>

Hi Kai, thanks so much!

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

```python
# 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

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [December 12, 2023, 7:57am UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/6 "2023-12-12T07:57:22Z")

</div>

@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

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [December 12, 2023, 9:31am UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/7 "2023-12-12T09:31:10Z")

</div>

3 posts were split to a new topic: [Selecting projection](https://openradar.discourse.group/t/selecting-projection/372)

---

<div class="post-metadata">

**Author:** ![atmoboran](https://avatars.discourse-cdn.com/v4/letter/a/ecccb3/32.png) [@atmoboran](https://openradar.discourse.group/u/atmoboran)\
**Post date:** [June 7, 2024, 9:10am UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/8 "2024-06-07T09:10:30Z")

</div>

@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

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [June 8, 2024, 8:43am UTC](https://openradar.discourse.group/t/georeferencing-single-sweep-data-and-converting-to-numpy-array/353/9 "2024-06-08T08:43:45Z")

</div>

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