Problems with georeferencing

[Jetstream]

Hello everyone,
I have problems with georeferencing with the HX composite of the German Weather Service. When I convert the data to WGS84, the individual pixels are slightly distorted (see images). Sorry - as a new user I can’t upload more than 1x picture…

Is there a way to avoid or correct these ‘distortions’? Or are they negligible for analyses?

I actually want to plot an exactly symmetrical radar pixel that lies exactly at the correct lat/long coordinate. I tried this once, but then the lat coordinates in particular don’t fit at all (image3 / image4) -

I hope you can help me and thank you in advance. The forum is really great and the websites with all the ideas are really useful.

Hier ist der Code zu Bild 3:

import matplotlib
matplotlib.use("TkAgg")

import h5py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from pyproj import Transformer
import wradlib as wrl

# -----------------------------------------------------
# 1. RADOLAN HX-Datei laden
# -----------------------------------------------------
file = "/Pfad/zur/Datei/composite_hx_20250331_1225-hd5"
with h5py.File(file, 'r') as f:
    dbz_raw = f["dataset1/data1/data"][:]
    gain = f["dataset1/data1/what"].attrs["gain"]
    offset = f["dataset1/data1/what"].attrs["offset"]

dbz = dbz_raw * gain + offset
dbz = np.ma.masked_outside(dbz, 10, 70)

# -----------------------------------------------------
# 2. Farbskala: Gradient mit viridis (oder turbo)
# -----------------------------------------------------
cmap = plt.get_cmap("viridis")  # alternativ: "turbo"
norm = Normalize(vmin=10, vmax=70)

# -----------------------------------------------------
# 3. Original-Koordinaten (DE4800, stereographisch)
# -----------------------------------------------------
nrows, ncols = dbz.shape
radolan_xy = wrl.georef.get_radolan_grid(nrows, ncols, center=True, proj="dwd-radolan")
x = radolan_xy[..., 0]
y = radolan_xy[..., 1]

# -----------------------------------------------------
# 4. Reprojektion nach EPSG:25832 (UTM Zone 32N)
# -----------------------------------------------------
proj_radolan = "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 +a=6370040 +b=6370040 +no_defs"
transformer = Transformer.from_crs(proj_radolan, "epsg:25832", always_xy=True)
x_utm, y_utm = transformer.transform(x, y)

# Flip Radarbild passend zur Koordinatenachse (Nord oben)
dbz = np.flipud(dbz)

# -----------------------------------------------------
# 5. Plot mit schlanker Gradient-Legende
# -----------------------------------------------------
fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.UTM(zone=32))
ax.set_title("RADOLAN HX (EPSG:25832) mit Gradient-Legende")

ax.set_extent([x_utm.min(), x_utm.max(), y_utm.min(), y_utm.max()], crs=ccrs.UTM(32))

# Basisfeatures (optional)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
ax.add_feature(cfeature.STATES, linewidth=0.4, edgecolor='gray')

# Gitterlinien in Lat/Lon
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=0.5, color='gray', alpha=0.5)
gl.top_labels = False
gl.right_labels = False

# Radarbild zeichnen
pm = ax.pcolormesh(x_utm, y_utm, dbz, transform=ccrs.UTM(32),
                   cmap=cmap, norm=norm, shading="auto")

# -----------------------------------------
# Gradient-Farbleiste (wie im Screenshot)
# -----------------------------------------
cb = plt.colorbar(pm, ax=ax, orientation="vertical", pad=0.01, shrink=0.75)
cb.set_label("Radarreflektivität (dBZ)", labelpad=15)
cb.set_ticks(np.arange(10, 75, 5))  # Ticks alle 5 dBZ

plt.tight_layout()
plt.show()

[kmuehlbauer]

Hi @Jetstream,

It’s not really clear to me what you are trying to achieve. But some of your assumptions don’t seem to be correct.

