Scu format into a utm projection (EPSG:25832)

Python
4

@brunoparente93

You have the coordinates in RADOLAN polar stereographic projection but write UTM Zone 32N projection information. You would need to reproject the coordinates to UTM Zone 32N before.

But there might be a even more convenient solution to your task (wradlib >=2.0).

Please note, for the below code you’ll need also rioxarray for the raster output.

import wradlib as wrl
import xarray as xr
import rioxarray
import pyproj
import xradar as xd

# load dataset
filename = "raa01-sf_10000-2311271150-dwd---bin"
ds = xr.open_dataset(filename, engine="radolan")

# get radolan crs (conforms to coordinates in ds) and add to dataset
# wradlib should really attach the crs info to the dataset on load ;-)
if ds.attrs.get("formatversion", 3) >= 5:
    radolan = "dwd-radolan-wgs84"
else:
    radolan = "dwd-radolan-sphere"
proj_stereo = wrl.georef.create_osr(radolan)
crs = pyproj.CRS.from_wkt(proj_stereo.ExportToWkt(["FORMAT=WKT2_2018"]))
ds = xd.georeference.add_crs(ds, crs=crs)

# this reprojection step is normally not needed
# the GIS system will make the transformation anyway
# reproject dataset to wanted output projection
# this also changes crs_wkt to dataset
proj_out = wrl.georef.epsg_to_osr(25832)
ds = ds.wrl.georef.reproject(trg_crs=proj_out, coords=dict(x="x", y="y"),
)

# rename crs_wkt to spatial_ref to work with rasterio
ds = ds.rename(crs_wkt="spatial_ref")

# only one _FillValue, use last given
ds.SF.encoding["_FillValue"] = ds.SF.encoding["_FillValue"][-1]

# export to GeoTiff using rioxarray/rasterio
ds.SF.rio.to_raster(filename+".tiff", driver="COG")

At least this fit’s perfectly within QGIS.

Qgis
Qgis1536×931 214 KB

HTH,
Kai