xarray RADOLAN backend¶
In this example, we read RADOLAN data files using the xarray radolan
backend.
[1]:
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¶
[2]:
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¶
[3]:
display(comp)
<xarray.Dataset> Dimensions: (time: 1, x: 900, y: 900) Coordinates: * time (time) datetime64[ns] 2014-08-03T09:50:00 * y (y) float64 -4.659e+03 -4.658e+03 ... -3.761e+03 -3.76e+03 * x (x) float64 -523.5 -522.5 -521.5 -520.5 ... 372.5 373.5 374.5 375.5 Data variables: RW (y, x) float32 ... Attributes: radarid: 10000 radolanversion: 2.13.1 radarlocations: ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu',...
xarray.Dataset
- time: 1
- x: 900
- y: 900
- 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.659e+03 -4.658e+03 ... -3.76e+03
- units :
- km
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
array([-4658.644724, -4657.644724, -4656.644724, ..., -3761.644724, -3760.644724, -3759.644724])
- x(x)float64-523.5 -522.5 ... 374.5 375.5
- units :
- km
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
array([-523.462167, -522.462167, -521.462167, ..., 373.537833, 374.537833, 375.537833])
- 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]
- radarid :
- 10000
- radolanversion :
- 2.13.1
- radarlocations :
- ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu', 'nhb', 'oft', 'tur', 'isn', 'fbg', 'mem']
Inspect RADOLAN moments¶
The DataArrays can be accessed by key or by attribute. Each DataArray has dimensions and coordinates of it’s parent dataset.
[5]:
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.659e+03 -4.658e+03 ... -3.761e+03 -3.76e+03 * x (x) float64 -523.5 -522.5 -521.5 -520.5 ... 372.5 373.5 374.5 375.5 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)
- y(y)float64-4.659e+03 -4.658e+03 ... -3.76e+03
- units :
- km
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
array([-4658.644724, -4657.644724, -4656.644724, ..., -3761.644724, -3760.644724, -3759.644724])
- x(x)float64-523.5 -522.5 ... 374.5 375.5
- units :
- km
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
array([-523.462167, -522.462167, -521.462167, ..., 373.537833, 374.537833, 375.537833])
- 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.
[6]:
comp.RW.plot(x="x", y="y", add_labels=False)
[6]:
<matplotlib.collections.QuadMesh at 0x7fc8ba42f790>
[7]:
fig = pl.figure(figsize=(10,10))
ax = fig.add_subplot(111)
comp.RW.plot(x="x", y="y", ax=ax)
[7]:
<matplotlib.collections.QuadMesh at 0x7fc8ba36afd0>
Mask some values¶
[8]:
ds = comp['RW'].where(comp['RW'] >= 1)
ds.plot()
[8]:
<matplotlib.collections.QuadMesh at 0x7fc8ba2500d0>
Export to NetCDF4¶
[9]:
# fix _FillValue
comp.RW.encoding["_FillValue"] = 65535
[10]:
comp.to_netcdf("test_radolan.nc")
Import again¶
[11]:
comp1 = xr.open_dataset("test_radolan.nc")
comp1
[11]:
<xarray.Dataset> Dimensions: (time: 1, x: 900, y: 900) Coordinates: * time (time) datetime64[ns] 2014-08-03T09:50:00 * y (y) float64 -4.659e+03 -4.658e+03 ... -3.761e+03 -3.76e+03 * x (x) float64 -523.5 -522.5 -521.5 -520.5 ... 372.5 373.5 374.5 375.5 Data variables: RW (y, x) float32 ... Attributes: radarid: 10000 radolanversion: 2.13.1 radarlocations: ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu',...
xarray.Dataset
- time: 1
- x: 900
- y: 900
- 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.659e+03 -4.658e+03 ... -3.76e+03
- units :
- km
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
array([-4658.644724, -4657.644724, -4656.644724, ..., -3761.644724, -3760.644724, -3759.644724])
- x(x)float64-523.5 -522.5 ... 374.5 375.5
- units :
- km
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
array([-523.462167, -522.462167, -521.462167, ..., 373.537833, 374.537833, 375.537833])
- 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]
- radarid :
- 10000
- radolanversion :
- 2.13.1
- radarlocations :
- ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu', 'nhb', 'oft', 'tur', 'isn', 'fbg', 'mem']
Check equality¶
[12]:
xr.testing.assert_equal(comp, comp1)
More RADOLAN loading mechanisms¶
Use xr.open_dataset
¶
[13]:
comp2 = xr.open_dataset(f, engine="radolan")
display(comp2)
<xarray.Dataset> Dimensions: (time: 1, x: 900, y: 900) Coordinates: * time (time) datetime64[ns] 2014-08-03T09:50:00 * y (y) float64 -4.659e+03 -4.658e+03 ... -3.761e+03 -3.76e+03 * x (x) float64 -523.5 -522.5 -521.5 -520.5 ... 372.5 373.5 374.5 375.5 Data variables: RW (y, x) float32 ... Attributes: radarid: 10000 radolanversion: 2.13.1 radarlocations: ['boo', 'ros', 'emd', 'hnr', 'pro', 'ess', 'asd', 'neu',...
xarray.Dataset
- time: 1
- x: 900
- y: 900
- 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.659e+03 -4.658e+03 ... -3.76e+03
- units :
- km
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
array([-4658.644724, -4657.644724, -4656.644724, ..., -3761.644724, -3760.644724, -3759.644724])
- x(x)float64-523.5 -522.5 ... 374.5 375.5
- units :
- km
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
array([-523.462167, -522.462167, -521.462167, ..., 373.537833, 374.537833, 375.537833])
- 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]
- radarid :
- 10000
- 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¶
[14]:
#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')
[15]:
comp3 = xr.open_mfdataset(f, engine="radolan")
display(comp3)
<xarray.Dataset> Dimensions: (time: 2, x: 900, y: 900) Coordinates: * time (time) datetime64[ns] 2013-05-27T00:50:00 2013-05-28T00:50:00 * y (y) float64 -4.659e+03 -4.658e+03 ... -3.761e+03 -3.76e+03 * x (x) float64 -523.5 -522.5 -521.5 -520.5 ... 372.5 373.5 374.5 375.5 Data variables: SF (time, y, x) float32 dask.array<chunksize=(1, 900, 900), meta=np.ndarray> Attributes: radarid: 10000 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
- time: 2
- x: 900
- y: 900
- 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.659e+03 -4.658e+03 ... -3.76e+03
- units :
- km
- long_name :
- y coordinate of projection
- standard_name :
- projection_y_coordinate
array([-4658.644724, -4657.644724, -4656.644724, ..., -3761.644724, -3760.644724, -3759.644724])
- x(x)float64-523.5 -522.5 ... 374.5 375.5
- units :
- km
- long_name :
- x coordinate of projection
- standard_name :
- projection_x_coordinate
array([-523.462167, -522.462167, -521.462167, ..., 373.537833, 374.537833, 375.537833])
- SF(time, y, x)float32dask.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
- radarid :
- 10000
- 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']