The DWD HX Composite Format isn’t really well described as of now (https://www.dwd.de/DE/leistungen/radarprodukte/radarkomposit_hx_hdf5.pdf). But we can at least get something meaningful out of it:

with h5py.File(file, 'r') as f:
    dbz_raw = f["dataset1/data1/data"][:]
    gain = f["dataset1/data1/what"].attrs["gain"]
    offset = f["dataset1/data1/what"].attrs["offset"]
    where_attrs = dict(f["where"].attrs.items())

In the global where-group we find all necessary projection information to build our needed coordinate arrays and CRS:

print(where_attrs)
{
    'LL_lat': np.float64(45.696425377390064),
    'LL_lon': np.float64(3.5669946350078914),
    'LR_lat': np.float64(45.68460578137082),
    'LR_lon': np.float64(16.580869348598274),
    'UL_lat': np.float64(55.862087108249824),
    'UL_lon': np.float64(1.463301510256666),
    'UR_lat': np.float64(55.845438563255755),
    'UR_lon': np.float64(18.73161645466747),
    'projdef': np.bytes_(b'+proj=stere +lat_ts=60 +lat_0=90 +lon_0=10 +x_0=543571.83521776402       +y_0=3622213.8619310022 +units=m +a=6378137 +b=6356752.3142451802 +no_defs'),
    'xscale': np.float64(250.0),
    'xsize': np.uint64(4400),
    'yscale': np.float64(250.0),
    'ysize': np.uint64(4800)
}

Let’s take a look at the projdef:

print(where_attrs["projdef"].decode("utf-8")
'
 +proj=stere +lat_ts=60 +lat_0=90 +lon_0=10 
 +x_0=543571.83521776402 +y_0=3622213.8619310022 
 +units=m 
 +a=6378137 +b=6356752.3142451802 
 +no_defs
'

The values of +a/+b indicate a WGS84 ellipsoid (used in new DWD products, instead of the DWD Erdkugel). Now we can create the x, y coordinate arrays. Note, that there is one inconsistency in the given information, yscale needs to be negative, to create it proper:

xscale = where_attrs["xscale"]
yscale = -where_attrs["yscale"]
xsize = where_attrs["xsize"]
ysize = where_attrs["ysize"]
x = np.arange(xsize) * xscale
y = np.arange(ysize) * yscale
x, y = np.meshgrid(x, y)
projdef = where_attrs["projdef"].decode("utf8")

The UTM coordinates can be transformed according to the given projdef:

transformer = Transformer.from_crs(projdef, "epsg:25832", always_xy=True)
x_utm, y_utm = transformer.transform(x, y)

The nice thing is we do not have to flip the source data, since everything is already aligned. The complete script, including download of the most recent HX composit:

import matplotlib

import h5py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from pyproj import Transformer
import wradlib as wrl

# -----------------------------------------------------
# 1. RADOLAN HX-Datei laden
# -----------------------------------------------------
# only works in IPython/Jupyter
!wget -N https://opendata.dwd.de/weather/radar/composite/hx/composite_hx_LATEST-hd5

file = "composite_hx_LATEST-hd5"

with h5py.File(file, 'r') as f:
    dbz_raw = f["dataset1/data1/data"][:]
    gain = f["dataset1/data1/what"].attrs["gain"]
    offset = f["dataset1/data1/what"].attrs["offset"]
    where_attrs = dict(f["where"].attrs.items())

dbz = dbz_raw * gain + offset
# just for better representation, uncomment this as needed
# dbz = np.ma.masked_outside(dbz, 10, 70)

# -----------------------------------------------------
# 2. Farbskala: Gradient mit viridis (oder turbo)
# -----------------------------------------------------
cmap = plt.get_cmap("viridis")  # alternativ: "turbo"
norm = Normalize(vmin=10, vmax=70)

# -----------------------------------------------------
# 3. Original-Koordinaten (DE4800, stereographisch)
# -----------------------------------------------------
xscale = where_attrs["xscale"]
yscale = -where_attrs["yscale"]
xsize = where_attrs["xsize"]
ysize = where_attrs["ysize"]
x = np.arange(xsize) * xscale
y = np.arange(ysize) * yscale
x, y = np.meshgrid(x, y)
projdef = where_attrs["projdef"].decode("utf8")

# -----------------------------------------------------
# 4. Reprojektion nach EPSG:25832 (UTM Zone 32N)
# -----------------------------------------------------
transformer = Transformer.from_crs(projdef, "epsg:25832", always_xy=True)
x_utm, y_utm = transformer.transform(x, y)

# -----------------------------------------------------
# 5. Plot mit schlanker Gradient-Legende
# -----------------------------------------------------
fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.UTM(zone=32))
ax.set_title("RADOLAN HX (EPSG:25832) mit Gradient-Legende")

ax.set_extent([x_utm.min(), x_utm.max(), y_utm.min(), y_utm.max()], crs=ccrs.UTM(32))

# Basisfeatures (optional)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
ax.add_feature(cfeature.STATES, linewidth=0.4, edgecolor='gray')

# Gitterlinien in Lat/Lon
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color='gray', alpha=0.5)
gl.top_labels = False
gl.right_labels = False

