# Reading Radolan RW files

**URL:** <https://openradar.discourse.group/t/reading-radolan-rw-files/39>\
**Category:** Openradar\
**Tags:** wradlib\
**Created:** [October 8, 2022, 5:53pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39 "2022-10-08T17:53:06Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Muru](https://avatars.discourse-cdn.com/v4/letter/m/779978/32.png) [@Muru](https://openradar.discourse.group/u/Muru)\
**Post date:** [October 8, 2022, 5:53pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/1 "2022-10-08T17:53:07Z")

</div>

Hi Guys,

```
     Wradlib Version: 1.16.2
     Python Version: 3.10.5

```

I am able to read and plot the graph with raa01-rw\_10000-2208301350-dwd—bin.gz  
file with xy co-ordinates.

The below code works fine. But I want your help to plot the graph with longitude and latitude. Could you please send me some working code/links for the same product type.

Moreover I need to know if there is any function which can tell me what is the rainfall rate value [data] for a particular longitude and latitude.  
Or how should I have to find it only with data without ploting it?

```auto
if __name__ == ' __main__':
   data, attrs = read_radolan("raa01-rw_10000-2208301350-dwd---bin.gz")
   data = np.ma.masked_equal(data, -9999)
   radolan_grid_xy = wrl.georef.get_radolan_grid(900, 900)
   plot_radolan(data, attrs, radolan_grid_xy, clabel="mm * h-1")
   pl.show()

def read_radolan(rad_file):
     rad_data_file = wrl.util.get_wradlib_data_file(rad_file)
     print_header(rad_data_file)
     return wrl.io.read_radolan_composite(rad_data_file)

def plot_radolan(data, attrs, grid, clabel=None):
   fig = pl.figure(figsize=(10, 8))
   ax = fig.add_subplot(111, aspect="equal")
   x = grid[:, :, 0]
   y = grid[:, :, 1]
   pm = ax.pcolormesh(x, y, data, cmap="viridis")
   cb = fig.colorbar(pm, shrink=0.75)
   cb.set_label(clabel)
   pl.xlabel("x [km]")
   pl.ylabel("y [km]")
   pl.title(
    "{0} Product\n{1}".format(attrs["producttype"], attrs["datetime"].isoformat())
    )
    pl.xlim((x[0, 0], x[-1, -1]))
    pl.ylim((y[0, 0], y[-1, -1]))
    pl.grid(color="r")

```

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [October 8, 2022, 7:32pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/2 "2022-10-08T19:32:01Z")

</div>

Hi Muru,

welcome to our openradar discourse and thanks for the question. I’ll move the discussion to Openradar until we have decided on the final structure.

I’ve created a jupyter notebook gist with a code example for your use case here → [radolan-rw.ipynb](https://gist.github.com/kmuehlbauer/83f5014a43e88dc467d3236b8da1d1c3). One more word of clarification. In the notebook I’ve used wradlib’s RADOLAN xarray backend to conveniently have data and coordinates merged in one Dataset/DataArray. Further reading in the [wradlib backend docs](https://docs.wradlib.org/en/stable/notebooks/fileio/wradlib_radolan_backend.html).

But you would need to be careful when using the DWD grids, they have changed their grid projection base ellipsoid from sphere to wgs84-ellipsoid for their latest products. But as you are using up-to-date wradlib you should have all the bits and pieces available to be prepared.

A write-up on DWD grids is available in this gist → [analyze\_dwd\_grids.ipynb](https://gist.github.com/kmuehlbauer/e050772fd0d5b7f88c71fff6a326f858). Please also refer to the [wradlib docs on RADOLAN](https://docs.wradlib.org/en/stable/notebooks/radolan.html), they are a bit aged and need refreshing, though.

---

<div class="post-metadata">

**Author:** ![Muru](https://avatars.discourse-cdn.com/v4/letter/m/779978/32.png) [@Muru](https://openradar.discourse.group/u/Muru)\
**Post date:** [October 9, 2022, 10:29am UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/3 "2022-10-09T10:29:05Z")

</div>

Hi Kai,  
Thank you for your support. It works very well. It prints like below,

![Figure_lat_long_grid](https://global.discourse-cdn.com/free1/uploads/openradar/original/1X/4073ff3beaae0ee5e8680c8ce68a30f96118ae42.png)

```auto
<xarray.DataArray 'RW' ()>
array(0., dtype=float32)
Coordinates:
    y float64 -3.759e+06
    x float64 37.83
    lon float64 9.993
    lat float64 54.9
Attributes:
    valid_min: 0
    valid_max: 4095
    standard_name: rainfall_rate
    long_name: RW
    unit: mm h-1

```

I would like to understand about the rain-fall rate values. It shows as

```auto
valid_min 0
valid_max 4095
Color values shows as 0 to 60.

```

Could you please help me on  
How can we decide **whether rainfall happened or not** in that location?

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [October 9, 2022, 11:16am UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/4 "2022-10-09T11:16:06Z")

</div>

Hi Muru,

the RADOLAN RW data is taken verbatim from the datafile as `uint16`. The xarray reader is applying CF convention based `scale_factor` and `add_offset` (see at `ds.RW.encoding`) and let xarray do the CF decoding. Please follow-up with the [netcdf-c attribute conventions](https://docs.unidata.ucar.edu/netcdf-c/current/attribute_conventions.html) for `valid_min` and `valid_max`.

So `valid_min` and `valid_max` are the packed values (`uint16`) of the data, which in this case evaluate to 0 and 409.5. If you have a value of `0` than no rain was measured at that location.

---

<div class="post-metadata">

**Author:** ![Muru](https://avatars.discourse-cdn.com/v4/letter/m/779978/32.png) [@Muru](https://openradar.discourse.group/u/Muru)\
**Post date:** [October 9, 2022, 12:48pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/5 "2022-10-09T12:48:11Z")

</div>

Hi Kai,  
Thank you, understood. So rainfall has occurred at this location.

```auto
    lon float64 9.993
    lat float64 54.9

Attributes:
    valid_min: 0
    valid_max: 4095
    standard_name: rainfall_rate
    long_name: RW
    unit: mm h-1

```

But if you look at the plotted map, it does not look like showing rain color there.  
is it possible to take the color value as well for that location?

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [October 9, 2022, 1:11pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/6 "2022-10-09T13:11:28Z")

</div>

Hi Muru,

maybe I was too sparse with my explanations. The `valid_min`/`valid_max` are two dataset-wide values and not point-specific.

The correct value in the above example is `0`:

```auto
<xarray.DataArray 'RW' ()>
array(0., dtype=float32)

```

---

<div class="post-metadata">

**Author:** ![Muru](https://avatars.discourse-cdn.com/v4/letter/m/779978/32.png) [@Muru](https://openradar.discourse.group/u/Muru)\
**Post date:** [October 9, 2022, 1:35pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/7 "2022-10-09T13:35:20Z")

</div>

Hi Kai, Ohh Ok.

```auto
def print_nearest_ds(lat, lon):
    # find nearest xy-grid point to a specific latlon-coordinate
    abslat = np.abs(dsg.lat - lat)
    abslon = np.abs(dsg.lon - lon)
    c = np.maximum(abslon, abslat)
    ([xloc], [yloc]) = np.where(c == np.min(c))

    # Select index location at the x/y dimension
    point_ds = dsg.sel(x=xloc, y=yloc, method="nearest")
    print(point_ds.RW)

print(print_nearest_ds(lat = 47.190504802, lon = 6.510980544))

```

As per the plot, I could see some rain color available in that area. [lat = 47.190504802, lon = 6.510980544]. But still it prints its value as 0 only. Is that function “print\_nearest\_ds” correct? No idea why it is returning nearest as [lon 9.993, lat 54.9].

```auto
<xarray.DataArray 'RW' ()>
array(0., dtype=float32)
Coordinates:
    y float64 -3.759e+06
    x float64 37.83
    lon float64 9.993
    lat float64 54.9
Attributes:
    valid_min: 0
    valid_max: 4095
    standard_name: rainfall_rate
    long_name: RW
    unit: mm h-1

```

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [October 9, 2022, 3:38pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/8 "2022-10-09T15:38:50Z")

</div>

Hi Muru,

sorry that was my bad, I’ve made a mistake when adapting the solution from SO. I’ve fixed the notebook gist (see link above).

For corrected code to determine the wanted pixel see below. There are two changes. First, y/lat is the first dim of our 2d array. Second, we need to use `.isel` to select from the given index. The SO solution worked since in that example x and y didn’t have coordinates and a zero-based index is used with `.sel`.

```python
lat = 50.0
lon = 10.0

# find nearest xy-grid point to a specific latlon-coordinate   
abslat = np.abs(dsg.lat-lat)
abslon = np.abs(dsg.lon-lon)
c = np.maximum(abslon, abslat)

# Attention: y/lat is first dim, get
([yidx], [xidx]) = np.where(c == np.min(c))
    
# Select index location at the x/y dimension
# use isel as we select with index 
point_ds = dsg.isel(x=xidx, y=yidx)

display(point_ds.RW)

```

```auto
<xarray.DataArray 'RW' ()>
array(1.1, dtype=float32)
Coordinates:
    y float64 -4.326e+06
    x float64 37.83
    lon float64 9.994
    lat float64 50.0
Attributes:
    valid_min: 0
    valid_max: 4095
    standard_name: rainfall_rate
    long_name: RW
    unit: mm h-1    

```

As you can see, the lon/lat coordinates are now giving back our selected values (“nearest”).

---

<div class="post-metadata">

**Author:** ![Muru](https://avatars.discourse-cdn.com/v4/letter/m/779978/32.png) [@Muru](https://openradar.discourse.group/u/Muru)\
**Post date:** [October 9, 2022, 4:29pm UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/9 "2022-10-09T16:29:34Z")

</div>

Hi Kai,  
Thank you so much. It works like gem. I could now see the proper rainfall value [14.7].

```auto
<xarray.DataArray 'RW' ()>
array(14.7, dtype=float32)
Coordinates:
    y float64 -4.65e+06
    x float64 -2.83e+05
    lon float64 6.512
    lat float64 47.19
Attributes:
    valid_min: 0
    valid_max: 4095
    standard_name: rainfall_rate
    long_name: RW
    unit: mm h-1

```

---

<div class="post-metadata">

**Author:** ![kmuehlbauer](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/kmuehlbauer/32/439_2.png) [@kmuehlbauer](https://openradar.discourse.group/u/kmuehlbauer)\
**Post date:** [October 10, 2022, 11:29am UTC](https://openradar.discourse.group/t/reading-radolan-rw-files/39/10 "2022-10-10T11:29:07Z")

</div>

A post was split to a new topic: [How to plot the graph to a specific location within a certain radius](https://openradar.discourse.group/t/how-to-plot-the-graph-to-a-specific-location-within-a-certain-radius/42)