Computing cartesian and geographical coordinates for polar data

import warnings

import wradlib as wrl
import wradlib_data
import xradar as xd
from IPython.display import display

warnings.filterwarnings("ignore")

Read the data

Here, we use an OPERA hdf5 dataset.

filename = "hdf5/20130429043000.rad.bewid.pvol.dbzh.scan1.hdf"
filename = wradlib_data.DATASETS.fetch(filename)
pvol = xd.io.open_odim_datatree(filename)
display(pvol)

<xarray.DataTree>
Group: /
│   Dimensions:              (sweep: 5)
│   Coordinates:
│       latitude             float64 8B ...
│       longitude            float64 8B ...
│       altitude             float64 8B ...
│   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'
│       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
├── Group: /sweep_0
│       Dimensions:            (azimuth: 360, range: 960)
│       Coordinates:
│         * azimuth            (azimuth) float32 1kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
│           elevation          (azimuth) float64 3kB ...
│           time               (azimuth) datetime64[ns] 3kB 2013-04-29T04:30:00.02777...
│         * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
│       Data variables:
│           DBZH               (azimuth, range) float64 3MB ...
│           sweep_mode         <U20 80B ...
│           sweep_number       int64 8B ...
│           prt_mode           <U7 28B ...
│           follow_mode        <U7 28B ...
│           sweep_fixed_angle  float64 8B ...
│           nyquist_velocity   float64 8B ...
├── Group: /sweep_1
│       Dimensions:            (azimuth: 360, range: 960)
│       Coordinates:
│         * azimuth            (azimuth) float32 1kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
│           elevation          (azimuth) float64 3kB ...
│           time               (azimuth) datetime64[ns] 3kB 2013-04-29T04:30:20.02777...
│         * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
│       Data variables:
│           DBZH               (azimuth, range) float64 3MB ...
│           sweep_mode         <U20 80B ...
│           sweep_number       int64 8B ...
│           prt_mode           <U7 28B ...
│           follow_mode        <U7 28B ...
│           sweep_fixed_angle  float64 8B ...
│           nyquist_velocity   float64 8B ...
├── Group: /sweep_2
│       Dimensions:            (azimuth: 360, range: 960)
│       Coordinates:
│         * azimuth            (azimuth) float32 1kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
│           elevation          (azimuth) float64 3kB ...
│           time               (azimuth) datetime64[ns] 3kB 2013-04-29T04:30:40.02777...
│         * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
│       Data variables:
│           DBZH               (azimuth, range) float64 3MB ...
│           sweep_mode         <U20 80B ...
│           sweep_number       int64 8B ...
│           prt_mode           <U7 28B ...
│           follow_mode        <U7 28B ...
│           sweep_fixed_angle  float64 8B ...
│           nyquist_velocity   float64 8B ...
├── Group: /sweep_3
│       Dimensions:            (azimuth: 360, range: 960)
│       Coordinates:
│         * azimuth            (azimuth) float32 1kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
│           elevation          (azimuth) float64 3kB ...
│           time               (azimuth) datetime64[ns] 3kB 2013-04-29T04:31:00.02777...
│         * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
│       Data variables:
│           DBZH               (azimuth, range) float64 3MB ...
│           sweep_mode         <U20 80B ...
│           sweep_number       int64 8B ...
│           prt_mode           <U7 28B ...
│           follow_mode        <U7 28B ...
│           sweep_fixed_angle  float64 8B ...
│           nyquist_velocity   float64 8B ...
└── Group: /sweep_4
        Dimensions:            (azimuth: 360, range: 960)
        Coordinates:
          * azimuth            (azimuth) float32 1kB 0.5 1.5 2.5 ... 357.5 358.5 359.5
            elevation          (azimuth) float64 3kB ...
            time               (azimuth) datetime64[ns] 3kB 2013-04-29T04:31:20.02777...
          * range              (range) float32 4kB 125.0 375.0 ... 2.396e+05 2.399e+05
        Data variables:
            DBZH               (azimuth, range) float64 3MB ...
            sweep_mode         <U20 80B ...
            sweep_number       int64 8B ...
            prt_mode           <U7 28B ...
            follow_mode        <U7 28B ...
            sweep_fixed_angle  float64 8B ...
            nyquist_velocity   float64 8B ...

xarray.DataTree

• /sweep_0(11)

  • Dimensions:
    • sweep: 5
    • azimuth: 360
    • range: 960
  • Coordinates: (4)
    • azimuth
      (azimuth)
      float32
      0.5 1.5 2.5 ... 357.5 358.5 359.5

      standard_name :

      ray_azimuth_angle

      long_name :

      azimuth_angle_from_true_north

      units :

      degrees

      axis :

      radial_azimuth_coordinate

      array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5],shape=(360,), dtype=float32)

    • elevation
      (azimuth)
      float64
      ...

      standard_name :

      ray_elevation_angle

      long_name :

      elevation_angle_from_horizontal_plane

      units :

      degrees

      axis :

      radial_elevation_coordinate

      [360 values with dtype=float64]

    • time
      (azimuth)
      datetime64[ns]
      2013-04-29T04:30:00.027777792 .....

      standard_name :

      time

      array(['2013-04-29T04:30:00.027777792', '2013-04-29T04:30:00.083333376','2013-04-29T04:30:00.138888960', ..., '2013-04-29T04:30:19.861120512','2013-04-29T04:30:19.916676096', '2013-04-29T04:30:19.972231680'],shape=(360,), dtype='datetime64[ns]')

    • range
      (range)
      float32
      125.0 375.0 ... 2.396e+05 2.399e+05

      units :

      meters

      standard_name :

      projection_range_coordinate

      long_name :

      range_to_measurement_volume

      axis :

      radial_range_coordinate

      meters_between_gates :

      250.0

      spacing_is_constant :

      true

      meters_to_center_of_first_gate :

      125.0

      array([1.25000e+02, 3.75000e+02, 6.25000e+02, ..., 2.39375e+05, 2.39625e+05,2.39875e+05], shape=(960,), dtype=float32)

  • Data variables: (7)
    • DBZH
      (azimuth, range)
      float64
      ...

      _Undetect :

      0.0

      long_name :

      Equivalent reflectivity factor H

      units :

      dBZ

      standard_name :

      radar_equivalent_reflectivity_factor_h

      [345600 values with dtype=float64]

    • sweep_mode
      ()
      <U20
      ...

      [1 values with dtype=<U20]

    • sweep_number
      ()
      int64
      ...

      [1 values with dtype=int64]

    • prt_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • follow_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • sweep_fixed_angle
      ()
      float64
      ...

      [1 values with dtype=float64]

    • nyquist_velocity
      ()
      float64
      ...

      standard_name :

      nyquist_velocity

      units :

      m s-1

      [1 values with dtype=float64]

  /sweep_1(11)

  • Dimensions:
    • sweep: 5
    • azimuth: 360
    • range: 960
  • Coordinates: (4)
    • azimuth
      (azimuth)
      float32
      0.5 1.5 2.5 ... 357.5 358.5 359.5

      standard_name :

      ray_azimuth_angle

      long_name :

      azimuth_angle_from_true_north

      units :

      degrees

      axis :

      radial_azimuth_coordinate

      array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5],shape=(360,), dtype=float32)

    • elevation
      (azimuth)
      float64
      ...

      standard_name :

      ray_elevation_angle

      long_name :

      elevation_angle_from_horizontal_plane

      units :

      degrees

      axis :

      radial_elevation_coordinate

      [360 values with dtype=float64]

    • time
      (azimuth)
      datetime64[ns]
      2013-04-29T04:30:20.027777792 .....

      standard_name :

      time

      array(['2013-04-29T04:30:20.027777792', '2013-04-29T04:30:20.083333376','2013-04-29T04:30:20.138888960', ..., '2013-04-29T04:30:39.861120512','2013-04-29T04:30:39.916676096', '2013-04-29T04:30:39.972231680'],shape=(360,), dtype='datetime64[ns]')

    • range
      (range)
      float32
      125.0 375.0 ... 2.396e+05 2.399e+05

      units :

      meters

      standard_name :

      projection_range_coordinate

      long_name :

      range_to_measurement_volume

      axis :

      radial_range_coordinate

      meters_between_gates :

      250.0

      spacing_is_constant :

      true

      meters_to_center_of_first_gate :

      125.0

      array([1.25000e+02, 3.75000e+02, 6.25000e+02, ..., 2.39375e+05, 2.39625e+05,2.39875e+05], shape=(960,), dtype=float32)

  • Data variables: (7)
    • DBZH
      (azimuth, range)
      float64
      ...

      _Undetect :

      0.0

      long_name :

      Equivalent reflectivity factor H

      units :

      dBZ

      standard_name :

      radar_equivalent_reflectivity_factor_h

      [345600 values with dtype=float64]

    • sweep_mode
      ()
      <U20
      ...

      [1 values with dtype=<U20]

    • sweep_number
      ()
      int64
      ...

      [1 values with dtype=int64]

    • prt_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • follow_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • sweep_fixed_angle
      ()
      float64
      ...

      [1 values with dtype=float64]

    • nyquist_velocity
      ()
      float64
      ...

      standard_name :

      nyquist_velocity

      units :

      m s-1

      [1 values with dtype=float64]

  /sweep_2(11)

  • Dimensions:
    • sweep: 5
    • azimuth: 360
    • range: 960
  • Coordinates: (4)
    • azimuth
      (azimuth)
      float32
      0.5 1.5 2.5 ... 357.5 358.5 359.5

      standard_name :

      ray_azimuth_angle

      long_name :

      azimuth_angle_from_true_north

      units :

      degrees

      axis :

      radial_azimuth_coordinate

      array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5],shape=(360,), dtype=float32)

    • elevation
      (azimuth)
      float64
      ...

      standard_name :

      ray_elevation_angle

      long_name :

      elevation_angle_from_horizontal_plane

      units :

      degrees

      axis :

      radial_elevation_coordinate

      [360 values with dtype=float64]

    • time
      (azimuth)
      datetime64[ns]
      2013-04-29T04:30:40.027777792 .....

      standard_name :

      time

      array(['2013-04-29T04:30:40.027777792', '2013-04-29T04:30:40.083333376','2013-04-29T04:30:40.138888960', ..., '2013-04-29T04:30:59.861120512','2013-04-29T04:30:59.916676096', '2013-04-29T04:30:59.972231680'],shape=(360,), dtype='datetime64[ns]')

    • range
      (range)
      float32
      125.0 375.0 ... 2.396e+05 2.399e+05

      units :

      meters

      standard_name :

      projection_range_coordinate

      long_name :

      range_to_measurement_volume

      axis :

      radial_range_coordinate

      meters_between_gates :

      250.0

      spacing_is_constant :

      true

      meters_to_center_of_first_gate :

      125.0

      array([1.25000e+02, 3.75000e+02, 6.25000e+02, ..., 2.39375e+05, 2.39625e+05,2.39875e+05], shape=(960,), dtype=float32)

  • Data variables: (7)
    • DBZH
      (azimuth, range)
      float64
      ...

      _Undetect :

      0.0

      long_name :

      Equivalent reflectivity factor H

      units :

      dBZ

      standard_name :

      radar_equivalent_reflectivity_factor_h

      [345600 values with dtype=float64]

    • sweep_mode
      ()
      <U20
      ...

      [1 values with dtype=<U20]

    • sweep_number
      ()
      int64
      ...

      [1 values with dtype=int64]

    • prt_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • follow_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • sweep_fixed_angle
      ()
      float64
      ...

      [1 values with dtype=float64]

    • nyquist_velocity
      ()
      float64
      ...

      standard_name :

      nyquist_velocity

      units :

      m s-1

      [1 values with dtype=float64]

  /sweep_3(11)

  • Dimensions:
    • sweep: 5
    • azimuth: 360
    • range: 960
  • Coordinates: (4)
    • azimuth
      (azimuth)
      float32
      0.5 1.5 2.5 ... 357.5 358.5 359.5

      standard_name :

      ray_azimuth_angle

      long_name :

      azimuth_angle_from_true_north

      units :

      degrees

      axis :

      radial_azimuth_coordinate

      array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5],shape=(360,), dtype=float32)

    • elevation
      (azimuth)
      float64
      ...

      standard_name :

      ray_elevation_angle

      long_name :

      elevation_angle_from_horizontal_plane

      units :

      degrees

      axis :

      radial_elevation_coordinate

      [360 values with dtype=float64]

    • time
      (azimuth)
      datetime64[ns]
      2013-04-29T04:31:00.027777792 .....

      standard_name :

      time

      array(['2013-04-29T04:31:00.027777792', '2013-04-29T04:31:00.083333376','2013-04-29T04:31:00.138888960', ..., '2013-04-29T04:31:19.861120512','2013-04-29T04:31:19.916676096', '2013-04-29T04:31:19.972231680'],shape=(360,), dtype='datetime64[ns]')

    • range
      (range)
      float32
      125.0 375.0 ... 2.396e+05 2.399e+05

      units :

      meters

      standard_name :

      projection_range_coordinate

      long_name :

      range_to_measurement_volume

      axis :

      radial_range_coordinate

      meters_between_gates :

      250.0

      spacing_is_constant :

      true

      meters_to_center_of_first_gate :

      125.0

      array([1.25000e+02, 3.75000e+02, 6.25000e+02, ..., 2.39375e+05, 2.39625e+05,2.39875e+05], shape=(960,), dtype=float32)

  • Data variables: (7)
    • DBZH
      (azimuth, range)
      float64
      ...

      _Undetect :

      0.0

      long_name :

      Equivalent reflectivity factor H

      units :

      dBZ

      standard_name :

      radar_equivalent_reflectivity_factor_h

      [345600 values with dtype=float64]

    • sweep_mode
      ()
      <U20
      ...

      [1 values with dtype=<U20]

    • sweep_number
      ()
      int64
      ...

      [1 values with dtype=int64]

    • prt_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • follow_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • sweep_fixed_angle
      ()
      float64
      ...

      [1 values with dtype=float64]

    • nyquist_velocity
      ()
      float64
      ...

      standard_name :

      nyquist_velocity

      units :

      m s-1

      [1 values with dtype=float64]

  /sweep_4(11)

  • Dimensions:
    • sweep: 5
    • azimuth: 360
    • range: 960
  • Coordinates: (4)
    • azimuth
      (azimuth)
      float32
      0.5 1.5 2.5 ... 357.5 358.5 359.5

      standard_name :

      ray_azimuth_angle

      long_name :

      azimuth_angle_from_true_north

      units :

      degrees

      axis :

      radial_azimuth_coordinate

      array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5],shape=(360,), dtype=float32)

    • elevation
      (azimuth)
      float64
      ...

      standard_name :

      ray_elevation_angle

      long_name :

      elevation_angle_from_horizontal_plane

      units :

      degrees

      axis :

      radial_elevation_coordinate

      [360 values with dtype=float64]

    • time
      (azimuth)
      datetime64[ns]
      2013-04-29T04:31:20.027777792 .....

      standard_name :

      time

      array(['2013-04-29T04:31:20.027777792', '2013-04-29T04:31:20.083333376','2013-04-29T04:31:20.138888960', ..., '2013-04-29T04:31:39.861120512','2013-04-29T04:31:39.916676096', '2013-04-29T04:31:39.972231680'],shape=(360,), dtype='datetime64[ns]')

    • range
      (range)
      float32
      125.0 375.0 ... 2.396e+05 2.399e+05

      units :

      meters

      standard_name :

      projection_range_coordinate

      long_name :

      range_to_measurement_volume

      axis :

      radial_range_coordinate

      meters_between_gates :

      250.0

      spacing_is_constant :

      true

      meters_to_center_of_first_gate :

      125.0

      array([1.25000e+02, 3.75000e+02, 6.25000e+02, ..., 2.39375e+05, 2.39625e+05,2.39875e+05], shape=(960,), dtype=float32)

  • Data variables: (7)
    • DBZH
      (azimuth, range)
      float64
      ...

      _Undetect :

      0.0

      long_name :

      Equivalent reflectivity factor H

      units :

      dBZ

      standard_name :

      radar_equivalent_reflectivity_factor_h

      [345600 values with dtype=float64]

    • sweep_mode
      ()
      <U20
      ...

      [1 values with dtype=<U20]

    • sweep_number
      ()
      int64
      ...

      [1 values with dtype=int64]

    • prt_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • follow_mode
      ()
      <U7
      ...

      [1 values with dtype=<U7]

    • sweep_fixed_angle
      ()
      float64
      ...

      [1 values with dtype=float64]

    • nyquist_velocity
      ()
      float64
      ...

      standard_name :

      nyquist_velocity

      units :

      m s-1

      [1 values with dtype=float64]

