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?
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. 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.
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.
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 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.
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].
<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
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.
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)
<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”).
