xarray RADOLAN backend

In this example, we read RADOLAN data files using the xarray radolan backend.

import glob
import os
import wradlib as wrl
import warnings
warnings.filterwarnings('ignore')
import matplotlib.pyplot as pl
import numpy as np
import xarray as xr
try:
    get_ipython().magic("matplotlib inline")
except:
    pl.ion()

Load RADOLAN Data

fpath = 'radolan/misc/raa01-rw_10000-1408030950-dwd---bin.gz'
f = wrl.util.get_wradlib_data_file(fpath)
comp = wrl.io.open_radolan_dataset(f)

Inspect Data

display(comp)

<xarray.Dataset>
Dimensions:  (y: 900, x: 900, time: 1)
Coordinates:
  * time     (time) datetime64[ns] 2014-08-03T09:50:00
  * y        (y) float64 -4.658e+06 -4.657e+06 ... -3.76e+06 -3.759e+06
  * x        (x) float64 -5.23e+05 -5.22e+05 -5.21e+05 ... 3.75e+05 3.76e+05
Data variables:
    RW       (y, x) float32 ...
Attributes:
    radarid:         10000
    formatversion:   3
    radolanversion:  2.13.1
    radarlocations:  ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu',...

xarray.Dataset

• Dimensions:
  • y: 900
  • x: 900
  • time: 1
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2014-08-03T09:50:00

    standard_name :

    time

    array(['2014-08-03T09:50:00.000000000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    -4.658e+06 ... -3.759e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([-4658144.724266, -4657144.724266, -4656144.724266, ..., -3761144.724266,
           -3760144.724266, -3759144.724266])

  • x
    (x)
    float64
    -5.23e+05 -5.22e+05 ... 3.76e+05

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-522962.166922, -521962.166922, -520962.166922, ...,  374037.833078,
            375037.833078,  376037.833078])

• Data variables: (1)
  • RW
    (y, x)
    float32
    ...

    valid_min :

    0

    valid_max :

    4095

    standard_name :

    rainfall_rate

    long_name :

    RW

    unit :

    mm h-1

    [810000 values with dtype=float32]

• Attributes: (4)

  radarid :

  10000

  formatversion :

  3

  radolanversion :

  2.13.1

  radarlocations :

  ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu', 'nhb', 'oft', 'tur', 'isn', 'fbg', 'mem']

Plotting

comp.RW.plot.pcolormesh(x='x', y='y')
pl.gca().set_aspect('equal')

Inspect RADOLAN moments

The DataArrays can be accessed by key or by attribute. Each DataArray has dimensions and coordinates of it’s parent dataset.

display(comp.RW)

<xarray.DataArray 'RW' (y: 900, x: 900)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
Coordinates:
  * y        (y) float64 -4.658e+06 -4.657e+06 ... -3.76e+06 -3.759e+06
  * x        (x) float64 -5.23e+05 -5.22e+05 -5.21e+05 ... 3.75e+05 3.76e+05
Attributes:
    valid_min:      0
    valid_max:      4095
    standard_name:  rainfall_rate
    long_name:      RW
    unit:           mm h-1

xarray.DataArray

'RW'

• y: 900
• x: 900

• nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

  array([[nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         ...,
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

• Coordinates: (2)
  • y
    (y)
    float64
    -4.658e+06 ... -3.759e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([-4658144.724266, -4657144.724266, -4656144.724266, ..., -3761144.724266,
           -3760144.724266, -3759144.724266])

  • x
    (x)
    float64
    -5.23e+05 -5.22e+05 ... 3.76e+05

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-522962.166922, -521962.166922, -520962.166922, ...,  374037.833078,
            375037.833078,  376037.833078])

• Attributes: (5)

  valid_min :

  0

  valid_max :

  4095

  standard_name :

  rainfall_rate

  long_name :

  RW

  unit :

  mm h-1

Create simple plot

Using xarray features a simple plot can be created like this.

comp.RW.plot(x="x", y="y", add_labels=False)

<matplotlib.collections.QuadMesh at 0x7fcf4f901750>

fig = pl.figure(figsize=(10,10))
ax = fig.add_subplot(111)
comp.RW.plot(x="x", y="y", ax=ax)

<matplotlib.collections.QuadMesh at 0x7fcf4f7bcca0>

Mask some values

ds = comp['RW'].where(comp['RW'] >= 1)
ds.plot()

<matplotlib.collections.QuadMesh at 0x7fcf4f8a83a0>

Export to NetCDF4

# fix _FillValue
comp.RW.encoding["_FillValue"] = 65535

comp.to_netcdf("test_radolan.nc")

Import again

comp1 = xr.open_dataset("test_radolan.nc")
comp1

<xarray.Dataset>
Dimensions:  (y: 900, x: 900, time: 1)
Coordinates:
  * time     (time) datetime64[ns] 2014-08-03T09:50:00
  * y        (y) float64 -4.658e+06 -4.657e+06 ... -3.76e+06 -3.759e+06
  * x        (x) float64 -5.23e+05 -5.22e+05 -5.21e+05 ... 3.75e+05 3.76e+05
Data variables:
    RW       (y, x) float32 ...
Attributes:
    radarid:         10000
    formatversion:   3
    radolanversion:  2.13.1
    radarlocations:  ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu',...

xarray.Dataset

• Dimensions:
  • y: 900
  • x: 900
  • time: 1
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2014-08-03T09:50:00

    standard_name :

    time

    array(['2014-08-03T09:50:00.000000000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    -4.658e+06 ... -3.759e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([-4658144.724266, -4657144.724266, -4656144.724266, ..., -3761144.724266,
           -3760144.724266, -3759144.724266])

  • x
    (x)
    float64
    -5.23e+05 -5.22e+05 ... 3.76e+05

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-522962.166922, -521962.166922, -520962.166922, ...,  374037.833078,
            375037.833078,  376037.833078])

