Reprojection of RADOLAN RY binary files to UTM

[Tim Schneider]

Hello Everybody,

first of all thanks a lot for this great library! I have a question concerning the conversion of DWD RADOLAN RY binary data to a GeoTIF in UTM (EPSG 25832). I managed to write a script that does the conversion. As a plausibility check, i loaded the same bin file using the radolan2map plugin from QGIS. The results from wradlib and from radolan2map plugin do not overlap, they differ in extent with the radolan2map tif covering a wider area. Without reprojection, both files seem to overlap fine.

What am i missing here?

Here is a link to a folder containing my script, the binary file, the two resulting TIF files and a QGIS project:

https://drive.google.com/drive/folders/1ObgljDBgQmnGYncvsf6HyGngNAyZtAXX?usp=sharing

Best regards

Tim

And here is my code:

----------------

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
convert RADOLAN ry data from binary polar stereographic to GeoTif UTM files
and clip to an extent given by shapefile
@author: tim
"""

import os
import wradlib as wrl
from osgeo import osr
import numpy as np
os.environ["WRADLIB_DATA"] = "D:/RADOLAN_Reprojection_Example"

# radolan ry grid
xy_stereo = wrl.georef.get_radolan_grid(900,900)

# radolan projection
proj_stereo = wrl.georef.create_osr("dwd-radolan")

# utm projection
proj_utm = osr.SpatialReference()
proj_utm.ImportFromEPSG(25832)

# radolan ry grid reprojected to EPSG 25832
xy_utm = wrl.georef.reproject(xy_stereo,
                              projection_source=proj_stereo,
                              projection_target=proj_utm,
                              spacing=1000)

binfile = 'raa01-ry_10000-2203232330-dwd---bin'

# open and read binary file
f = wrl.util.get_wradlib_data_file(binfile)
data_raw, attrs = wrl.io.read_radolan_composite(f, missing=np.nan)

# process raster origin. I read that upper is used for many datasets, however
# it seems to me that lower is the right choice for ry datasets?
data, xy_utm = wrl.georef.set_raster_origin(data_raw, xy_utm, 'lower')

# write outfile
outf = 'result_wradlib.tif'
ds = wrl.georef.create_raster_dataset(data, xy_utm, projection=proj_utm)
wrl.io.write_raster_dataset(outf, ds, 'GTiff')

[Kai Muehlbauer]

Hello Tim,

you are missing just one step to resample the curved grid data onto
cartesian grid. See below:

import os
import wradlib as wrl
from osgeo import osr
import numpy as np

# radolan ry grid
xy_stereo = wrl.georef.get_radolan_grid(900,900)

# radolan projection
proj_stereo = wrl.georef.create_osr("dwd-radolan")

# utm projection
proj_utm = osr.SpatialReference()
proj_utm.ImportFromEPSG(25832)

# open and read binary file

binfile = 'raa01-ry_10000-2203232330-dwd---bin'
data_raw, attrs = wrl.io.read_radolan_composite(binfile, missing=np.nan)

# process raster origin, move to origin upper, as tif is using it
data, xy_stereo = wrl.georef.set_raster_origin(data_raw, xy_stereo, 'upper')

#create raster dataset in original cartesian coordinates
ds = wrl.georef.create_raster_dataset(data, xy_stereo,
projection=proj_stereo)

# resample/reproject to UTM
ds_out = wrl.georef.reproject_raster_dataset(ds, spacing=1000,
resample=0, projection_target=proj_utm)

# write outfile
outf = 'result_wradlib.tif'

wrl.io.write_raster_dataset(outf, ds_out, 'GTiff')

There are still some differences between the wradlib approach and the
radolan2map-approach. Although both are using GDAL under the hood, the
outputs still differ a bit. It might be due to both using different
approaches (wradlib - gdal.ReprojectImage, radolan2map - gdal.Warp) but
it could also be due to using different parameterizations. Good chance
that there might be slight differences in other parts of the machinery.

If someone finds an inaccurracy within wradlib, please report here or
open an issue on github.

Best,
Kai



Am 28.03.22 um 12:26 schrieb Tim Schneider:

> --
> You received this message because you are subscribed to the Google
> Groups "wradlib-users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to wradlib-user...@googlegroups.com
> <mailto:wradlib-user...@googlegroups.com>.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/wradlib-users/2e7fe1cf-f156-4cb2-b207-7427b50c04fdn%40googlegroups.com
> <https://groups.google.com/d/msgid/wradlib-users/2e7fe1cf-f156-4cb2-b207-7427b50c04fdn%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
Kai Muehlbauer
Institute of Geosciences, Meteorology Section, University of Bonn
Auf dem Huegel 20 | +49 228 739083
D-53121 Bonn | kai.mue...@uni-bonn.de

[Schneider, Tim]

Hello Kai,

thank you, a lot, for the quick help and the elaboration on the wradlib and the radolan2map approach, it works well now!

Best regards
Tim

-----Ursprüngliche Nachricht-----
Von: wradli...@googlegroups.com <wradli...@googlegroups.com> Im Auftrag von Kai Muehlbauer
Gesendet: Montag, 28. März 2022 15:32
An: wradli...@googlegroups.com
Betreff: Re: Reprojection of RADOLAN RY binary files to UTM

To unsubscribe from this group and stop receiving emails from it, send an email to wradlib-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/wradlib-users/3b2178eb-7cb0-f771-3f1a-ad6bc856824c%40uni-bonn.de.