Computing cartesian and geographical coordinates for polar data#

[1]:

import numpy as np
import wradlib as wrl
import xradar as xd
import warnings

warnings.filterwarnings("ignore")

Read the data#

Here, we use an OPERA hdf5 dataset.

[2]:

filename = "hdf5/20130429043000.rad.bewid.pvol.dbzh.scan1.hdf"
filename = wrl.util.get_wradlib_data_file(filename)
pvol = xd.io.open_odim_datatree(filename)
display(pvol)

<xarray.DatasetView>
Dimensions:              ()
Data variables:
    volume_number        int64 0
    platform_type        <U5 'fixed'
    instrument_type      <U5 'radar'
    time_coverage_start  <U20 '2013-04-29T04:30:00Z'
    time_coverage_end    <U20 '2013-04-29T04:31:39Z'
    longitude            float64 5.506
    altitude             float64 592.0
    latitude             float64 49.91
Attributes:
    Conventions:      ODIM_H5/V2_1
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using xradar
    instrument_name:  None

Retrieve azimuthal equidistant coordinates and projection#

[3]:

for key in list(pvol.children):
    if "sweep" in key:
        pvol[key].ds = pvol[key].ds.wrl.georef.georeference()

[4]:

pvol["sweep_0"].ds.DBZH.plot(x="x", y="y")

[4]:

<matplotlib.collections.QuadMesh at 0x7fa2341cd610>
../../_images/notebooks_georeferencing_coords_7_1.png

Retrieve geographic coordinates (longitude and latitude)#

Using crs-keyword argument.#

[5]:

for key in list(pvol.children):
    if "sweep" in key:
        pvol[key].ds = pvol[key].ds.wrl.georef.georeference(
            crs=wrl.georef.get_default_projection()
        )

[6]:

ds1 = pvol["sweep_0"].ds.wrl.georef.georeference(
    crs=wrl.georef.get_default_projection()
)
ds1.DBZH.plot(x="x", y="y")

[6]:

<matplotlib.collections.QuadMesh at 0x7fa227fc6850>
../../_images/notebooks_georeferencing_coords_11_1.png

Using reproject#

[7]:

ds2 = pvol["sweep_0"].ds.wrl.georef.reproject(
    trg_crs=wrl.georef.epsg_to_osr(32632),
)
ds2.DBZH.plot(x="x", y="y")

[7]:

<matplotlib.collections.QuadMesh at 0x7fa2250e1d90>
../../_images/notebooks_georeferencing_coords_13_1.png