Scu format into a utm projection (EPSG:25832)

Python
1

I have this file format ‘hd2107132355.scu’ from Radolan and I would like to convert it into a geotiff file.
So far, I could just read the file (see the output below) and extract the numpy array. Is it any possibility that I can convert this file into a raster file using wradlib?

h5_gen = wrl.io.read_generic_hdf5(h5_filename)

output:
{‘dataset_DXk/image’: {‘attrs’: {‘CLASS’: b’IMAGE’,
‘IMAGE_BACKGROUNDINDEX’: 0,
‘IMAGE_COLORMODEL’: b’RGB’,
‘IMAGE_MINMAXRANGE’: array([0, 0], dtype=uint8),
‘IMAGE_SUBCLASS’: b’BITMAP’,
‘IMAGE_VERSION’: b’1.2’,
‘PALETTE’: },
‘data’: array([[0., 0., 0., …, 0., 0., 0.],
[0., 0., 0., …, 0., 0., 0.],
[0., 0., 0., …, 0., 0., 0.],
…,
[0., 0., 0., …, 0., 0., 0.],
[0., 0., 0., …, 0., 0., 0.],
[0., 0., 0., …, 0., 0., 0.]], dtype=float32)},
‘dataset_DXk/what’: {‘attrs’: {‘corr_code’: b’INTT’,
‘corr_no’: 1,
‘corr_perc’: array([0], dtype=int32),
‘enddate [YYYYMMDD]’: b’20210713’,
‘endtime [HHmmss]’: b’235500’,
‘gain’: 0.1,
‘nodata’: -999.0,
‘ori_format’: b’RY’,
‘ori_name’: b’HD2107132355.scu’,
‘product’: b’COMP ‘,
‘quantity’: b’RATE ‘,
‘rad’: 427,
‘undetect’: -999.0}},
‘dataset_DXk/where’: {‘attrs’: {‘xsize’: 900, ‘ysize’: 900}},
‘how’: {‘attrs’: {‘WMO’: 10000,
‘adjustment’: 0.99099994,
‘maxlevel [dBZ]’: 9.69,
‘minlevel [dBZ]’: 0.01,
‘nodes’: b’asb,boo,ros,hnr,umd,pro,ess,fld,drs,neu,nhb,oft,eis,tur,isn,fbg,mem’,
‘nodesn’: 17,
‘place’: b’dwd’,
‘wavelength [cm]’: 5.0,
‘zr_a’: 256.0,
‘zr_b’: 1.42}},
‘what’: {‘attrs’: {‘date [YYYYMMDD]’: b’20210713’,
‘object’: b’COMP ‘,
‘sets’: 1,
‘time [HHmmss]’: b’235500’,
‘version [H5Vol 0.2]’: b’SCP 5.2.4.13 ATDZ '}},
‘where’: {‘attrs’: {‘LL_lat [deg]’: 46.9526,
‘LL_lon [deg]’: 3.5889,
‘UL_ipixel’: 3072,
‘UL_jpixel’: 6068,
‘UR_lat [deg]’: 54.7405,
‘UR_lon [deg]’: 15.7208,
‘lat [deg]’: 51.0,
‘lon [deg]’: 9.0,
‘start_lat [deg]’: 6068.895,
‘start_lon [deg]’: 3072.5002,
‘xscale [m]’: 950.0,
‘xsize’: 900,
‘yscale [m]’: 950.0,
‘ysize’: 900}}}

2

Hi @brunoparente93, welcome to openradar discourse!

I’ve thought that I’ve seen almost all shades of RADOLAN, but this scu-format is obviously some handcrafted thing.

If it’s true what’s given as metadata (“RY”), you might want to look at RADKLIM (Download RADKLIM products). In RADKLIM-project the RY data has been reprocessed as YW in a consistent manner.

If you use original RADOLAN data provided by DWD, you’d be able to work along Export a dataset in GIS-compatible format — wradlib to convert to any raster GDAL can think of.

If you need to use this scu-data you can still use the above example, but would need to create the xarray Dataset by yourself from the source data.

HTH,

Kai

3

Thank you for the reply, I tried your example and I could generate a readable shapefile for QGIS (code 1). But what I could not figure out is how I can export the raster in EPSG:25832 already. I tried to convert the raster generated by code 1 using code 2, and it worked, but the radolan raster is always being displayed in the middle of the ocean, totally way from the correct position. Can you please give me a hint about what I am doing wrong? Thank you in advance.

Code 1

import matplotlib.pyplot as plt

import os
import wradlib as wrl
import xarray as xr
import numpy as np
import warnings
from pyproj.crs import CRS
warnings.filterwarnings("ignore")
try:
get_ipython().run_line_magic("matplotlib inline")
except:
plt.ion()

# We will export this RADOLAN dataset to a GIS compatible format

os.environ["WRADLIB_DATA"] = r"C:\Users\thomas\Desktop\python_codes\wradlib_radar"

wdir = wrl.util.get_wradlib_data_path() + "/radolan/grid/"

if not os.path.exists(wdir):

os.makedirs(wdir)

filename = "raa01-sf_10000-2311272050-dwd---bin"

filename = wrl.util.get_wradlib_data_file(filename)

ds = xr.open_dataset(filename, engine="radolan")

display(ds)

# This is the RADOLAN projection

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

crs = CRS.from_wkt(proj_osr.ExportToWkt(["FORMAT=WKT2_2018"]))

print(proj_osr)

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

data, xy = wrl.georef.set_raster_origin(ds.SF.values, xy_raw, "upper")
print(data.shape)

# create 3 bands
#data = np.stack((data, data + 100, data + 1000), axis=0)
print(data.shape)
gds = wrl.georef.create_raster_dataset(data, xy, crs=proj_osr)
wrl.io.write_raster_dataset(wdir + "geotiff3.tif", gds, driver="GTiff")

Code 2

from osgeo import gdal, osr

# Path to your raster file
input_raster_path = r'C:\Users\thomas\Desktop\python_codes\wradlib_radar\radolan\grid\geotiff3.tif'

# Output raster file path with the new coordinate system
output_raster_path = r'C:\Users\thomas\Desktop\python_codes\wradlib_radar\radolan\grid\file_epsg25832.tif'

# Target EPSG code (EPSG:25832 for UTM Zone 32N, you can replace it with your target code)
target_epsg_code = 25832

# Open the input dataset
input_dataset = gdal.Open(input_raster_path)

# Check if the dataset was successfully opened
if input_dataset is not None:
    # Create a SpatialReference object for the target coordinate system
    target_srs = osr.SpatialReference()
    target_srs.ImportFromEPSG(target_epsg_code)

    # Set the target coordinate system to the input dataset
    input_dataset.SetProjection(target_srs.ExportToWkt())

    # Create an output dataset with the new coordinate system
    driver = gdal.GetDriverByName('GTiff')  # Choose the appropriate driver for your output format
    output_dataset = driver.CreateCopy(output_raster_path, input_dataset)

    # Close the datasets
    input_dataset = None
    output_dataset = None

    print(f"Conversion to EPSG:{target_epsg_code} complete. Output saved to {output_raster_path}")
else:
    print("Failed to open the input dataset.")
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

5

The script worked very well here. Thank you so much, Kai. I will think of possible ways to combine your script with this scu format. One thing that I thought of would be to replace the array from your xarray with the array from the SCU format file, and then generate the raster. It worked here, because the same dimensions, I just have to certify if they have the same coordinates, which I think they do, because they came from the same radolan product.

I wish you all the best! Best regards!
Bruno

6

Hi Bruno,

glad it worked. And yes you might just swap the arrays as we are only fiddling with the coordinates and do not interpolate grids. Just make sure, that the metadata is correct (data units etc.).

Best,
Kai

P.S. You might mark my above post as solution, Then this will be shown more prominent. Or you post your final workflow here and mark that. Up to you :slight_smile: