Radolan bin data export as Geotif and reproject to EPSG:3857 (WGS 84)

[till Kadabra]

Hi people,

i try my best to follow the documentation on wradlib docs but i get in trouble when i try to reporject a geotiff file created from RADOLAN SF data.

import wradlib as wrl
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# We will export this RADOLAN dataset to a GIS compatible format
data_raw, meta  = wrl.io.read_radolan_composite("/home/elisabeth/Dokumente/master/2014_07/raa01-sf_10000-1407282350-dwd---bin")

# This is the RADOLAN projection
proj_osr = wrl.georef.create_osr("dwd-radolan")

# Get projected RADOLAN coordinates for corner definition
xy_raw = wrl.georef.get_radolan_grid(900, 900)

data, xy = wrl.georef.set_raster_origin(data_raw, xy_raw, 'upper')

# create 3 bands
wdir = "/home/elisabeth/Dokumente/master/2014_01/out/"
data = np.stack((data, data+100, data+1000))
ds = wrl.georef.create_raster_dataset(data, xy, projection=proj_osr)
wrl.io.write_raster_dataset(wdir + "R_2014_07_28.tif", ds, 'GTiff')

I am right that the exported Geotiffs Coordinate system is the radolan projection?

Now i need to reporject (is it the right thing?) this to EPSG:3857
So i read the tif file into my script...

ds1 = wrl.io.open_raster(wdir + "R_2014_07_28.tif")
data1, xy1, proj1 = wrl.georef.extract_raster_dataset(ds1, nodata=-9999.)
np.testing.assert_array_equal(data1, data)
np.testing.assert_array_equal(xy1[:-1,:-1,:], xy)

but now i am stuck. I dont know how to reproject the coordinate system and export it again as geotiff. May somebody can help me ;)?

Best regards
Till

[till Kadabra]

I guess i have to take this function https://docs.wradlib.org/en/stable/generated/wradlib.georef.raster.reproject_raster_dataset.html

new_raster = wradlib.georef.raster.reproject_raster_dataset(ds1, **kwargs)

but i don't understand the Keyword Arguments which i have to choose to reproject  to WGS84 ;()

[Kai Muehlbauer]

Hi,

yes that would be the correct function.

What you have to keep in mind: While working with raster data we need to
interpolate between different projections (like UTM, RADOLAN polar
stereographic, Gauss-Krueger and the like) if we want to have raster
data as in- and output . If we just want to project coordinates to
another projection we will gain 2D-coordinates instead of 1D (cartesian).

If your dataset has already the correct projection attached (like in
your example), you can specify the target coordinate system using the
keyword `projection_target` which needs to be a OSR coordinate reference
object.

A working example using your code snippets would be:

import wradlib as wrl
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# We will export this RADOLAN dataset to a GIS compatible format
data_raw, meta =
wrl.io.read_radolan_composite("/home/elisabeth/Dokumente/master/2014_07/raa01-sf_10000-1407282350-dwd---bin")

# This is the RADOLAN projection
proj_osr = wrl.georef.create_osr("dwd-radolan")

# Get projected RADOLAN coordinates for corner definition
xy_raw = wrl.georef.get_radolan_grid(900, 900)

data, xy = wrl.georef.set_raster_origin(data_raw, xy_raw, 'upper')

# create 3 bands
wdir = "/home/elisabeth/Dokumente/master/2014_01/out/"
data = np.stack((data, data+100, data+1000))
ds = wrl.georef.create_raster_dataset(data, xy, projection=proj_osr)

# EPSG 3857
epsg_3857 = wrl.georef.epsg_to_osr(3857)
ds2 = wrl.georef.reproject_raster_dataset(ds,
projection_target=epsg_3857, spacing=1000)
# This is a bug I just discovered, we need to put the next line in any
case until fixed
ds2.SetProjection(epsg_3857.ExportToWkt())
wrl.io.write_raster_dataset(wdir + "R_2014_07_28.tif", ds2, 'GTiff')

There is some inconsistency in the wradlib function which will be fixed
with the next release. Follow up in
https://github.com/wradlib/wradlib/issues/368.

Cheers,
Kai



Am 26.07.19 um 21:14 schrieb till Kadabra:> I guess i have to take this

> --
> 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/88d3ab00-cbbb-472e-9707-3311955d68ab%40googlegroups.com
>
<https://groups.google.com/d/msgid/wradlib-users/88d3ab00-cbbb-472e-9707-3311955d68ab%40googlegroups.com?utm_medium=email&utm_source=footer>.


Am 26.07.19 um 20:53 schrieb till Kadabra:
> |
> importwradlib aswrl
> importnumpy asnp
> importwarnings

> warnings.filterwarnings('ignore')
>
> # We will export this RADOLAN dataset to a GIS compatible format
> data_raw,meta

>  =wrl.io.read_radolan_composite("/home/elisabeth/Dokumente/master/2014_07/raa01-sf_10000-1407282350-dwd---bin")

>
> # This is the RADOLAN projection

> proj_osr =wrl.georef.create_osr("dwd-radolan")

>
> # Get projected RADOLAN coordinates for corner definition
> xy_raw =wrl.georef.get_radolan_grid(900,900)
>

> data,xy =wrl.georef.set_raster_origin(data_raw,xy_raw,'upper')

[till Kadabra]

Thanks, the code works great. With the new updated package i tested the code without the line:

ds2.SetProjection(epsg_3857.ExportToWkt())

and it still works.

[till Kadabra]

Today i am plot the data and i got this ;))

Do you know what the hell is going on here? The higher the dates of the image the higher image is located on the map. What exactly means spacing ? i thought it is the length between two cells?

[till Kadabra]

I solve the problem by shifting the definition and creation of osr objects and coordinates inside the for loop. I don't know why but know every image has the same geometry

Umfang