Projections

Projections#

In this notebook \(\omega radlib\)’s projection capabilities are highlighted, from several helper functions to full featured coordinate reprojection.

import warnings

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pyproj
import wradlib as wrl
import xradar as xd

warnings.filterwarnings("ignore")

Projection Helpers#

get_default_projection#

def_proj = wrl.georef.get_default_projection()
print(type(def_proj))
display(def_proj)

<class 'pyproj.crs.crs.CRS'>

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

get_earth_projection#

Defaults to WGS84 ellipsoid.

earth = wrl.georef.get_earth_projection()
display(earth)

<Geographic 3D CRS: EPSG:4979>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

earth2 = wrl.georef.get_earth_projection(model="geoid")
display(earth2)

<Compound CRS: EPSG:9707>
Name: WGS 84 + EGM96 height
Axis Info [ellipsoidal|vertical]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- H[up]: Gravity-related height (metre)
Area of Use:
- undefined
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Sub CRS:
- WGS 84
- EGM96 height

get_earth_radius#

Return the radius of the earth for a given spheroid model at a given latitude.

\(R^2 = \frac{a^4 \cos(f)^2 + b^4 \sin(f)^2} {a^2 \cos(f)^2 + b^2 \sin(f)^2}\)

lat = np.arange(-90, 90)
er = wrl.georef.get_earth_radius(lat)

plt.plot(lat, er)

[<matplotlib.lines.Line2D at 0x76959d08ead0>]
../../_images/1485496c0bbe925f48c1cac03e52c34d8e984efafb556ea4873828b5cdeb0994.png

get_radar_projection#

Get the native radar projection which is an azimuthal equidistant projection centered at the site using WGS84.

site = (7, 53)
radar = wrl.georef.get_radar_projection(site)
display(radar)

<Projected CRS: PROJCRS["Unknown Azimuthal Equidistant",BASEGEOGCR ...>
Name: Unknown Azimuthal Equidistant
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Azimuthal Equidistant
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

create_crs#

Conveniently create RADOLAN projection objects.

radolan = wrl.georef.create_crs("dwd-radolan")
display(radolan.to_cf())

{'crs_wkt': 'PROJCRS["Radolan Projection",BASEGEOGCRS["Radolan Coordinate System",DATUM["Radolan_Kugel",ELLIPSOID["Erdkugel",6370040,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False easting",0,LENGTHUNIT["kilometre",1000],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["kilometre",1000],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",south,ORDER[1],LENGTHUNIT["kilometre",1000,ID["EPSG",9036]]],AXIS["northing",south,ORDER[2],LENGTHUNIT["kilometre",1000,ID["EPSG",9036]]]]',
 'semi_major_axis': 6370040.0,
 'semi_minor_axis': 6370040.0,
 'inverse_flattening': 0.0,
 'reference_ellipsoid_name': 'Erdkugel',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'Radolan Coordinate System',
 'horizontal_datum_name': 'Radolan_Kugel',
 'projected_crs_name': 'Radolan Projection',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': 10.0,
 'false_easting': 0.0,
 'false_northing': 0.0}

radolan = wrl.georef.create_crs("dwd-radolan-wgs84")
display(radolan.to_cf())

{'crs_wkt': 'PROJCRS["Radolan Projection",BASEGEOGCRS["Radolan Coordinate System",DATUM["Unknown based on WGS 84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.25722356301,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",south,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",south,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.31424518,
 'inverse_flattening': 298.25722356301,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'Radolan Coordinate System',
 'horizontal_datum_name': 'Unknown based on WGS 84 ellipsoid',
 'projected_crs_name': 'Radolan Projection',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': 10.0,
 'false_easting': 0.0,
 'false_northing': 0.0}

radolan = wrl.georef.create_crs("dwd-radolan-wgs84-de1200")
display(radolan.to_cf())

{'crs_wkt': 'PROJCRS["Radolan Projection",BASEGEOGCRS["Radolan Coordinate System",DATUM["Unknown based on WGS 84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.25722356301,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False easting",543196.835217764,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",3622588.861931,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",south,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",south,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.31424518,
 'inverse_flattening': 298.25722356301,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'Radolan Coordinate System',
 'horizontal_datum_name': 'Unknown based on WGS 84 ellipsoid',
 'projected_crs_name': 'Radolan Projection',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': 10.0,
 'false_easting': 543196.835217764,
 'false_northing': 3622588.861931002}

radolan = wrl.georef.create_crs("dwd-radolan-wgs84-rx")
display(radolan.to_cf())

{'crs_wkt': 'PROJCRS["Radolan Projection",BASEGEOGCRS["Radolan Coordinate System",DATUM["Unknown based on WGS 84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.25722356301,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unnamed",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False easting",523196.835217778,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",3772588.86193113,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",south,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",south,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.31424518,
 'inverse_flattening': 298.25722356301,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'Radolan Coordinate System',
 'horizontal_datum_name': 'Unknown based on WGS 84 ellipsoid',
 'projected_crs_name': 'Radolan Projection',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': 10.0,
 'false_easting': 523196.83521777776,
 'false_northing': 3772588.861931134}

ensure_crs#

Helper function to always return correct projection objects. Defaults to pyproj.crs.CRS. For easy conversion from/to pyproj/cartopy/osgeo.osr.

radar_pyproj = wrl.georef.ensure_crs(radar, trg="pyproj")
print(type(radar_pyproj))
display(radar_pyproj)

<class 'pyproj.crs.crs.CRS'>

<Projected CRS: PROJCRS["Unknown Azimuthal Equidistant",BASEGEOGCR ...>
Name: Unknown Azimuthal Equidistant
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Azimuthal Equidistant
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

radar_cartopy = wrl.georef.ensure_crs(radar, trg="cartopy")
print(type(radar_cartopy))
display(radar_cartopy)

<class 'cartopy.crs.CRS'>

<Projected CRS: PROJCRS["Unknown Azimuthal Equidistant",BASEGEOGCR ...>
Name: Unknown Azimuthal Equidistant
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Azimuthal Equidistant
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

radar_osr = wrl.georef.ensure_crs(radar, trg="osr")
print(type(radar_osr))
print(radar_osr)

<class 'osgeo.osr.SpatialReference'>
PROJCS["Unknown Azimuthal Equidistant",
    GEOGCS["unknown",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Azimuthal_Equidistant"],
    PARAMETER["latitude_of_center",53],
    PARAMETER["longitude_of_center",7],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

Reproject#

\(\omega radlib\)’s power work horse for coordinate transformation. Relying on pyproj.transformer.Transformer and pyproj.crs.CRS.

Can be attributed with Coordinate Reference Systems (CRS) of any provenience.

swp = (
    xd.model.create_sweep_dataset(rng=150)
    .swap_dims(time="azimuth")
)
swp = swp.assign_coords(sweep_mode="azimuthal_surveillance")
display(swp)

<xarray.Dataset> Size: 17kB
Dimensions:     (azimuth: 360, range: 1000)
Coordinates:
  * azimuth     (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    time        (azimuth) datetime64[ns] 3kB 2022-08-27T10:00:00 ... 2022-08-...
    elevation   (azimuth) float64 3kB 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
  * range       (range) float64 8kB 75.0 225.0 375.0 ... 1.498e+05 1.499e+05
    longitude   float64 8B 8.788
    latitude    float64 8B 46.17
    altitude    int64 8B 375
    sweep_mode  <U22 88B 'azimuthal_surveillance'
Data variables:
    *empty*
swp = swp.wrl.georef.georeference()
display(swp)

<xarray.Dataset> Size: 17MB
Dimensions:     (azimuth: 360, range: 1000)
Coordinates: (12/15)
  * azimuth     (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    time        (azimuth) datetime64[ns] 3kB 2022-08-27T10:00:00 ... 2022-08-...
    elevation   (azimuth) float64 3kB 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
  * range       (range) float64 8kB 75.0 225.0 375.0 ... 1.498e+05 1.499e+05
    x           (azimuth, range) float64 3MB 0.6544 1.963 ... -1.308e+03
    y           (azimuth, range) float64 3MB 74.98 224.9 ... 1.497e+05 1.498e+05
    ...          ...
    bins        (azimuth, range) float64 3MB 75.0 225.0 ... 1.498e+05 1.499e+05
    longitude   float64 8B 8.788
    latitude    float64 8B 46.17
    altitude    int64 8B 375
    sweep_mode  <U22 88B 'azimuthal_surveillance'
    crs_wkt     int64 8B 0
Data variables:
    *empty*
crs = pyproj.CRS.from_epsg(2056)
swp2 = swp.wrl.georef.reproject(trg_crs=crs)
display(swp2)

<xarray.Dataset> Size: 17MB
Dimensions:     (azimuth: 360, range: 1000)
Coordinates: (12/15)
  * azimuth     (azimuth) float64 3kB 0.5 1.5 2.5 3.5 ... 357.5 358.5 359.5
    time        (azimuth) datetime64[ns] 3kB 2022-08-27T10:00:00 ... 2022-08-...
    elevation   (azimuth) float64 3kB 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
  * range       (range) float64 8kB 75.0 225.0 375.0 ... 1.498e+05 1.499e+05
    x           (azimuth, range) float64 3MB 2.704e+06 2.704e+06 ... 2.7e+06
    y           (azimuth, range) float64 3MB 1.114e+06 1.115e+06 ... 1.264e+06
    ...          ...
    bins        (azimuth, range) float64 3MB 75.0 225.0 ... 1.498e+05 1.499e+05
    longitude   float64 8B 8.788
    latitude    float64 8B 46.17
    altitude    int64 8B 375
    sweep_mode  <U22 88B 'azimuthal_surveillance'
    crs_wkt     int64 8B 0
Data variables:
    *empty*
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

swp.z.plot(x="x", y="y", ax=ax1)
ax1.set_title(swp.xradar.get_crs().name)
swp2.z.plot(x="x", y="y", ax=ax2)
ax2.set_title(swp2.xradar.get_crs().name)
fig.tight_layout()
../../_images/434ac847a45ccf37e8f1f1d49d4afd0a5a6989494e716108d5363f70f325d98e.png