• Data variables: (1)
  • RW
    (y, x)
    float32
    ...

    valid_min :

    0

    valid_max :

    4095

    standard_name :

    rainfall_rate

    long_name :

    RW

    unit :

    mm h-1

    [810000 values with dtype=float32]

• Attributes: (4)

  radarid :

  10000

  formatversion :

  3

  radolanversion :

  2.13.1

  radarlocations :

  ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu', 'nhb', 'oft', 'tur', 'isn', 'fbg', 'mem']

Check equality

xr.testing.assert_equal(comp, comp1)

More RADOLAN loading mechanisms

Use xr.open_dataset

comp2 = xr.open_dataset(f, engine="radolan")
display(comp2)

<xarray.Dataset>
Dimensions:  (y: 900, x: 900, time: 1)
Coordinates:
  * time     (time) datetime64[ns] 2014-08-03T09:50:00
  * y        (y) float64 -4.658e+06 -4.657e+06 ... -3.76e+06 -3.759e+06
  * x        (x) float64 -5.23e+05 -5.22e+05 -5.21e+05 ... 3.75e+05 3.76e+05
Data variables:
    RW       (y, x) float32 ...
Attributes:
    radarid:         10000
    formatversion:   3
    radolanversion:  2.13.1
    radarlocations:  ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu',...

xarray.Dataset

• Dimensions:
  • y: 900
  • x: 900
  • time: 1
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2014-08-03T09:50:00

    standard_name :

    time

    array(['2014-08-03T09:50:00.000000000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    -4.658e+06 ... -3.759e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([-4658144.724266, -4657144.724266, -4656144.724266, ..., -3761144.724266,
           -3760144.724266, -3759144.724266])

  • x
    (x)
    float64
    -5.23e+05 -5.22e+05 ... 3.76e+05

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-522962.166922, -521962.166922, -520962.166922, ...,  374037.833078,
            375037.833078,  376037.833078])

• Data variables: (1)
  • RW
    (y, x)
    float32
    ...

    valid_min :

    0

    valid_max :

    4095

    standard_name :

    rainfall_rate

    long_name :

    RW

    unit :

    mm h-1

    [810000 values with dtype=float32]

• Attributes: (4)

  radarid :

  10000

  formatversion :

  3

  radolanversion :

  2.13.1

  radarlocations :

  ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu', 'nhb', 'oft', 'tur', 'isn', 'fbg', 'mem']

Use xr.open_mfdataset to retrieve timeseries

#fpath = 'radolan/misc/raa01-rw_10000-1408030950-dwd---bin.gz'
fpath = wrl.util.get_wradlib_data_path()
f = os.path.join(fpath, 'radolan/misc/raa01-sf_10000-1305*.gz')

comp3 = xr.open_mfdataset(f, engine="radolan")
display(comp3)

<xarray.Dataset>
Dimensions:  (time: 2, y: 900, x: 900)
Coordinates:
  * time     (time) datetime64[ns] 2013-05-27T00:50:00 2013-05-28T00:50:00
  * y        (y) float64 -4.658e+06 -4.657e+06 ... -3.76e+06 -3.759e+06
  * x        (x) float64 -5.23e+05 -5.22e+05 -5.21e+05 ... 3.75e+05 3.76e+05
Data variables:
    SF       (time, y, x) float32 dask.array<chunksize=(1, 900, 900), meta=np.ndarray>
Attributes:
    radarid:         10000
    formatversion:   3
    radolanversion:  2.12.0
    radarlocations:  ['ham', 'ros', 'emd', 'han', 'bln', 'ess', 'fld', 'drs',...
    radardays:       ['bln 24', 'drs 24', 'eis 24', 'emd 24', 'ess 24', 'fld ...

xarray.Dataset

• Dimensions:
  • time: 2
  • y: 900
  • x: 900
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2013-05-27T00:50:00 2013-05-28T0...

    standard_name :

    time

    array(['2013-05-27T00:50:00.000000000', '2013-05-28T00:50:00.000000000'],
          dtype='datetime64[ns]')

  • y
    (y)
    float64
    -4.658e+06 ... -3.759e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([-4658144.724266, -4657144.724266, -4656144.724266, ..., -3761144.724266,
           -3760144.724266, -3759144.724266])

  • x
    (x)
    float64
    -5.23e+05 -5.22e+05 ... 3.76e+05

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-522962.166922, -521962.166922, -520962.166922, ...,  374037.833078,
            375037.833078,  376037.833078])

• Data variables: (1)
  • SF
    (time, y, x)
    float32
    dask.array<chunksize=(1, 900, 900), meta=np.ndarray>

    valid_min :

    0

    valid_max :

    4095

    standard_name :

    rainfall_amount

    long_name :

    SF

    unit :

    mm

    Array Chunk Bytes 6.18 MiB 3.09 MiB Shape (2, 900, 900) (1, 900, 900) Count 8 Tasks 2 Chunks Type float32 numpy.ndarray

• Attributes: (5)

  radarid :

  10000

  formatversion :

  3

  radolanversion :

  2.12.0

  radarlocations :

  ['ham', 'ros', 'emd', 'han', 'bln', 'ess', 'fld', 'drs', 'neu', 'nhb', 'oft', 'eis', 'muc', 'mem']

  radardays :

  ['bln 24', 'drs 24', 'eis 24', 'emd 24', 'ess 24', 'fld 24', 'ham 24', 'han 24', 'mem 24', 'muc 24', 'neu 24', 'nhb 24', 'oft 24', 'ros 24']