# Radarbild zeichnen
pm = ax.pcolormesh(x_utm, y_utm, dbz, cmap=cmap, norm=norm, shading="auto")

# -----------------------------------------
# Gradient-Farbleiste (wie im Screenshot)
# -----------------------------------------
cb = plt.colorbar(pm, ax=ax, orientation="vertical", pad=0.01, shrink=0.75)
cb.set_label("Radarreflektivität (dBZ)", labelpad=15)
cb.set_ticks(np.arange(10, 75, 5))  # Ticks alle 5 dBZ

plt.tight_layout()
plt.show()

Same functionality is also possible using wradlib radolan grid functionality, but you would need to use "dwd-radolan-wgs84" as projection. I did not check, if the grid-reference is correct (lower-left vs center).

HTH,
Kai

P.S. As DWD is in the process to move to HDF5 for their RADOLAN products, wradlib will surely make a move too and add those to it’s composite reader.

[Jetstream]

Hi Kai,
Thank you for your solution. I am currently testing several radar products from the German Weather Service. Especially the composites (HX / DMAX). Unfortunately, as you say, the information base at the German Weather Service is thin.

The aim is to filter out severe thunderstorms with high radar intensity dbz and reconstruct train paths. The next step is to take a look at the probability of hail.
In my first attempt, I noticed a distortion of the radar pixels, which can probably lead to a deviation of 20 to 30 metres (but these deviations are probably negligible).
I had solved this here:

# GUI-Backend setzen (für interaktive Anzeige im Skript)
import matplotlib
matplotlib.use("TkAgg")

import h5py
import numpy as np
import matplotlib.pyplot as plt
import wradlib as wrl
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# -----------------------------------------------------
# 1. Datei einlesen (Pfad ggf. anpassen)
# -----------------------------------------------------
file = "/Pfad/zur/Datei/composite_hx_20250331_1225-hd5"
with h5py.File(file, 'r') as f:
    dbz_raw = f["dataset1/data1/data"][:]
    attrs = f["dataset1/data1/what"].attrs
    gain = attrs["gain"]
    offset = attrs["offset"]
    nodata = attrs["nodata"]

# Debugging: Ausgabe der Attribute
print(f"gain: {gain}")
print(f"offset: {offset}")
print(f"nodata: {nodata}")

# -----------------------------------------------------
# 2. dBZ skalieren und maskieren
# -----------------------------------------------------
dbz_scaled = dbz_raw * gain + offset
dbz_final = np.ma.masked_outside(dbz_scaled, 10, 70)

print(f"dBZ min: {dbz_final.min()}")
print(f"dBZ max: {dbz_final.max()}")
print(f"Gültige Werte: {dbz_final.count()}")

# -----------------------------------------------------
# 3. RADOLAN-Gitterkoordinaten berechnen (WGS84)
# -----------------------------------------------------
nrows, ncols = dbz_final.shape
radolan_grid_ll = wrl.georef.get_radolan_grid(nrows, ncols, wgs84=True)
lon = radolan_grid_ll[:, :, 0]
lat = radolan_grid_ll[:, :, 1]

# -----------------------------------------------------
# 4. Flip (Spiegelung Y-Achse korrigieren)
# -----------------------------------------------------
dbz_final_flipped = np.flipud(dbz_final)

# -----------------------------------------------------
# 5. Plotfunktion mit kontinuierlicher Farbskala (Gradient)
# -----------------------------------------------------
def create_radar_fig_cartopy(lon, lat, data):
    fig = plt.figure(figsize=(10, 10))
    ax = plt.axes(projection=ccrs.Mercator())

    # Kartenausschnitt: ganz Deutschland
    ax.set_extent([4.5, 16.5, 46.5, 55.5], crs=ccrs.PlateCarree())

    # Kartenfeatures hinzufügen
    ax.add_feature(cfeature.COASTLINE, linewidth=0.6)
    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
    ax.add_feature(cfeature.STATES, linewidth=0.4, edgecolor='gray')

    # Kontinuierliche Colormap verwenden
    cmap = plt.get_cmap("viridis")  # oder z. B. "turbo", "plasma", "magma"
    pm = ax.pcolormesh(lon, lat, data, transform=ccrs.PlateCarree(),
                       cmap=cmap, shading="auto", vmin=10, vmax=70)

    # Gradient-Farblegende
    cb = plt.colorbar(pm, ax=ax, orientation="vertical", pad=0.03, shrink=0.75)
    cb.set_label("Radarreflektivität (dBZ)")

    ax.set_title("RADOLAN HX Radar-Komposit (Gradient)", fontsize=14)
    plt.tight_layout()
    plt.show(block=True)

