Converting DWD RADOLAN data results in smaller tile size

[Nicolas T.]

Hello everyone,

my name is Nicolas and I am currently using the DWD Radolan format to use the precipitation forecast data (WN). Wradlib is used to convert the binary data which has worked great so far.

Right now I am struggling to get the tile borders of a coordinate. 

To describe my problem a little more detailed I will provide a small walkthrough of my code:

1. Create a GeoTIFF file of the converted WN data:

f = wrl.util.get_wradlib_data_file(path)
data_raw, meta = wrl.io.read_radolan_composite(f)
proj_osr = wrl.georef.create_osr("dwd-radolan")
xy_raw = wrl.georef.get_radolan_grid(1200, 1100) # see the format description https://www.dwd.de/DE/leistungen/radarprodukte/formatbeschreibung_wndaten.pdf?__blob=publicationFile&v=6
data, xy = wrl.georef.set_raster_origin(data_raw, xy_raw, 'upper')
# change no data value (-9999) to 0
data[data < -1000] = 0
# transform to dbz
data_dbz = np.where(data > 0, data / 2 - 32.5, data)
# calculate indizes for x y data of tile relative to grid

# this probably the most interesting part for my question
indizes = np.indices(data_dbz.shape)
indizes = np.fliplr(indizes)
indizes = np.flipud(indizes)

# now i have the correct order for the relative coordinates
data = np.stack((data_dbz, indizes[0], indizes[1]))

ds = wrl.georef.create_raster_dataset(data, xy, projection=proj_osr)
wrl.io.write_raster_dataset(path + '/' + tifname, ds, 'GTiff')

After running this method I have a GeoTIFF file which is used later on and my data consists of three bands (1: precipitation value, 2: x grid coordinate, 3: y grid coordinate)

Now I use another tool (PyQGIS) to retrieve the band values for a specific coordinate.

The next step would be to get the tile borders of the specific coordinates. For that I used the wradlib documentation for grid reprojection (https://docs.wradlib.org/en/stable/notebooks/radolan/radolan_grid.html#Grid-Reprojection) and modified it for my use case

2. Get grid border

radolan_grid_ll = wrl.georef.get_radolan_grid(1200,1100, wgs84=True) 

y0 = 788 #example values of grid coordinates
x0 = 775
# print(radolan_grid_ll.shape)
print("ll: {:f} {:f} ".format(radolan_grid_ll[x0][y0][0], radolan_grid_ll[x0][y0][1]))
print("lr: {:f} {:10.4f} ".format(radolan_grid_ll[x0][y0+1][0], radolan_grid_ll[x0][y0+1][1]))
print("ur: {:10.4f} {:10.4f} ".format(radolan_grid_ll[x0+1][y0+1][0], radolan_grid_ll[x0+1][y0+1][1]))
print("ul: {:10.4f} {:10.4f} ".format(radolan_grid_ll[x0+1][y0][0], radolan_grid_ll[x0+1][y0][1]))

Now to my problem. In most of the cases the coordinate is not located in the rectangle I create from the four border coordinates. For some coordinates I would have to increase the x or y tile value by one while for other coordinates its fine. 

I struggled a long time to find an explanation for this issue and I realized that the created tiles do not have a resolution of 1 km x 1km as described in the format specification (https://www.dwd.de/DE/leistungen/radarprodukte/formatbeschreibung_wndaten.pdf?__blob=publicationFile&v=6). The resolution of the tiles are approximately 960 m x 960 m and that is the reason why the coordinates aren`t located in the border rectangle sometimes.

My question now is whether it is an issue converting the 1200 x 1100 km grid as I just read about five grid sizes in (https://docs.wradlib.org/en/stable/notebooks/radolan/radolan_grid.html#Standard-Formats) or if the problem is coming from the DWD and the resolution is not exactly 1 km x 1 km.

I hope my question is understandable, if not please say that and I will try to be more precise, and someone can help me.

Thanks in advance and kind regards

Nicolas

[Kai Mühlbauer]

Hi Nicolas,

Thanks for the extensive roundup.

From first glance, I would pinpoint the problem in bullet point 2.

As you retrieve the RADOLAN coordinates in wgs84 they do not represent a cartesian grid anymore. So your tiles won't be squares, but some kind of distorted trapezoid.

There might be more involved, but I need to play with it first.

Best,
Kai

--
Diese Nachricht wurde von meinem Android-Gerät mit K-9 Mail gesendet.

[Nicolas T.]

Hi Kai,

thanks for your fast reply and your first guess.

I think you might be right that the retrieval of the tile border coordinates is different to the result you get when using wrl.georef.get_radolan_grid without WGS84 set to true.

When I run this for a specific tile I get output that looks like this:

       x [km]      y [km]
ll: 244.537833  -4025.644724
lr: 245.537833  -4025.6447
ur:   245.5378  -4024.6447
ul:   244.5378  -4024.6447

There the borders seem to be calculated correctly.

However I looked at the output TIFF file generated at the end of step 1:

ds = wrl.georef.create_raster_dataset(data, xy, projection=proj_osr)
wrl.io.write_raster_dataset(path + '/' + tifname, ds, 'GTiff')

When importing this file as a raster layer into QGIS I still receive the  tile size of approximately 960 m x 960 m.

The same thing happens when exporting the data as a ESRI ASCII file instead of GTiff (https://docs.wradlib.org/en/stable/notebooks/fileio/wradlib_gis_export_example.html) or using a reprojected WGS84 TIFF file.

So all the different outputs result in a tile size of 960 m x 960 m. 

If you need the tiff files or anything else, please let me know. 

Thanks and kind regards

Nicolas