• Dimensions:
  • sweep: 5
• Coordinates: (3)
  • latitude
    ()
    float64
    ...

    long_name :

    latitude

    units :

    degrees_north

    positive :

    up

    standard_name :

    latitude

    [1 values with dtype=float64]

  • longitude
    ()
    float64
    ...

    long_name :

    longitude

    units :

    degrees_east

    standard_name :

    longitude

    [1 values with dtype=float64]

  • altitude
    ()
    float64
    ...

    long_name :

    altitude

    units :

    meters

    standard_name :

    altitude

    [1 values with dtype=float64]

• Data variables: (7)
  • volume_number
    ()
    int64
    0

    array(0)

  • platform_type
    ()
    <U5
    'fixed'

    array('fixed', dtype='<U5')

  • instrument_type
    ()
    <U5
    'radar'

    array('radar', dtype='<U5')

  • time_coverage_start
    ()
    <U20
    '2013-04-29T04:30:00Z'

    array('2013-04-29T04:30:00Z', dtype='<U20')

  • time_coverage_end
    ()
    <U20
    '2013-04-29T04:31:39Z'

    array('2013-04-29T04:31:39Z', dtype='<U20')

  • sweep_fixed_angle
    (sweep)
    float64
    0.3 0.9 1.8 3.3 6.0

    array([0.3, 0.9, 1.8, 3.3, 6. ])

  • sweep_group_name
    (sweep)
    int64
    0 1 2 3 4

    array([0, 1, 2, 3, 4])

• Attributes: (9)

  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

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

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

<matplotlib.collections.QuadMesh at 0x7f977b0edd30>

Retrieve geographic coordinates (longitude and latitude)

Using crs-keyword argument.

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()
        )

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

<matplotlib.collections.QuadMesh at 0x7f9772fa02d0>

Using reproject

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

<matplotlib.collections.QuadMesh at 0x7f97797f7390>