# -----------------------------------------------------
# 6. Plot aufrufen
# -----------------------------------------------------
create_radar_fig_cartopy(lon, lat, dbz_final_flipped)

After all, I wanted to make an attempt to obtain a symmetrical shape of the radar pixels. But I failed because of the projection. Unfortunately, the information available from the German Weather Service is also sparse.
Your solution definitely works. Thank you so much :-).

Unfortunately, the projections drive me crazy. If I want to use an Openstreet map, the pixels are in the wrong place again… : The OPM cards in the following image are incorrect:

Do you have any ideas? That would help me a lot.

That was my attempt:

import warnings
import numpy as np
import h5py
from matplotlib import pyplot as plt
from matplotlib.colors import Normalize
from pyproj import Transformer
from PIL import Image
import folium
import os
import webbrowser

# --------------------------------------------
# 1. Pfade & Datei
# --------------------------------------------
output_folder = "/Users/Downloads/"
png_filename = "radar_turbo_overlay.png"
html_filename = "radar_turbo_map.html"

output_png = os.path.join(output_folder, png_filename)
output_html = os.path.join(output_folder, html_filename)
file = os.path.join(output_folder, "composite_hx_20250331_1225-hd5")

warnings.filterwarnings("ignore", category=RuntimeWarning)

# --------------------------------------------
# 2. Daten & Attribute laden
# --------------------------------------------
with h5py.File(file, 'r') as f:
    dbz_raw = f["dataset1/data1/data"][:]
    gain = f["dataset1/data1/what"].attrs["gain"]
    offset = f["dataset1/data1/what"].attrs["offset"]
    where_attrs = dict(f["where"].attrs.items())
    projdef = where_attrs["projdef"].decode("utf8")

# --------------------------------------------
# 3. dBZ skalieren und maskieren
# --------------------------------------------
dbz = dbz_raw * gain + offset
dbz = np.ma.masked_outside(dbz, 10, 70)

# --------------------------------------------
# 4. Farbskala: Turbo-Gradiert (wie Cartopy-Plot)
# --------------------------------------------
cmap = plt.get_cmap("turbo")
norm = Normalize(vmin=10, vmax=70)

# --------------------------------------------
# 5. Zellkantenraster berechnen (metrisch)
# --------------------------------------------
xscale = where_attrs["xscale"]
yscale = -where_attrs["yscale"]
xsize = where_attrs["xsize"]
ysize = where_attrs["ysize"]
x_edges = (np.arange(xsize + 1) - 0.5) * xscale
y_edges = (np.arange(ysize + 1) - 0.5) * yscale
x_grid, y_grid = np.meshgrid(x_edges, y_edges)

# --------------------------------------------
# 6. Zellkanten reprojizieren → EPSG:4326
# --------------------------------------------
transformer = Transformer.from_crs(projdef, "EPSG:4326", always_xy=True)
lon_edges, lat_edges = transformer.transform(x_grid, y_grid)
lat_min, lat_max = lat_edges.min(), lat_edges.max()
lon_min, lon_max = lon_edges.min(), lon_edges.max()
bounds = [[lat_min, lon_min], [lat_max, lon_max]]

# --------------------------------------------
# 7. PNG erzeugen mit Padding (4801 x 4401)
# --------------------------------------------
dbz_padded = np.pad(dbz.filled(0), ((0, 1), (0, 1)), mode="edge")
rgba = cmap(norm(dbz_padded))
rgba = (rgba * 255).astype(np.uint8)
Image.fromarray(rgba).save(output_png)

# --------------------------------------------
# 8. Folium-Karte erzeugen
# --------------------------------------------
m = folium.Map(location=[51, 10], zoom_start=6, tiles="OpenStreetMap")
folium.raster_layers.ImageOverlay(
    name="RADOLAN HX (Turbo-Legende)",
    image=output_png,
    bounds=bounds,
    opacity=0.7,
    interactive=True,
    cross_origin=False,
    zindex=1
).add_to(m)

folium.LayerControl().add_to(m)
m.save(output_html)

# --------------------------------------------
# 9. Öffnen im Browser
# --------------------------------------------
if os.path.exists(output_html):
    print(f"✅ HTML-Karte gespeichert unter:\n{output_html}")
    webbrowser.open(f"file://{output_html}")

if os.path.exists(output_png):
    print(f"✅ PNG-Overlay gespeichert unter:\n{output_png}")

Here is the test file if anyone wants to try / copy it.

composite_hx_20250331_1225-hd5

Have a nice weekend and thank you for your support!