Computing cartesian and geographical coordinates for polar data

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> Size: 312B
Dimensions:              (sweep: 5)
Dimensions without coordinates: sweep
Data variables:
    volume_number        int64 8B 0
    platform_type        <U5 20B 'fixed'
    instrument_type      <U5 20B 'radar'
    time_coverage_start  <U20 80B '2013-04-29T04:30:00Z'
    time_coverage_end    <U20 80B '2013-04-29T04:31:39Z'
    longitude            float64 8B 5.506
    altitude             float64 8B 592.0
    latitude             float64 8B 49.91
    sweep_fixed_angle    (sweep) float64 40B 0.3 0.9 1.8 3.3 6.0
    sweep_group_name     (sweep) int64 40B 0 1 2 3 4
Attributes:
    Conventions:      ODIM_H5/V2_2
    instrument_name:  None
    version:          None
    title:            None
    institution:      None
    references:       None
    source:           None
    history:          None
    comment:          im/exported using xradar

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 0x7fcbc373f110>
../../_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 0x7fcbc4c1dcd0>
../../_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 0x7fcbb98ad250>
../../_images/notebooks_georeferencing_coords_13_1.png