I am using the python library wradlib to read radolan radar datasets provided by the opendata service of the Deutscher Wetterdienst (DWD). I want to crop the data by providing maximum and minimum longitude and latitude coordinates so that just a selected area remains which I can perform further operations on.
My code looks like this:
import wradlib as wrl
rv_file = "DE1200_RV2201271935_000"
ds = wrl.io.open_radolan_dataset(rv_file)
data = ds.RV
print(data)
The output shows that the array only has x and y coordinates:
<xarray.DataArray 'RV' (y: 1200, x: 1100)>
[1320000 values with dtype=float32]
Coordinates:
* y (y) float64 -4.809e+03 -4.808e+03 ... -3.611e+03 -3.61e+03
* x (x) float64 -543.5 -542.5 -541.5 -540.5 ... 552.5 553.5 554.5 555.5
Attributes:
valid_min: 0
valid_max: 4095
standard_name: rainfall_rate
long_name: RV
unit: mm h-1
I know that I need to use a reprojection and add a geographic coordinate system to the array. I tried using the following functions:
proj_stereo = wrl.georef.create_osr("dwd-radolan")
proj_wgs = osr.SpatialReference()
proj_wgs.ImportFromEPSG(4326)
radolan_grid_xy = wrl.georef.get_radolan_grid(1200, 1100)
But I am not sure how to apply the reprojection and then crop the data.
The projections for the radar data and the target reference system need to be defined. Then the lon-lat-grid can be created by taking the radolan grid and the source and target projections.
proj_stereo = wrl.georef.create_osr("dwd-radolan")
proj_wgs = osr.SpatialReference()
proj_wgs.ImportFromEPSG(4326)
radolan_grid_x_rv, radolan_grid_y_rv = 1200, 1100
radolan_grid_xy = wrl.georef.get_radolan_grid(radolan_grid_x_rv, radolan_grid_y_rv)
radolan_grid_ll = wrl.georef.reproject(radolan_grid_xy, projection_source=proj_stereo, projection_target=proj_wgs)
The new coordinates can be assigned to the data as a dictionary. The data can be cropped by filtering by lon and lat coordinates.
lonlat_coords = dict(lon=(["y", "x"], radolan_grid_ll[:, :, 0]), lat=(["y", "x"], radolan_grid_ll[:, :, 1]))
new_data = data.assign_coords(lonlat_coords)
lon_min, lon_max = 8.2, 8.9
lat_min, lat_max = 53.2, 53.8
selected_area = new_data.where((new_data.lon > lon_min) & (new_data.lon < lon_max) & (new_data.lat > lat_min) & (new_data.lat < lat_max), drop=True)
print(selected_area)