Incorrect coordinates when transforming from RADOLAN to WGS84 using pyproj

General
1

Hello everyone,

I’m trying to process DWD’s RADOLAN-RW data and display precipitation data on a map. Unfortunately, I’m getting incorrect geographic coordinates after the transformation and need help troubleshooting.

My goal:

  • Download DWD RADOLAN-RW data (HDF5 format)

  • Extract precipitation values

  • Transform coordinates from RADOLAN projection to WGS84

  • Display heatmap on OpenStreetMap

My code:

I’m using this PROJ4 definition for the RADOLAN projection:

radolan_proj4 = (
    "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 "
    "+x_0=543196.83521776402 +y_0=3622588.8619310018 "
    "+a=6378137 +b=6356752.3142451802 +units=m +no_defs"
)

And I’m transforming the coordinates as follows:

radolan_crs = CRS.from_proj4(radolan_proj4)
wgs84 = CRS.from_epsg(4326)
transformer = Transformer.from_crs(radolan_crs, wgs84, always_xy=True)

# For each pixel in the raster:
x_rad = col * DX  # DX = 1000 (1km grid)
y_rad = (nrows - 1 - row) * DY  # DY = 1000, invert Y axis
lon, lat = transformer.transform(x_rad, y_rad)

The problem:

The transformed coordinates don’t seem to be correct. The resulting points are not positioned over Germany but at incorrect locations (mostly far to the north).

What I’ve tried already:

  1. Different PROJ4 definitions for RADOLAN

  2. Various transformation approaches (always_xy=True/False)

  3. Alternative calculations of RADOLAN coordinates

  4. Comparison with a script that uses GK3 transformation (which works correctly)

My question:

Can anyone help me understand why the transformation with the RADOLAN stereographic projection isn’t working correctly? Is the issue with:

  • The PROJ4 definition?

  • The calculation of RADOLAN coordinates?

  • The transformation itself?

  • Something else?

Is there an authoritative source for the correct PROJ4 definition of the RADOLAN projection?

here my code:

import h5py
import numpy as np
from pyproj import CRS, Transformer
import json
import requests
from tqdm import tqdm
import logging
import argparse
import folium
from folium.plugins import HeatMap
import os
import pandas as pd

# === Konfiguration ===
URL = "https://opendata.dwd.de/weather/radar/radolan/rw/raa01-rw_10000-latest-dwd---bin.hdf5"
LOCAL_FILE = "raa01-rw_latest.hdf5"
OUTPUT_JSON = "rain_rasters.json"
OUTPUT_HTML = "rain_map.html"

# Rasterparameter
NROWS = 1200
NCOLS = 1100
DX = 1000   # 1 km Rasterweite
DY = 1000

# === CRS/Transformator ===
radolan_proj4 = (
    "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 "
    "+x_0=543196.83521776402 +y_0=3622588.8619310018 "
    "+a=6378137 +b=6356752.3142451802 +units=m +no_defs"
)
radolan_crs = CRS.from_proj4(radolan_proj4)
wgs84 = CRS.from_epsg(4326)
transformer = Transformer.from_crs(radolan_crs, wgs84, always_xy=True)

# === CLI Argumente ===
parser = argparse.ArgumentParser(description="DWD RADOLAN RW für OpenHAB")
parser.add_argument("-v", "--verbose", action="store_true", help="Logs anzeigen")
args = parser.parse_args()

log_level = logging.INFO if args.verbose else logging.WARNING
logging.basicConfig(level=log_level, format="%(asctime)s [%(levelname)s] %(message)s")

# === Download ===
def download_latest_file():
    logging.info("Starte Download...")
    r = requests.get(URL, stream=True, timeout=30)
    r.raise_for_status()
    total_size = int(r.headers.get("content-length", 0))
    block_size = 1024
    if args.verbose:
        t = tqdm(total=total_size, unit="B", unit_scale=True, desc="Download")
    with open(LOCAL_FILE, "wb") as f:
        for data in r.iter_content(block_size):
            f.write(data)
            if args.verbose:
                t.update(len(data))
    if args.verbose:
        t.close()
    logging.info(f"Download abgeschlossen: {LOCAL_FILE}")

# === RADOLAN Bounding Box berechnen ===
def calculate_radolan_bbox():
    """Berechnet die exakten Eckpunkte des RADOLAN-Gebiets in WGS84"""
    # RADOLAN-Koordinaten der Eckpunkte
    corners_radolan = [
        (0, 0),                    # Südwest
        (NCOLS * DX, 0),           # Südost
        (NCOLS * DX, NROWS * DY),  # Nordost
        (0, NROWS * DY)            # Nordwest
    ]
    
    corners_wgs84 = []
    for x, y in corners_radolan:
        lon, lat = transformer.transform(x, y)
        corners_wgs84.append((lat, lon))  # Folium erwartet (lat, lon)
    
    return corners_wgs84

