RADOLAN Grid#
Polar Stereographic Projection#
The projected composite raster is equidistant with a grid-spacing of 1.0 km in most cases. There are composites which have 2.0 km grid-spacing (e.g. PC).
There are five different grid sizes, the well-known 900 rows by 900 columns (normal, 1km grid spacinf), 1500 rows by 1400 columns (extended, 1km grid spacing), 460 rows by 460 columns (small, 2km grid spacing) and the legacy 450 rows by 450 rows (2km grid spacing). Since the RADSYS-E project is finalized an extended national composite with 1100 rows by 900 columns (normal_wx, 1km grid spacing) is available, too.
Common to all is that the plane of projection intersects the earth sphere at \(\phi_0 = 60.0^{\circ}N\). The cartesian co-ordinate system is aligned parallel to the \(\lambda_0 = 10.0^{\circ}E\) meridian.
The reference point (\(\lambda_m\), \(\phi_m\)) is \(9.0^{\circ}E\) and \(51.0^{\circ}N\), which is the center of the two smaller grids. The extended european grid has an offset in respect to this reference point of 350km by 150km, the extended national grid 100km by -80km.
The earth as sphere with an radius of 6370.04 km is used for all calculations.
With formulas (1), (2) and (3) the geographic reference points (\(\lambda\), \(\phi\)) can be converted to projected cartesian coordinates. The calculated (x y) is the distance vector to the origign of the cartesian coordinate system (north pole).
\(\begin{equation} x = R * M(\phi) * cos(\phi) * sin(\lambda - \lambda_0) \tag{1} \end{equation}\)
\(\begin{equation} y = -R * M(\phi) * cos(\phi) * cos(\lambda - \lambda_0) \tag{2} \end{equation}\)
\(\begin{equation} M(\phi) = \frac {1 + sin(\phi_0)} {1 + sin(\phi)} \tag{3} \end{equation}\)
Assumed the point (\(10.0^{\circ}E\), \(90.0^{\circ}N\)) is defined as coordinate system origin. Then all ccordinates can be calculated with the known grid-spacing d as:
\(\begin{equation} x = x_0 + d * (j - j_0) \tag{4} \end{equation}\)
\(\begin{equation} y = y_0 + d * (i - i_0) \tag{5} \end{equation}\)
with i, j as cartesian indices.
\(\omega radlib\) provides the convenience function wradlib.georef.get_radolan_grid() which returns the radolan grid for further processing. It takes nrows and ncols as parameters and returns the projected cartesian coordinates or the wgs84 coordinates (keyword arg wgs84=True) as numpy ndarray (nrows x ncols x 2).
[1]:
import wradlib as wrl
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
try:
get_ipython().run_line_magic("matplotlib inline")
except:
plt.ion()
import numpy as np
900x900, 1km grid box#
[2]:
radolan_grid_xy = wrl.georef.get_radolan_grid(900, 900)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_xy.shape, *radolan_grid_xy[0, 0, :])
)
radolan_grid_ll = wrl.georef.get_radolan_grid(900, 900, wgs84=True)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_ll.shape, *radolan_grid_ll[0, 0, :])
)
(900, 900, 2), (-523.4622, -4658.6447)
(900, 900, 2), (3.5889, 46.9526)
1100x900, 1km grid box#
[3]:
radolan_grid_xy = wrl.georef.get_radolan_grid(1100, 900)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_xy.shape, *radolan_grid_xy[0, 0, :])
)
radolan_grid_ll = wrl.georef.get_radolan_grid(1100, 900, wgs84=True)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_ll.shape, *radolan_grid_ll[0, 0, :])
)
(1100, 900, 2), (-443.4622, -4758.6447)
(1100, 900, 2), (4.6759, 46.1929)
1500x1400, 1km grid box#
[4]:
radolan_grid_xy = wrl.georef.get_radolan_grid(1500, 1400)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_xy.shape, *radolan_grid_xy[0, 0, :])
)
radolan_grid_ll = wrl.georef.get_radolan_grid(1500, 1400, wgs84=True)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_ll.shape, *radolan_grid_ll[0, 0, :])
)
(1500, 1400, 2), (-673.4622, -5008.6447)
(1500, 1400, 2), (2.3419, 43.9336)
460x460, 2km grid box#
[5]:
radolan_grid_xy = wrl.georef.get_radolan_grid(460, 460)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_xy.shape, *radolan_grid_xy[0, 0, :])
)
radolan_grid_ll = wrl.georef.get_radolan_grid(460, 460, wgs84=True)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_ll.shape, *radolan_grid_ll[0, 0, :])
)
(460, 460, 2), (-533.4622, -4668.6447)
(460, 460, 2), (3.4814, 46.8603)
450x450, 2km grid box#
[6]:
radolan_grid_xy = wrl.georef.get_radolan_grid(450, 450)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_xy.shape, *radolan_grid_xy[0, 0, :])
)
radolan_grid_ll = wrl.georef.get_radolan_grid(450, 450, wgs84=True)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_ll.shape, *radolan_grid_ll[0, 0, :])
)
(450, 450, 2), (-523.4622, -4658.6447)
(450, 450, 2), (3.5889, 46.9526)
Inverse Polar Stereographic Projection#
The geographic coordinates of specific datapoints can be calculated by using the cartesian coordinates (x,y) and the following formulas:
\(\begin{equation} \lambda = \arctan\left(\frac {-x} {y}\right) + \lambda_0 \tag{6} \end{equation}\)
\(\begin{equation} \phi = \arcsin\left(\frac {R^2 * \left(1 + \sin\phi_0\right)^2 - \left(x^2 + y^2\right)} {R^2 * \left(1 + \sin\phi_0\right)^2 + \left(x^2 + y^2\right)}\right) \tag{7} \end{equation}\)
Standard Formats#
WKT-String#
The German Weather Service provides a WKT-string. This WKT (well known text) is used to create the osr-object representation of the radolan projection.
For the scale_factor the intersection of the projection plane with the earth sphere at \(60.0^{\circ}N\) has to be taken into account:
\(\begin{equation} scale\_factor = \frac {1 + \sin\left(60.^{\circ}\right)} {1 + \sin\left(90.^{\circ}\right)} = 0.93301270189 \tag{8} \end{equation}\)
Also, the PROJECTION["Stereographic_North_Pole"] isn’t known within GDAL/OSR. It has to be changed to the known PROJECTION["polar_stereographic"].
With these adaptions we finally yield the Radolan Projection as WKT-string. This WKT-string is used within \(\omega radlib\) to create the osr-object by using the helper-function wradlib.georef.create_osr().
[7]:
proj_stereo = wrl.georef.create_osr("dwd-radolan")
print(proj_stereo)
PROJCS["Radolan Projection",
GEOGCS["Radolan Coordinate System",
DATUM["Radolan_Kugel",
SPHEROID["Erdkugel",6370040,0]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",60],
PARAMETER["central_meridian",10],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["kilometre",1000,
AUTHORITY["EPSG","9036"]],
AXIS["Easting",SOUTH],
AXIS["Northing",SOUTH]]
PROJ#
Using the above WKT-String the PROJ representation can be derived as:
+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 +k=0.93301270189
+x_0=0 +y_0=0 +a=6370040 +b=6370040 +to_meter=1000 +no_defs
This PROJ-string can also be used to create the osr-object by using the helper-function wradlib.georef.projstr_to_osr():
[8]:
# create radolan projection osr object
dwd_string = (
"+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 +k=0.93301270189 "
"+x_0=0 +y_0=0 +a=6370040 +b=6370040 +to_meter=1000 +no_defs"
)
proj_stereo = wrl.georef.projstr_to_osr(dwd_string)
print(proj_stereo)
PROJCS["unknown",
GEOGCS["unknown",
DATUM["unknown",
SPHEROID["unknown",6370040,0]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",90],
PARAMETER["central_meridian",10],
PARAMETER["scale_factor",0.93301270189],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["kilometre",1000,
AUTHORITY["EPSG","9036"]],
AXIS["Easting",SOUTH],
AXIS["Northing",SOUTH]]
Grid Reprojection#
Within \(\omega radlib\) the wradlib.georef.reproject() function can be used to convert the radolan grid data from xy-space to lonlat-space and back. First, we need to create the necessary Spatial Reference Objects for the RADOLAN-projection and wgs84.
[9]:
from osgeo import osr
proj_stereo = wrl.georef.create_osr("dwd-radolan")
print(proj_stereo)
proj_wgs = osr.SpatialReference()
proj_wgs.ImportFromEPSG(4326)
print(proj_wgs)
PROJCS["Radolan Projection",
GEOGCS["Radolan Coordinate System",
DATUM["Radolan_Kugel",
SPHEROID["Erdkugel",6370040,0]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",60],
PARAMETER["central_meridian",10],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["kilometre",1000,
AUTHORITY["EPSG","9036"]],
AXIS["Easting",SOUTH],
AXIS["Northing",SOUTH]]
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AXIS["Latitude",NORTH],
AXIS["Longitude",EAST],
AUTHORITY["EPSG","4326"]]
Then, we call reproject with the osr-objects as projection_source and projection_target parameters.
[10]:
radolan_grid_xy = wrl.georef.get_radolan_grid(900, 900)
radolan_grid_ll = wrl.georef.reproject(
radolan_grid_xy, src_crs=proj_stereo, trg_crs=proj_wgs
)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_ll.shape, *radolan_grid_ll[0, 0, :])
)
(900, 900, 2), (3.5889, 46.9526)
And the other way round.
[11]:
radolan_grid_xy = wrl.georef.reproject(
radolan_grid_ll, src_crs=proj_wgs, trg_crs=proj_stereo
)
print(
"{0}, ({1:.4f}, {2:.4f})".format(radolan_grid_xy.shape, *radolan_grid_xy[0, 0, :])
)
(900, 900, 2), (-523.4622, -4658.6447)
In the following example the RADOLAN grid is projected to wgs84 and GaussKrüger Zone3.
[12]:
# create UTM zone 32 projection osr object
proj_utm32 = osr.SpatialReference()
proj_utm32.ImportFromEPSG(32632)
# transform radolan polar stereographic projection to wgs84 and then to utm zone 32
radolan_grid_ll = wrl.georef.reproject(
radolan_grid_xy, src_crs=proj_stereo, trg_crs=proj_wgs
)
radolan_grid_utm32 = wrl.georef.reproject(
radolan_grid_ll, src_crs=proj_wgs, trg_crs=proj_utm32
)
lon_wgs0 = radolan_grid_ll[:, :, 0]
lat_wgs0 = radolan_grid_ll[:, :, 1]
x_utm32 = radolan_grid_utm32[:, :, 0]
y_utm32 = radolan_grid_utm32[:, :, 1]
x_rad = radolan_grid_xy[:, :, 0]
y_rad = radolan_grid_xy[:, :, 1]
print("\n------------------------------")
print("source radolan x,y-coordinates")
print(" {0} {1} ".format("x [km]", "y [km]"))
print("ll: {:10.4f} {:10.3f} ".format(x_rad[0, 0], y_rad[0, 0]))
print("lr: {:10.4f} {:10.3f} ".format(x_rad[0, -1], y_rad[0, -1]))
print("ur: {:10.4f} {:10.3f} ".format(x_rad[-1, -1], y_rad[-1, -1]))
print("ul: {:10.4f} {:10.3f} ".format(x_rad[-1, 0], y_rad[-1, 0]))
print("\n--------------------------------------")
print("transformed radolan lonlat-coordinates")
print(" {0} {1} ".format("lon [degE]", "lat [degN]"))
print("ll: {:10.4f} {:10.4f} ".format(lon_wgs0[0, 0], lat_wgs0[0, 0]))
print("lr: {:10.4f} {:10.4f} ".format(lon_wgs0[0, -1], lat_wgs0[0, -1]))
print("ur: {:10.4f} {:10.4f} ".format(lon_wgs0[-1, -1], lat_wgs0[-1, -1]))
print("ul: {:10.4f} {:10.4f} ".format(lon_wgs0[-1, 0], lat_wgs0[-1, 0]))
print("\n-----------------------------------")
print("transformed radolan utm32-coordinates")
print(" {0} {1} ".format("easting [m]", "northing [m]"))
print("ll: {:10.0f} {:10.0f} ".format(x_utm32[0, 0], y_utm32[0, 0]))
print("lr: {:10.0f} {:10.0f} ".format(x_utm32[0, -1], y_utm32[0, -1]))
print("ur: {:10.0f} {:10.0f} ".format(x_utm32[-1, -1], y_utm32[-1, -1]))
print("ul: {:10.0f} {:10.0f} ".format(x_utm32[-1, 0], y_utm32[-1, 0]))
------------------------------
source radolan x,y-coordinates
x [km] y [km]
ll: -523.4622 -4658.645
lr: 375.5378 -4658.645
ur: 375.5378 -3759.645
ul: -523.4622 -3759.645
--------------------------------------
transformed radolan lonlat-coordinates
lon [degE] lat [degN]
ll: 3.5889 46.9526
lr: 14.6087 47.0711
ur: 15.7042 54.7327
ul: 2.0736 54.5790
-----------------------------------
transformed radolan utm32-coordinates
easting [m] northing [m]
ll: 88298 5214122
lr: 925789 5228352
ur: 931372 6085693
ul: 52659 6070030