# === Daten extrahieren (RADOLAN-Raster) ===
def extract_rain_coordinates(path, rain_threshold_mm=0.1):
    rain_points = []

    with h5py.File(path, "r") as f:
        ds = f["/dataset1/data1"]
        what = ds["what"].attrs

        gain = float(what.get("gain", 0.1))
        offset = float(what.get("offset", 0.0))
        nodata = float(what.get("nodata", 65535.0))

        raw = ds["data"][:].astype(float)
        raw[raw == nodata] = np.nan
        data = raw * gain + offset

        nrows, ncols = data.shape
        rain_mask = data > rain_threshold_mm
        logging.info(f"Niederschlagspunkte >{rain_threshold_mm}mm: {np.sum(rain_mask)}")

        for row in range(nrows):
            for col in range(ncols):
                if data[row, col] > rain_threshold_mm:
                    # RADOLAN-Koordinaten
                    x_rad = col * DX
                    y_rad = (nrows - 1 - row) * DY  # Y invertieren
                    lon, lat = transformer.transform(x_rad, y_rad)

                    rain_points.append({
                        "lon": float(lon),
                        "lat": float(lat),
                        "rain_mm": float(data[row, col])
                    })

    return rain_points

# === Karte erstellen ===
def create_rain_map(rain_points, output_html):
    if not rain_points:
        logging.warning("Keine Regenpunkte für Karte")
        return

    heatmap_data = [[p["lat"], p["lon"], p["rain_mm"]] for p in rain_points]

    # RADOLAN-Bounding Box berechnen
    radolan_bbox = calculate_radolan_bbox()
    
    # Kartenmittelpunkt (Deutschland)
    center_lat = 51.1657
    center_lon = 10.4515
    
    m = folium.Map(
        location=[center_lat, center_lon],
        zoom_start=6,
        tiles='OpenStreetMap',
        width='100%',
        height='100%'
    )

    # Heatmap hinzufügen
    HeatMap(
        heatmap_data,
        name='Niederschlag',
        min_opacity=0.4,
        max_opacity=0.8,
        radius=12,
        blur=8,
        gradient={0.1: 'blue', 0.3: 'cyan', 0.5: 'lime', 0.7: 'yellow', 0.9: 'red'}
    ).add_to(m)

    # RADOLAN-Bounding Box als Overlay einzeichnen
    folium.Polygon(
        locations=radolan_bbox,
        color="red",
        weight=2,
        fill=False,
        tooltip="RADOLAN Abdeckung",
        opacity=0.8
    ).add_to(m)

    # === No-Cache + Auto-Refresh ergänzen ===
    meta_html = """
    <head>
      <meta http-equiv="Cache-Control" content="no-cache, no-store, must-revalidate"/>
      <meta http-equiv="Pragma" content="no-cache"/>
      <meta http-equiv="Expires" content="0"/>
      <meta http-equiv="refresh" content="300"/> <!-- alle 5 Minuten -->
    </head>
    """
    m.get_root().html.add_child(folium.Element(meta_html))

    title_html = f'''
    <div style="position: fixed; top: 10px; left: 50px; z-index: 1000; background: white; padding: 10px; border-radius: 5px;">
        <h4 style="margin: 0;">🌧 DWD RADOLAN Niederschlag</h4>
        <p style="margin: 5px 0 0 0; font-size: 12px;">Aktualisiert: {pd.Timestamp.now().strftime('%Y-%m-%d %H:%M')}</p>
    </div>
    '''
    m.get_root().html.add_child(folium.Element(title_html))

    # Karte auf Deutschland zentrieren
    m.fit_bounds([[47.0, 5.0], [55.0, 15.0]])

    m.save(output_html)
    os.chmod(output_html, 0o644)
    logging.info(f"Karte gespeichert: {output_html}")

# === Hauptprogramm ===
if __name__ == "__main__":
    try:
        download_latest_file()
        rain_coords = extract_rain_coordinates(LOCAL_FILE, rain_threshold_mm=0.1)

        if rain_coords:
            with open(OUTPUT_JSON, "w") as f:
                json.dump(rain_coords, f, indent=2)
            print(f"✅ JSON gespeichert: {OUTPUT_JSON}")

            create_rain_map(rain_coords, OUTPUT_HTML)
            print(f"🗺️  Karte gespeichert: {OUTPUT_HTML}")

            example = rain_coords[0]
            print(f"📋 Beispiel-Punkt: lon={example['lon']:.6f}, lat={example['lat']:.6f}, rain_mm={example['rain_mm']:.1f}")
        else:
            print("☀ Kein Niederschlag erkannt.")

    except Exception as e:
        logging.exception("Fehler während des Skriptlaufs.")
        print(f"❌ Fehler: {e}")

Thank you for any help!