Georeferencing

Georeferencing#

GeorefMethods

wradlib xarray SubAccessor methods for Georef.

Miscellaneous#

`bin_altitude`	Calculates the height of a radar bin taking the refractivity of the atmosphere into account.
`bin_distance`	Calculates great circle distance from radar site to radar bin over spherical earth, taking the refractivity of the atmosphere into account.
`site_distance`	Calculates great circle distance from bin at certain altitude to the radar site over spherical earth, taking the refractivity of the atmosphere into account.
`GeorefMiscMethods`	wradlib xarray SubAccessor methods for Georef Misc Methods.

wradlib.georef.misc.bin_altitude(r, theta, sitealt, *, re=6371000, ke=1.3333333333333333)[source]#

wradlib.georef.misc.bin_altitude(obj: DataArray, **kwargs)

wradlib.georef.misc.bin_altitude(obj: Dataset, **kwargs)

Calculates the height of a radar bin taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the bin altitude is calculated as

\[h = \sqrt{r^2 + (k_e r_e)^2 + 2 r k_e r_e \sin\theta} - k_e r_e\]

Parameters:

r (numpy.ndarray) – Array of ranges [m]
theta (scalar or numpy.ndarray) – Array broadcastable to the shape of r elevation angles in degrees with 0° at horizontal and +90° pointing vertically upwards from the radar
sitealt (float) – Altitude in [m] a.s.l. of the referencing radar site
re (float, optional) – earth’s radius [m], defaults to 6371000.
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

altitude (numpy.ndarray) – Array of heights of the radar bins in [m]

wradlib.georef.misc.bin_distance(r, theta, sitealt, *, re=6371000, ke=1.3333333333333333)[source]#

wradlib.georef.misc.bin_distance(obj: DataArray, **kwargs)

wradlib.georef.misc.bin_distance(obj: Dataset, **kwargs)

Calculates great circle distance from radar site to radar bin over spherical earth, taking the refractivity of the atmosphere into account.

\[s = k_e r_e \arctan\left( \frac{r \cos\theta}{r \cos\theta + k_e r_e + h}\right)\]

where \(h\) would be the radar site altitude amsl.

Parameters:

r (numpy.ndarray) – Array of ranges [m]
theta (scalar or numpy.ndarray) – Array broadcastable to the shape of r elevation angles in degrees with 0° at horizontal and +90° pointing vertically upwards from the radar
sitealt (float) – site altitude [m] amsl.
re (float) – earth’s radius [m]
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

distance (numpy.ndarray) – Array of great circle arc distances [m]

wradlib.georef.misc.site_distance(r, theta, binalt, *, re=6371000, ke=1.3333333333333333)[source]#

wradlib.georef.misc.site_distance(obj: DataArray, **kwargs)

wradlib.georef.misc.site_distance(obj: Dataset, **kwargs)

Calculates great circle distance from bin at certain altitude to the radar site over spherical earth, taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the site distance may be calculated as

\[s = k_e r_e \arcsin\left( \frac{r \cos\theta}{k_e r_e + h_n(r, \theta, r_e, k_e)}\right)\]

where \(h_n\) would be provided by bin_altitude.

Parameters:

r (numpy.ndarray) – Array of ranges [m]
theta (scalar or numpy.ndarray) – Array broadcastable to the shape of r elevation angles in degrees with 0° at horizontal and +90° pointing vertically upwards from the radar
binalt (numpy.ndarray) – site altitude [m] amsl. same shape as r.
re (float, optional) – earth’s radius [m], defaults to 6371000.
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

distance (numpy.ndarray) – Array of great circle arc distances [m]

class wradlib.georef.misc.GeorefMiscMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Misc Methods.

bin_altitude(**kwargs)[source]#

Calculates the height of a radar bin taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the bin altitude is calculated as

\[h = \sqrt{r^2 + (k_e r_e)^2 + 2 r k_e r_e \sin\theta} - k_e r_e\]

Parameters:: obj (xarray.DataArray | xarray.Dataset) – DataArray
Returns:: z (xarray.DataArray) – DataArray

bin_distance(**kwargs)[source]#

Calculates great circle distance from radar site to radar bin over spherical earth, taking the refractivity of the atmosphere into account.

\[s = k_e r_e \arctan\left( \frac{r \cos\theta}{r \cos\theta + k_e r_e + h}\right)\]

where \(h\) would be the radar site altitude amsl.

Parameters:: obj (xarray.DataArray | xarray.Dataset) – DataArray or Dataset
Returns:: bin_distance (xarray.DataArray) – DataArray

site_distance(**kwargs)[source]#

Calculates great circle distance from bin at certain altitude to the radar site over spherical earth, taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the site distance may be calculated as

\[s = k_e r_e \arcsin\left( \frac{r \cos\theta}{k_e r_e + h_n(r, \theta, r_e, k_e)}\right)\]

where \(h_n\) would be provided by bin_altitude.

Parameters:: obj (xarray.DataArray | xarray.Dataset) – DataArray or Dataset
Returns:: z (xarray.DataArray) – DataArray

Polar Grid Functions#

`georeference`	Georeference Dataset/DataArray.
`spherical_to_xyz`	Transforms spherical coordinates (r, phi, theta) to cartesian coordinates (x, y, z) centered at site (aeqd).
`spherical_to_proj`	Transforms spherical coordinates (r, phi, theta) to projected coordinates centered at site in given projection.
`spherical_to_polyvert`	Generate 3-D polygon vertices directly from spherical coordinates (r, phi, theta).
`spherical_to_centroids`	Generate 3-D centroids of the radar bins from the sperical coordinates (r, phi, theta).
`centroid_to_polyvert`	Calculates the 2-D Polygon vertices necessary to form a rectangular polygon around the centroid's coordinates.
`sweep_centroids`	Construct sweep centroids native coordinates.
`maximum_intensity_projection`	Computes the maximum intensity projection along an arbitrary cut through the ppi from polar data.
`GeorefPolarMethods`	wradlib xarray SubAccessor methods for Georef Polar Methods.

wradlib.georef.polar.georeference(obj, **kwargs)[source]#

Georeference Dataset/DataArray.

Added in version 1.5.

This function adds georeference data to xarray Dataset/DataArray obj.

Parameters:

obj (xarray.Dataset or xarray.DataArray)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
re (float) – earth’s radius [m]
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

obj (xarray.Dataset or xarray.DataArray)

wradlib.georef.polar.spherical_to_xyz(r, phi, theta, site, *, re=None, ke=1.3333333333333333, **kwargs)[source]#

wradlib.georef.polar.spherical_to_xyz(obj: DataArray, **kwargs)

wradlib.georef.polar.spherical_to_xyz(obj: Dataset, **kwargs)

Transforms spherical coordinates (r, phi, theta) to cartesian coordinates (x, y, z) centered at site (aeqd).

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

r (numpy.ndarray) – Contains the radial distances in meters.
phi (numpy.ndarray) – Contains the azimuthal angles in degree.
theta (numpy.ndarray) – Contains the elevation angles in degree.
site (sequence) – the lon / lat / alt coordinates of the radar location and its altitude a.m.s.l. (in meters)
re (float, optional) – earth’s radius [m], defaults to None (calculating from given latitude).
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Keyword Arguments:

squeeze (bool, optional) – If True, returns squeezed array. Defaults to False.
strict_dims (bool, optional) – If True, generates output of (theta, phi, r, 3) in any case. If False, dimensions with same length are “merged”. Defaults to False.

Returns:

xyz (numpy.ndarray) – Array of shape (…, 3). Contains cartesian coordinates.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD.

wradlib.georef.polar.spherical_to_proj(r, phi, theta, site, *, crs=None, re=None, ke=1.3333333333333333)[source]#

wradlib.georef.polar.spherical_to_proj(obj: DataArray, **kwargs)

wradlib.georef.polar.spherical_to_proj(obj: Dataset, **kwargs)

Transforms spherical coordinates (r, phi, theta) to projected coordinates centered at site in given projection.

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

r (numpy.ndarray) – Contains the radial distances.
phi (numpy.ndarray) – Contains the azimuthal angles.
theta (numpy.ndarray) – Contains the elevation angles.
site (sequence) – the lon / lat coordinates of the radar location and its altitude a.m.s.l. (in meters) if site is of length two, altitude is assumed to be zero
crs – Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
re (float, optional) – earth’s radius [m], defaults to None (calculating from given latitude).
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

coords (numpy.ndarray) – Array of shape (…, 3). Contains projected map coordinates.

Examples

A few standard directions (North, South, North, East, South, West) with different distances (amounting to roughly 1°) from a site located at 48°N 9°E

>>> r  = np.array([0.,   0., 111., 111., 111., 111.,])*1000
>>> az = np.array([0., 180.,   0.,  90., 180., 270.,])
>>> th = np.array([0.,   0.,   0.,   0.,   0.,  0.5,])
>>> csite = (9.0, 48.0, 0)
>>> coords = spherical_to_proj(r, az, th, csite)
>>> for coord in coords:
...     print( '{0:7.4f}, {1:7.4f}, {2:7.4f}'.format(*coord))
...
 9.0000, 48.0000,  0.0000
 9.0000, 48.0000,  0.0000
 9.0000, 48.9981, 725.7160
10.4872, 47.9904, 725.7160
 9.0000, 47.0017, 725.7160
 7.5131, 47.9904, 1694.2234

Here, the coordinates of the east and west directions won’t come to lie on the latitude of the site because the beam doesn’t travel along the latitude circle but along a great circle.

See wradlib_workflow:georeferencing-and-projection.

wradlib.georef.polar.spherical_to_polyvert(r, phi, theta, site, *, crs=None)[source]#

wradlib.georef.polar.spherical_to_polyvert(obj: DataArray, **kwargs)

wradlib.georef.polar.spherical_to_polyvert(obj: Dataset, **kwargs)

Generate 3-D polygon vertices directly from spherical coordinates (r, phi, theta).

This is an alternative to centroid_to_polyvert which does not use centroids, but generates the polygon vertices by simply connecting the corners of the radar bins.

Both azimuth and range arrays are assumed to be equidistant and to contain only unique values. For further information refer to the documentation of spherical_to_xyz.

Parameters:

r (numpy.ndarray) – Array of ranges [m]; r defines the exterior boundaries of the range bins! (not the centroids). Thus, values must be positive!
phi (numpy.ndarray) – Array of azimuth angles containing values between 0° and 360°. The angles are assumed to describe the pointing direction fo the main beam lobe! The first angle can start at any values, but make sure the array is sorted continuously positively clockwise and the angles are equidistant. An angle if 0 degree is pointing north.
theta (float) – Elevation angle of scan
site (sequence) – the lon/lat/alt coordinates of the radar location
crs – Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj:pyproj.CRS instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.

Returns:

output (numpy.ndarray) – A 3-d array of polygon vertices with shape(num_vertices, num_vertex_nodes, 2). The last dimension carries the xyz-coordinates either in aeqd or given crs.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD, only returned if crs is None.

Examples

>>> import wradlib.georef as georef  # noqa
>>> import numpy as np
>>> from matplotlib import collections
>>> import matplotlib.pyplot as plt
>>> # define the polar coordinates and the site coordinates in lat/lon
>>> r = np.array([50., 100., 150., 200.]) * 1000
>>> az = np.array([0., 45., 90., 135., 180., 225., 270., 315.])
>>> el = 1.0
>>> site = (9.0, 48.0, 0)
>>> polygons, aeqd = georef.spherical_to_polyvert(r, az, el, site)
>>> # plot the resulting mesh
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> polycoll = collections.PolyCollection(polygons[...,:2], closed=True, facecolors='None')  # noqa
>>> ret = ax.add_collection(polycoll, autolim=True)
>>> plt.autoscale()
>>> plt.show()

wradlib.georef.polar.spherical_to_centroids(r, phi, theta, site, *, crs=None)[source]#

wradlib.georef.polar.spherical_to_centroids(obj: Dataset, **kwargs)

Generate 3-D centroids of the radar bins from the sperical coordinates (r, phi, theta).

Both azimuth and range arrays are assumed to be equidistant and to contain only unique values. The ranges are assumed to define the exterior boundaries of the range bins (thus they must be positive). The angles are assumed to describe the pointing direction fo the main beam lobe.

For further information refer to the documentation of spherical_to_xyz.

Parameters:

r (numpy.ndarray) – Array of ranges [m]; r defines the exterior boundaries of the range bins! (not the centroids). Thus, values must be positive!
phi (numpy.ndarray) – Array of azimuth angles containing values between 0° and 360°. The angles are assumed to describe the pointing direction fo the main beam lobe! The first angle can start at any values, but make sure the array is sorted continuously positively clockwise and the angles are equidistant. An angle if 0 degree is pointing north.
theta (float) – Elevation angle of scan
site (sequence) – the lon/lat/alt coordinates of the radar location
crs – Coordinate Reference System (CRS). Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.

Returns:

output (numpy.ndarray) – A 3-d array of bin centroids with shape(num_rays, num_bins, 3). The last dimension carries the xyz-coordinates either in aeqd or given crs.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD, only returned if crs is None.

Note

Azimuth angles of 360 deg are internally converted to 0 deg.

wradlib.georef.polar.centroid_to_polyvert(centroid, delta, /)[source]#

Calculates the 2-D Polygon vertices necessary to form a rectangular polygon around the centroid’s coordinates.

The vertices order will be clockwise, as this is the convention used by ESRI’s shapefile format for a polygon.

Parameters:

centroid (array-like) – List of 2-D coordinates of the center point of the rectangle.
delta (scalar or numpy.ndarray) – Symmetric distances of the vertices from the centroid in each direction. If delta is scalar, it is assumed to apply to both dimensions.

Returns:

vertices (numpy.ndarray) – An array with 5 vertices per centroid.

Note

The function can currently only deal with 2-D data (If you come up with a higher dimensional version of ‘clockwise’ you’re welcome to add it). The data is then assumed to be organized within the centroid array with the last dimension being the 2-D coordinates of each point.

Examples

>>> centroid_to_polyvert([0., 1.], [0.5, 1.5])
array([[-0.5, -0.5],
       [-0.5,  2.5],
       [ 0.5,  2.5],
       [ 0.5, -0.5],
       [-0.5, -0.5]])
>>> centroid_to_polyvert(np.arange(4).reshape((2,2)), 0.5)
array([[[-0.5,  0.5],
        [-0.5,  1.5],
        [ 0.5,  1.5],
        [ 0.5,  0.5],
        [-0.5,  0.5]],

       [[ 1.5,  2.5],
        [ 1.5,  3.5],
        [ 2.5,  3.5],
        [ 2.5,  2.5],
        [ 1.5,  2.5]]])

wradlib.georef.polar.sweep_centroids(nrays, rscale, nbins, elangle, /)[source]#

Construct sweep centroids native coordinates.

Parameters:

nrays (int) – number of rays
rscale (float) – length [m] of a range bin
nbins (int) – number of range bins
elangle (float) – elevation angle [deg]

Returns:

coordinates (numpy.ndarray) – array of shape (nrays,nbins,3) containing native centroid radar coordinates (slant range, azimuth, elevation)

wradlib.georef.polar.maximum_intensity_projection(data, *, r=None, az=None, angle=None, elev=None, autoext=True)[source]#

Computes the maximum intensity projection along an arbitrary cut through the ppi from polar data.

Parameters:

data (numpy.ndarray) – Array containing polar data (azimuth, range)
r (numpy.ndarray, optional) – Array containing range data
az (numpy.ndarray, optional) – Array containing azimuth data
angle (float, optional) – angle of slice, Defaults to 0. Should be between 0 and 180. 0. means horizontal slice, 90. means vertical slice
elev (float, optional) – elevation angle of scan, Defaults to 0.
autoext (bool, optional) – This routine uses numpy.numpy.digitize to bin the data. As this function needs bounds, we create one set of coordinates more than would usually be provided by r and az. Defaults to True.

Returns:

xs (numpy.ndarray) – meshgrid x array
ys (numpy.ndarray) – meshgrid y array
mip (numpy.ndarray) – Array containing the maximum intensity projection (range, range*2)

class wradlib.georef.polar.GeorefPolarMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Polar Methods.

georeference(**kwargs)[source]#

Georeference Dataset/DataArray.

Added in version 1.5.

This function adds georeference data to xarray Dataset/DataArray obj.

Parameters:

obj (xarray.Dataset or xarray.DataArray)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
re (float) – earth’s radius [m]
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

obj (xarray.Dataset or xarray.DataArray)

spherical_to_xyz(**kwargs)[source]#

Transforms spherical coordinates (r, phi, theta) to cartesian coordinates (x, y, z) centered at site (aeqd).

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

re (float) – earth’s radius [m], defaults to None (calculating from given latitude)
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

xyz (xarray.DataArray) – Array of shape (…, 3). Contains cartesian coordinates.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD.

spherical_to_proj(**kwargs)[source]#

Transforms spherical coordinates (r, phi, theta) to projected coordinates centered at site in given projection.

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

coords (xarray.DataArray) – Array of shape (…, 3). Contains projected map coordinates.

spherical_to_polyvert(**kwargs)[source]#

Generate 3-D polygon vertices directly from spherical coordinates (r, phi, theta).

This is an alternative to centroid_to_polyvert which does not use centroids, but generates the polygon vertices by simply connecting the corners of the radar bins.

Both azimuth and range arrays are assumed to be equidistant and to contain only unique values. For further information refer to the documentation of spherical_to_xyz.

Currently only works for PPI.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

xyz (xarray.DataArray) – Array of shape (…, 3). Contains cartesian coordinates.
crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD, only returned if crs is None.

spherical_to_centroids(**kwargs)[source]#

Generate 3-D centroids of the radar bins from the sperical coordinates (r, phi, theta).

For further information refer to the documentation of spherical_to_xyz.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.
PPI. (Currently only works for)

Returns:

xyz (xarray.DataArray) – Array of shape (…, 3). Contains cartesian coordinates.
crs (pyproj.crs.CoordinateSystem) – Reference System (CRS), AEQD, only returned if crs is None.

Note

Azimuth angles of 360 deg are internally converted to 0 deg.

Projection Functions#

`ensure_crs`	Return CRS object from given entry.
`reproject`	Transform coordinates from a source projection to a target projection.
`create_crs`	Conveniently supports the construction of pyproj Coordinate Reference System (CRS)
`create_osr`	Conveniently supports the construction of osr spatial reference objects
`projstr_to_osr`	Transform a PROJ string to an osr spatial reference object
`epsg_to_osr`	Create osr spatial reference object from EPSG number
`wkt_to_osr`	Create osr spatial reference object from WKT string
`get_default_projection`	Create a default projection object (wgs84)
`get_earth_radius`	Get the radius of the Earth (in km) for a given Spheroid model (sr) at a given position.
`get_radar_projection`	Get the native radar projection which is an azimuthal equidistant projection centered at the site using WGS84.
`get_earth_projection`	Get a default earth projection based on WGS
`get_extent`	Get the extent of 2d coordinates
`project_bounds`	Get geographic bounds in projected coordinate system
`meters_to_degrees`	Converts a distance in meters to degrees of latitude and longitude using the WGS84 ellipsoid.
`GeorefProjectionMethods`	wradlib xarray SubAccessor methods for Georef Projection Methods.

wradlib.georef.projection.ensure_crs(crs, trg='pyproj')[source]#

Return CRS object from given entry. Default to pyproj CRS.

Parameters:

crs – Coordinate Reference System (CRS) of the coordinates. Must be given and can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
None

Keyword Arguments:

trg (str, {"pyproj", "cartopy", "osr"}) – Target Coordinate Reference System type

Returns:

crs (pyproj.crs.CoordinateSystem, cartopy.crs.CRS, osgeo.osr.SpatialReference or None)

wradlib.georef.projection.reproject(*args, **kwargs)[source]#

wradlib.georef.projection.reproject(obj: Dataset, **kwargs)

wradlib.georef.projection.reproject(obj: DataArray, **kwargs)

Transform coordinates from a source projection to a target projection.

Call signatures:

reproject(C, **kwargs)
reproject(X, Y, **kwargs)
reproject(X, Y, Z, **kwargs)

C is the np array of source coordinates. X, Y and Z specify arrays of x, y, and z coordinate values

Parameters:

C (numpy.ndarray) – Array of shape (…,2) or (…,3) with coordinates (x,y) or (x,y,z) respectively
X (numpy.ndarray) – Array of x coordinates
Y (numpy.ndarray) – Array of y coordinates
Z (numpy.ndarray) – Array of z coordinates

Keyword Arguments:

src_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
trg_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
area_of_interest (tuple) – floats (WestLongitudeDeg, SouthLatitudeDeg, EastLongitudeDeg, NorthLatitudeDeg), defaults to None (no area of interest)

Returns:

trans (numpy.ndarray) – Array of reprojected coordinates x,y (…,2) or x,y,z (…,3) depending on input array.
X, Y (numpy.ndarray) – Arrays of reprojected x,y coordinates, shape depending on input array
X, Y, Z (numpy.ndarray) – Arrays of reprojected x,y,z coordinates, shape depending on input array

Examples

See notebooks/georeferencing/georef.

wradlib.georef.projection.create_crs(projname, **kwargs)[source]#

Conveniently supports the construction of pyproj Coordinate Reference System (CRS)

Currently, the following projection names (argument projname) are supported:

“aeqd”: Azimuthal Equidistant

needs the following keyword arguments:

lat_0 (latitude at projection center),

lon_0 (longitude at projection center),

x_0 (false Easting, also known as x-offset),

y_0 (false Northing, also known as y-offset)

“dwd-radolan” : RADOLAN Composite Coordinate System

no additional arguments needed.

Polar stereographic projection used by the German Weather Service (DWD) for all Radar composite products. See the final report on the RADOLAN project [DrWeigl et al., 2004] for details.

Parameters:

projname (str) – “aeqd” or “dwd-radolan”
kwargs (dict) – depends on projname - see above!

Returns:

output (pyproj.crs.CoordinateSystem) – pyproj Coordinate Reference System (CRS)

Examples

See wradlib_workflow:georeferencing-and-projection.

wradlib.georef.projection.create_osr(projname, **kwargs)[source]#

Conveniently supports the construction of osr spatial reference objects

Currently, the following projection names (argument projname) are supported:

“aeqd”: Azimuthal Equidistant

needs the following keyword arguments:

lat_0 (latitude at projection center),

lon_0 (longitude at projection center),

x_0 (false Easting, also known as x-offset),

y_0 (false Northing, also known as y-offset)

“dwd-radolan” : RADOLAN Composite Coordinate System

no additional arguments needed.

Polar stereographic projection used by the German Weather Service (DWD) for all Radar composite products. See the final report on the RADOLAN project [DrWeigl et al., 2004] for details.

Parameters:

projname (str) – “aeqd” or “dwd-radolan”
kwargs (dict) – depends on projname - see above!

Returns:

output (osgeo.osr.SpatialReference) – GDAL/OSR object defining projection

Examples

See wradlib_workflow:georeferencing-and-projection.

Deprecated since version 2.4: This will be removed in 3.0. Use wradlib.georef.projection.create_crs instead.

wradlib.georef.projection.projstr_to_osr(projstr)[source]#

Transform a PROJ string to an osr spatial reference object

Parameters:: projstr (str) – PROJ string describing projection
Returns:: crs (osgeo.osr.SpatialReference) – GDAL OSR SRS object defining projection

Examples

See radolan_grid:proj.

Deprecated since version 2.4: This will be removed in 3.0. Use pyproj.CRS.from_proj4 instead.

wradlib.georef.projection.epsg_to_osr(epsg=None)[source]#

Create osr spatial reference object from EPSG number

Parameters:: epsg (int) – EPSG-Number defining the coordinate system
Returns:: crs (osgeo.osr.SpatialReference) – GDAL/OSR object defining projection

Deprecated since version 2.4: This will be removed in 3.0. Use pyproj.CRS.from_epsg instead.

wradlib.georef.projection.wkt_to_osr(wkt=None)[source]#

Create osr spatial reference object from WKT string

Parameters:: wkt (str) – WTK string defining the coordinate reference system
Returns:: crs (osgeo.osr.SpatialReference) – GDAL/OSR object defining projection

Deprecated since version 2.4: This will be removed in 3.0. Use pyproj.CRS.from_wkt instead.

wradlib.georef.projection.get_default_projection()[source]#: Create a default projection object (wgs84)

wradlib.georef.projection.get_earth_radius(latitude, *, crs=None)[source]#

wradlib.georef.projection.get_earth_radius(obj: DataArray, *, crs=None)

wradlib.georef.projection.get_earth_radius(obj: Dataset, *, crs=None)

Get the radius of the Earth (in km) for a given Spheroid model (sr) at a given position.

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

Parameters:

latitude (float) – geodetic latitude in degrees

Keyword Arguments:

crs –

Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)

Defaults to EPSG(4326).

Returns:

radius (float) – earth radius in meter

wradlib.georef.projection.get_radar_projection(site)[source]#

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

Parameters:: site (sequence) – the WGS84 lon / lat coordinates of the radar location
Returns:: crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS) - radar centric AEQD

wradlib.georef.projection.get_earth_projection(model='ellipsoid', arcsecond=False)[source]#

Get a default earth projection based on WGS

Parameters:

model (str) – earth model used, defaults to ellipsoid:
- ‘ellipsoid’ - WGS84 with ellipsoid heights -> EPSG 4979
- ‘geoid’ - WGS84 with egm96 geoid heights -> EPSG 4326 + 5773
- ‘sphere’ - GRS 1980 authalic sphere -> EPSG 4047
arcsecond (boolean) – true to use arcsecond as unit instead of degree, defaults to False

Returns:

crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS)

wradlib.georef.projection.get_extent(coords)[source]#

Get the extent of 2d coordinates

Parameters:: coords (numpy.ndarray) – coordinates array with shape (…,(x,y))
Returns:: extent (tuple) – (xmin, xmax, ymin, ymax)

wradlib.georef.projection.project_bounds(bounds, crs)[source]#

Get geographic bounds in projected coordinate system

Parameters:

bounds (tuple of float) – (lon_min, lon_max, lat_min, lat_max) geographic bounds
crs – Coordinate Reference System (CRS) to be used for projection. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)

Returns:

bounds (tuople of float) – (xmin, xmax, ymin, ymax) projected bounds

wradlib.georef.projection.meters_to_degrees(meters, longitude=0.0, latitude=0.0)[source]#

Converts a distance in meters to degrees of latitude and longitude using the WGS84 ellipsoid. If scalar, assumes equal east/north offset.

Parameters:

meters (float or tuple(float, float)) –
Distance in meters.
- If scalar: interpreted as [meters, meters] (diagonal NE).
- If 2D: interpreted as [east, north] in meters.
latitude (float) – Reference latitude in degrees.
longitude (float) – Reference longitude in degrees.

Returns:

tuple – (delta_latitude, delta_longitude) in degrees

class wradlib.georef.projection.GeorefProjectionMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Projection Methods.

get_earth_radius(*, crs=None)[source]#

Get the radius of the Earth (in km) for a given Spheroid model (sr) at a given position.

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

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –

Coordinate Reference System (CRS) of the coordinates. Can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)

Defaults to EPSG(4326).

Returns:

radius (float) – earth radius in meter

reproject(**kwargs)[source]#

Transform coordinates from current projection to a target projection.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

src_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to source data CRS.
trg_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to WGS84.
coords (dict, optional) – Mapping of coordinates. Defaults to None
area_of_interest (tuple) – floats (WestLongitudeDeg, SouthLatitudeDeg, EastLongitudeDeg, NorthLatitudeDeg), defaults to None (no area of interest).

Returns:

obj (xarray.DataArray | xarray.Dataset) – reprojected Dataset/DataArray

Examples

See notebooks/georeferencing/georef.

Raster Functions#

`read_gdal_values`	Read values from a gdal object.
`read_gdal_projection`	Get a projection (OSR object) from a GDAL dataset.
`read_gdal_coordinates`	Get the projected coordinates from a GDAL dataset.
`extract_raster_dataset`	Extract data, coordinates and projection information
`get_raster_extent`	Get the coordinates of the 4 corners of the raster dataset
`get_raster_elevation`	Return surface elevation corresponding to raster dataset
`reproject_raster_dataset`	Reproject/Resample given dataset according to keyword arguments
`merge_raster_datasets`	Merge rasters.
`create_raster_dataset`	Create In-Memory Raster Dataset
`set_raster_origin`	Converts Data and Coordinates Origin
`set_raster_indexing`	Sets Data and Coordinates Indexing Scheme
`set_coordinate_indexing`	Sets Coordinates Indexing Scheme
`raster_to_polyvert`	Get raster polygonal vertices from gdal dataset.
`snap_bounds`	Adjusts integer bounds so they align with the resolution grid.
`snap_resolution`	Snap an integer resolution so it evenly divides the extent, staying close to the target.
`create_raster_xarray`	Create empty raster as xarray dataset following CF conventions
`create_raster_geographic`	Create an empty geographic raster as xarray dataset following CF conventions.
`get_raster_coordinates`	Generate a 2D array of raster coordinates and return it as a DataArray.
`add_raster_grid_mapping`

wradlib.georef.raster.read_gdal_values(dataset, *, nodata=None)[source]#

Read values from a gdal object.

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing
nodata (float) – replace nodata values

Returns:

values (numpy.ndarray) – Array of shape (nrows, ncols) or (nbands, nrows, ncols) containing the data values.

Examples

See notebooks/classify/clutter_cloud.

wradlib.georef.raster.read_gdal_projection(dataset)[source]#

Get a projection (OSR object) from a GDAL dataset.

Parameters:: dataset (osgeo.gdal.Dataset) – raster image with georeferencing
Returns:: crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS)

Examples

See notebooks/classify/clutter_cloud.

wradlib.georef.raster.read_gdal_coordinates(dataset, *, mode='center')[source]#

Get the projected coordinates from a GDAL dataset.

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing
mode (str) – either ‘center’ or ‘edge’

Returns:

coordinates (numpy.ndarray) – Array of shape (nrows,ncols,2) containing xy coordinates. The array indexing follows image convention with origin at upper left pixel. The shape is (nrows+1,ncols+1,2) if mode == edge.

Examples

See notebooks/classify/clutter_cloud.

wradlib.georef.raster.extract_raster_dataset(dataset, *, mode='center', nodata=None)[source]#

Extract data, coordinates and projection information

Parameters:

dataset (osgeo.gdal.Dataset) – raster dataset
mode (str) – either ‘center’ or ‘edge’
nodata (float) – replace nodata values

Returns:

values (numpy.ndarray) – Array of shape (nrows, ncols) or (nbands, nrows, ncols) containing the data values.
coords (numpy.ndarray) – Array of shape (nrows,ncols,2) containing xy coordinates. The array indexing follows image convention with origin at the upper left pixel (northup). The shape is (nrows+1,ncols+1,2) if mode == edge.
crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS) of the used coordinates.

wradlib.georef.raster.get_raster_extent(dataset, *, geo=False, window=True)[source]#

Get the coordinates of the 4 corners of the raster dataset

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)
geo (bool) – True to get geographical coordinates
window (bool) – True to get the window containing the corners

Returns:

extent (numpy.ndarray) – corner coordinates [ul,ll,lr,ur] or window extent [xmin, xmax, ymin, ymax]

wradlib.georef.raster.get_raster_elevation(dataset, *, resample=None, **kwargs)[source]#

Return surface elevation corresponding to raster dataset: The resampling algorithm is chosen based on scale ratio

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)
resample (gdal:osgeo.gdalconst.ResampleAlg) – If None the best algorithm is chosen based on scales. GRA_NearestNeighbour = 0, GRA_Bilinear = 1, GRA_Cubic = 2, GRA_CubicSpline = 3, GRA_Lanczos = 4, GRA_Average = 5, GRA_Mode = 6, GRA_Max = 8, GRA_Min = 9, GRA_Med = 10, GRA_Q1 = 11, GRA_Q3 = 12
kwargs (dict) – keyword arguments passed to wradlib.io.dem.get_srtm

Returns:

elevation (numpy.ndarray) – Array of shape (rows, cols, 2) containing elevation

wradlib.georef.raster.reproject_raster_dataset(src_ds, **kwargs)[source]#

Reproject/Resample given dataset according to keyword arguments

Parameters:

src_ds (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)

Keyword Arguments:

spacing (float) – float or tuple of two floats pixel spacing of destination dataset, same unit as pixel coordinates
size (int) – tuple of two ints X/YRasterSize of destination dataset
resample (gdal:osgeo.gdalconst.ResampleAlg) – defaults to GRA_Bilinear GRA_NearestNeighbour = 0, GRA_Bilinear = 1, GRA_Cubic = 2, GRA_CubicSpline = 3, GRA_Lanczos = 4, GRA_Average = 5, GRA_Mode = 6, GRA_Max = 8, GRA_Min = 9, GRA_Med = 10, GRA_Q1 = 11, GRA_Q3 = 12
src_crs –
Coordinate Reference System (CRS) of source dataset. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None (get projection from source dataset)
trg_crs –
Coordinate Reference System (CRS) of source dataset. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None.
align (bool or tuple) – If False, there is no destination grid aligment. If True, aligns the destination grid to the next integer multiple of destination grid. If tuple (upper-left x,y-coordinate), the destination grid is aligned to this point.

Returns:

dst_ds (osgeo.gdal.Dataset) – reprojected/resampled raster dataset

wradlib.georef.raster.merge_raster_datasets(datasets, **kwargs)[source]#

Merge rasters.

Parameters:

datasets (list) – list of osgeo.gdal.Dataset raster images with georeferencing
kwargs (dict) – keyword arguments passed to gdal.Warp()

Returns:

dataset (osgeo.gdal.Dataset) – merged raster dataset

wradlib.georef.raster.create_raster_dataset(data, coords, *, crs=None, nodata=-9999)[source]#

Create In-Memory Raster Dataset

Parameters:

data (numpy.ndarray) – Array of shape (rows, cols) or (bands, rows, cols) containing the data values.
coords (numpy.ndarray) – Array of shape (nrows, ncols, 2) containing pixel center coordinates or Array of shape (nrows+1, ncols+1, 2) containing pixel edge coordinates
crs – Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None.
nodata (int) – Value of NODATA

Returns:

dataset (osgeo.gdal.Dataset) – In-Memory raster dataset

Note

The origin of the provided data and coordinates is UPPER LEFT.

wradlib.georef.raster.set_raster_origin(data, coords, direction)[source]#

Converts Data and Coordinates Origin

Parameters:

data (numpy.ndarray) – Array of shape (rows, cols) or (bands, rows, cols) containing the data values.
coords (numpy.ndarray) – Array of shape (rows, cols, 2) containing xy-coordinates.
direction (str) – ‘lower’ or ‘upper’, direction in which to convert data and coordinates.

Returns:

data (numpy.ndarray) – Array of shape (rows, cols) or (bands, rows, cols) containing the data values.
coords (numpy.ndarray) – Array of shape (rows, cols, 2) containing xy-coordinates.

wradlib.georef.raster.set_raster_indexing(data, coords, *, indexing='xy')[source]#

Sets Data and Coordinates Indexing Scheme

This converts data and coordinate layout from row-major to column major indexing.

Parameters:

data (numpy.ndarray) – Array of shape (…, M, N) containing the data values.
coords (numpy.ndarray) – Array of shape (…, M, N, 2) containing xy-coordinates.
indexing (str) – ‘xy’ or ‘ij’, indexing scheme in which to convert data and coordinates.

Returns:

data (numpy.ndarray) – Array of shape (…, N, M) containing the data values.
coords (numpy.ndarray) – Array of shape (…, N, M, 2) containing xy-coordinates.

wradlib.georef.raster.set_coordinate_indexing(coords, *, indexing='xy')[source]#

Sets Coordinates Indexing Scheme

This converts coordinate layout from row-major to column major indexing.

Parameters:

coords (numpy.ndarray) – Array of shape (…, M, N, 2) containing xy-coordinates.
indexing (str) – ‘xy’ or ‘ij’, indexing scheme in which to convert data and coordinates.

Returns:

coords (numpy.ndarray) – Array of shape (…, N, M, 2) containing xy-coordinates.

wradlib.georef.raster.raster_to_polyvert(dataset)[source]#

Get raster polygonal vertices from gdal dataset.

Parameters:: dataset (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)
Returns:: polyvert (numpy.ndarray) – A N-d array of polygon vertices with shape (…, 5, 2).

wradlib.georef.raster.snap_bounds(bounds, resolution)[source]#

Adjusts integer bounds so they align with the resolution grid.

This ensures that the width and height of the bounds are evenly divisible by the resolution, centering the snapped bounds around the original center.

Parameters:

bounds (tuple of int) – (minx, miny, maxx, maxy) in integer units.
resolution (int or tuple of int) – Desired resolution per axis.

Returns:

tuple of int – Snapped bounds aligned to the resolution grid.

wradlib.georef.raster.snap_resolution(extent, target)[source]#

Snap an integer resolution so it evenly divides the extent, staying close to the target.

This function finds the closest integer resolution that divides the given extent exactly, while minimizing the difference from the desired target resolution.

Parameters:

extent (int) – Total extent along one axis (e.g., width or height in arcseconds or pixels).
target (int) – Desired resolution (size of each cell), must be a positive integer.

Returns:

int – Snapped resolution that evenly divides the extent and is closest to the target.

wradlib.georef.raster.create_raster_xarray(crs, bounds, resolution, snap_bounds=False, snap_resolution=False)[source]#

Create empty raster as xarray dataset following CF conventions

Parameters:

crs – Coordinate Reference System (CRS) mapping geographic to x,y coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
bounds (tuple of int) – Bounding box as (min_x, min_y, max_x, max_y), as integer.
resolution (int or tuple of int) – Grid resolution in x and y directions. If a single int is provided, it applies to both axes.
snap_bounds (bool, optional) – If True, adjusts bounds to align with resolution grid.
snap_resolution (bool, optional) – If True, adjusts resolution to evenly divide the extent.

Returns:

xarray.Dataset – An xarray Dataset with ‘x’ and ‘y’ pixel center coordinates, ‘spatial_ref’ coordinate with crs_wkt and GeoTransform attributes

wradlib.georef.raster.create_raster_geographic(bounds, resolution, resolution_in_meters=False)[source]#

Create an empty geographic raster as xarray dataset following CF conventions. The arc-second is used as unit to avoid precision issues.

Parameters:

bounds (tuple of int) – Bounding box in degrees: (min_lon, min_lat, max_lon, max_lat).
resolution (int or tuple of int) – Resolution value. Interpreted as meters if resolution_in_meters=True, otherwise as arcseconds.
resolution_in_meters (bool, optional) – If True, converts resolution from meters to arcseconds and snaps it to evenly divide the bounds.

Returns:

xarray.Dataset – xarray dataset with WGS84 CRS and coordinates using arcsecond unit, following CF conventions.

wradlib.georef.raster.get_raster_coordinates(ds, mode='center')[source]#

Generate a 2D array of raster coordinates and return it as a DataArray.

Parameters:

ds (xarray.Dataset) – A dataset containing 1D coordinates ‘x’ and ‘y’.
mode (str, optional) – Determines how coordinates are computed: - “center”: uses the center points of each grid cell (default). - “edge” : computes coordinates at the edges between cells.

Returns:

xarray.DataArray – A 3D DataArray of shape (y, x, 2), where the last dimension contains the spatial coordinates [x, y] for each grid cell.

wradlib.georef.raster.add_raster_grid_mapping(ds)[source]#

Rectangular Grid Functions#

`get_radolan_coords`	Calculates x,y coordinates of radolan grid from lon, lat
`get_radolan_coordinates`	Calculates x/y coordinates of radolan projection of the German Weather Service
`get_radolan_grid`	Calculates x/y coordinates of radolan grid of the German Weather Service
`xyz_to_spherical`
`grid_to_polyvert`	Get polygonal vertices from rectangular grid coordinates.
`GeorefRectMethods`	wradlib xarray SubAccessor methods for Georef Rect Methods.

wradlib.georef.rect.get_radolan_coords(lon, lat, **kwargs)[source]#

Calculates x,y coordinates of radolan grid from lon, lat

Parameters:

lon (float or numpy.ndarray) – longitude
lat (float, numpy.ndarray) – latitude

Keyword Arguments:

crs –

Coordinate Reference System (CRS) of the DWD grid with spheroid model. Can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
string “trig” to use trigonometric formulas for calculation (only for earth model - sphere).

Defaults to None (earth model - sphere).

wradlib.georef.rect.get_radolan_coordinates(nrows=None, ncols=None, **kwargs)[source]#

Calculates x/y coordinates of radolan projection of the German Weather Service

Returns the 1D x,y coordinates of the radolan projection for the given grid layout.

Parameters:

nrows (int) – number of rows (460, 900 by default, 1100, 1500)
ncols (int) – number of columns (460, 900 by default, 1400)

Keyword Arguments:

wgs84 (bool) – if True, output coordinates are in wgs84 lonlat format (default: False)
mode (str) – ‘radolan’ - lower left pixel coordinates ‘center’ - pixel center coordinates ‘edge’ - pixel edge coordinates
crs –
Coordinate Reference System (CRS) of the DWD grid with spheroid model. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
- string “trig” to use trigonometric formulas for calculation (only for earth model - sphere).
Defaults to None (earth model - sphere).

Returns:

radolan_ccords (tuple) – x and y 1D coordinate numpy.ndarray shape is (nrows, ) and (ncols, ) if mode=’radolan’ shape is (nrows, ) and (ncols, ) if mode=’center’ shape is (nrows+1, ) and (ncols+1, ) if mode=’edge’

wradlib.georef.rect.get_radolan_grid(nrows=None, ncols=None, **kwargs)[source]#

Calculates x/y coordinates of radolan grid of the German Weather Service

Returns the x,y coordinates of the radolan grid positions (lower left corner of every pixel). The radolan grid is a polarstereographic projection, the projection information was taken from RADOLAN-RADVOR-OP Kompositformat_2.2.2 [DWD, 2009]

Coordinates for 900km x 900km grid#
Coordinate	lon	lat	x	y
LowerLeft	3.5889E	46.9526N	-523.4622	-4658.645
LowerRight	14.6209E	47.0705N	376.5378	-4658.645
UpperRight	15.7208E	54.7405N	376.5378	-3758.645
UpperLeft	2.0715E	54.5877N	-523.4622	-3758.645

Coordinates for 1100km x 900km grid#
Coordinate	lon	lat	x	y
LowerLeft	4.6759E	46.1929N	-443.4622	-4758.645
LowerRight	15.4801E	46.1827N	456.5378	-4758.645
UpperRight	17.1128E	55.5342N	456.5378	-3658.645
UpperLeft	3.0889E	55.5482N	-443.4622	-3658.645

Coordinates for 1500km x 1400km grid#
Coordinate	lon	lat	x	y
LowerLeft	2.3419E	43.9336N	-673.4622	-5008.645

Parameters:

nrows (int) – number of rows (460, 900 by default, 1100, 1500)
ncols (int) – number of columns (460, 900 by default, 1400)

Keyword Arguments:

wgs84 (bool) – if True, output coordinates are in wgs84 lonlat format (default: False)
mode (str) – ‘radolan’ - lower left pixel coordinates ‘center’ - pixel center coordinates ‘edge’ - pixel edge coordinates
crs –
Coordinate Reference System (CRS) of the DWD grid with spheroid model. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
- string “trig” to use trigonometric formulas for calculation (only for earth model - sphere).
Defaults to None (earth model - sphere).

Returns:

radolan_grid (numpy.ndarray) – Array of xy- or lonlat-grid. shape is (nrows, ncols, 2) if mode=’radolan’ shape is (nrows, ncols, 2) if mode=’center’ shape is (nrows+1, ncols+1, 2) if mode=’edge’

Examples

>>> # using osr spatial reference transformation
>>> import wradlib.georef as georef  # noqa
>>> radolan_grid = georef.get_radolan_grid()
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(900, 900, 2), (-523.4622, -4658.6447)

>>> # using pure trigonometric transformations
>>> import wradlib.georef as georef
>>> radolan_grid = georef.get_radolan_grid(crs="trig")
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(900, 900, 2), (-523.4622, -4658.6447)

>>> # using osr spatial reference transformation
>>> import wradlib.georef as georef
>>> radolan_grid = georef.get_radolan_grid(1500, 1400)
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(1500, 1400, 2), (-673.4622, -5008.6447)

>>> # using osr spatial reference transformation
>>> import wradlib.georef as georef
>>> radolan_grid = georef.get_radolan_grid(900, 900, wgs84=True)
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(900, 900, 2), (3.5889, 46.9526)

See radolan_grid:polar-stereographic-projection.

Raises:: TypeError –

wradlib.georef.rect.xyz_to_spherical(*args, **kwargs)[source]#
wradlib.georef.rect.xyz_to_spherical(xyz: ndarray, *, altitude=0, crs=None, ke=1.3333333333333333)
wradlib.georef.rect.xyz_to_spherical(obj: DataArray, **kwargs)
wradlib.georef.rect.xyz_to_spherical(obj: Dataset, **kwargs)

wradlib.georef.rect.grid_to_polyvert(grid, *, ravel=False)[source]#

Get polygonal vertices from rectangular grid coordinates.

Parameters:

grid (numpy.ndarray) – grid edge coordinates
ravel (bool, optional) – option to flatten the grid, defaults to False.

Returns:

polyvert (numpy.ndarray) – A 3-d array of polygon vertices with shape (…, 5, 2).

class wradlib.georef.rect.GeorefRectMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Rect Methods.

xyz_to_spherical(**kwargs)[source]#

Returns spherical representation (r, theta, phi) of given cartesian coordinates (x, y, z) with respect to the reference altitude (asl) considering earth’s geometry (crs).

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the source coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None.
ke (float) – Adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

obj (xarray.Dataset) – obj with added spherical coordinates.

Satellite Functions#

`correct_parallax`
`dist_from_orbit`
`GeorefSatelliteMethods`	wradlib xarray SubAccessor methods for Georef Satellite Methods.

wradlib.georef.satellite.correct_parallax(*args, **kwargs)[source]#
wradlib.georef.satellite.correct_parallax(sr_xy: ndarray, nbin, drt, alpha)
wradlib.georef.satellite.correct_parallax(obj: Dataset, drt, **kwargs)

wradlib.georef.satellite.dist_from_orbit(*args, **kwargs)[source]#
wradlib.georef.satellite.dist_from_orbit(sr_alt: float, alpha, beta, r_sr_inv, *, re=6371000)
wradlib.georef.satellite.dist_from_orbit(obj: Dataset, bw_sr, freq, re)

class wradlib.georef.satellite.GeorefSatelliteMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Satellite Methods.

correct_parallax(drt, **kwargs)[source]#

Adjust the geolocations of the SR pixels

With SR, we refer to precipitation radars based on space-born platforms such as TRMM or GPM.

The sr_xy coordinates of the SR beam footprints need to be in the azimuthal equidistant projection of the ground radar. This ensures that the ground radar is fixed at xy-coordinate (0, 0), and every SR bin has its relative xy-coordinates with respect to the ground radar site.

Parameters:

obj (xarray.Dataset)
drt (float) – Gate lenght of SR in meter.

Returns:

obj (xarray.Dataset) – obj with added coordinates in ground radar projection and range

dist_from_orbit(bw_sr, freq, re)[source]#

Returns range distances of SR bins (in meters) as seen from the orbit

With SR, we refer to precipitation radars based on space-born platforms such as TRMM or GPM.

Parameters:

obj (xarray.Dataset)
bw_sr (float) – Beam width of SR in degree.
freq (int) – Frequency index of PR.
re (float) – earth radius [m]

Returns:

obj (xarray.Dataset) – obj with added PR bin range distances from SR platform in orbit.

Vector Functions (GDAL)#

`get_vector_coordinates`	Function iterates over gdal ogr layer features and packs extracted vector coordinate points into nested ndarray
`get_vector_points`	Extract coordinate points from given ogr geometry as generator object
`transform_geometry`	Perform geotransformation to given destination SpatialReferenceSystem
`ogr_create_layer`	Creates OGR.Layer objects in gdal.Dataset object.
`ogr_copy_layer`	Copy OGR.Layer object.
`ogr_copy_layer_by_name`	Copy OGR.Layer object.
`ogr_reproject_layer`	Reproject src_lyr to dst_lyr.
`ogr_add_feature`	Creates OGR.Feature objects in OGR.Layer object.
`ogr_add_geometry`	Copies single OGR.Geometry object to an OGR.Layer object.
`numpy_to_ogr`	Convert a vertex array to gdal/ogr geometry.
`ogr_to_numpy`	Backconvert a gdal/ogr geometry to a numpy vertex array.
`ogr_geocol_to_numpy`	Backconvert a gdal/ogr geometry Collection to a numpy vertex array.
`get_centroid`	Return centroid of a polygon

wradlib.georef.vector.get_vector_coordinates(layer, **kwargs)[source]#

Function iterates over gdal ogr layer features and packs extracted vector coordinate points into nested ndarray

It transforms coordinates to a given destination osr spatial reference if trg_crs is given and a geotransform is necessary.

Parameters:

layer (osgeo.ogr.Layer)

Keyword Arguments:

src_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Source Projection
trg_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Destination Projection
key (str) – attribute key to extract from layer feature

Returns:

shp (numpy.ndarray) – Dimension of subarrays Nx2 extracted shape coordinate points
attrs (list) – List of attributes extracted from features

wradlib.georef.vector.get_vector_points(geom)[source]#

Extract coordinate points from given ogr geometry as generator object

If geometries are nested, function recurses.

Parameters:: geom (osgeo.ogr.Geometry)
Yields:: result (numpy.ndarray) – expands to Nx2 dimensional nested point arrays

wradlib.georef.vector.transform_geometry(geom, trg_crs, **kwargs)[source]#

Perform geotransformation to given destination SpatialReferenceSystem

It transforms coordinates to a given destination osr spatial reference if a geotransform is neccessary.

Parameters:

geom (osgeo.ogr.Geometry)
trg_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Destination Projection

Keyword Arguments:

src_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Source Projection

Returns:

geom (osgeo.ogr.Geometry) – Transformed Geometry

wradlib.georef.vector.ogr_create_layer(ds, name, *, crs=None, geom_type=None, fields=None)[source]#

Creates OGR.Layer objects in gdal.Dataset object.

Creates one OGR.Layer with given name in given gdal.Dataset object using given OGR.GeometryType and FieldDefinitions

Parameters:

ds (osgeo.gdal.Dataset) – object
name (str) – OGRLayer name
crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input object
geom_type (gdal:osgeo.ogr.GeometryType) – (e.g. ogr.wkbPolygon)
fields (list) – list of 2 element tuples (str, gdal:osgeo.ogr.DataType) field name, field type

Returns:

out (osgeo.ogr.Layer) – object

wradlib.georef.vector.ogr_copy_layer(src_ds, index, dst_ds, *, reset=True)[source]#

Copy OGR.Layer object.

Copy OGR.Layer object from src_ds gdal.Dataset to dst_ds gdal.Dataset

Parameters:

src_ds (osgeo.gdal.Dataset) – object
index (int) – layer index
dst_ds (osgeo.gdal.Dataset) – object
reset (bool) – if True resets src_layer

wradlib.georef.vector.ogr_copy_layer_by_name(src_ds, name, dst_ds, *, reset=True)[source]#

Copy OGR.Layer object.

Copy OGR.Layer object from src_ds gdal.Dataset to dst_ds gdal.Dataset

Parameters:

src_ds (osgeo.gdal.Dataset) – object
name (str) – layer name
dst_ds (osgeo.gdal.Dataset) – object
reset (bool) – if True resets src_layer

wradlib.georef.vector.ogr_reproject_layer(src_lyr, dst_lyr, trg_crs, *, src_crs=None)[source]#

Reproject src_lyr to dst_lyr.

Creates one OGR.Layer with given name in given gdal.Dataset object using given OGR.GeometryType and FieldDefinitions

Parameters:

src_lyr (osgeo.ogr.Layer) – OGRLayer source layer
dst_lyr (osgeo.ogr.Layer) – OGRLayer destination layer
trg_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Projection Target crs
src_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Projection Source crs

Returns:

dst_lyr (osgeo.ogr.Layer) – OGRLayer destination layer

wradlib.georef.vector.ogr_add_feature(ds, src, *, name=None)[source]#

Creates OGR.Feature objects in OGR.Layer object.

OGR.Features are built from numpy src points or polygons.

OGR.Features ‘FID’ and ‘index’ corresponds to source data element

Parameters:

ds (osgeo.gdal.Dataset) – object
src (numpy.ndarray) – source data
name (str) – name of wanted Layer

wradlib.georef.vector.ogr_add_geometry(layer, geom, attrs)[source]#

Copies single OGR.Geometry object to an OGR.Layer object.

Given OGR.Geometry is copied to new OGR.Feature and written to given OGR.Layer by given index. Attributes are attached.

Parameters:

layer (osgeo.ogr.Layer) – object
geom (osgeo.ogr.Geometry) – object
attrs (list) – attributes referring to layer fields

wradlib.georef.vector.numpy_to_ogr(vert, geom_name)[source]#

Convert a vertex array to gdal/ogr geometry.

Using JSON as a vehicle to efficiently deal with numpy arrays.

Parameters:

vert (array-like) – a numpy array of vertices of shape (num vertices, 2)
geom_name (str) – Name of Geometry

Returns:

out (osgeo.ogr.Geometry) – object of type geom_name

wradlib.georef.vector.ogr_to_numpy(ogrobj)[source]#

Backconvert a gdal/ogr geometry to a numpy vertex array.

Using JSON as a vehicle to efficiently deal with numpy arrays.

Parameters:: ogrobj (osgeo.ogr.Geometry) – object
Returns:: out (numpy.ndarray) – a nested ndarray of vertices of shape (num vertices, 2)

wradlib.georef.vector.ogr_geocol_to_numpy(ogrobj)[source]#

Backconvert a gdal/ogr geometry Collection to a numpy vertex array.

This extracts only Polygon geometries!

Using JSON as a vehicle to efficiently deal with numpy arrays.

Parameters:: ogrobj (osgeo.ogr.Geometry) – Collection object
Returns:: out (numpy.ndarray) – a nested ndarray of vertices of shape (num vertices, 2)

wradlib.georef.vector.get_centroid(polyg)[source]#

Return centroid of a polygon

Parameters:: polyg (numpy.ndarray) – of shape (num vertices, 2) or ogr.Geometry object
Returns:: out (tuple) – x and y coordinate of the centroid

Xarray Functions#

`as_xarray_dataarray`	Create Xarray DataArray from NumPy Array
`create_xarray_dataarray`	Create Xarray DataArray from Polar Radar Data

wradlib.georef.xarray.as_xarray_dataarray(data, dims, coords)[source]#

Create Xarray DataArray from NumPy Array

Added in version 1.3.

Parameters:

data (numpy.ndarray) – data array
dims (dict) – dictionary describing xarray dimensions
coords (dict) – dictionary describing xarray coordinates

Returns:

dataset (xarray.DataArray) – DataArray

wradlib.georef.xarray.create_xarray_dataarray(data, *, r=None, phi=None, theta=None, site=None, sweep_mode='azimuth_surveillance', rf=1.0, **kwargs)[source]#

Create Xarray DataArray from Polar Radar Data

Added in version 1.3.

Parameters:

data (numpy.ndarray) – The data array. It is assumed that the first dimension is over the azimuth angles, while the second dimension is over the range bins
r (numpy.ndarray) – The ranges. Units may be chosen arbitrarily, m preferred.
phi (numpy.ndarray) – The azimuth angles in degrees.
theta (numpy.ndarray) – The elevation angles in degrees.
site (tuple) – Coordinates of the radar site.
sweep_mode (str) – Defaults to ‘azimuth_surveillance’.
rf (float) – factor to scale range, defaults to 1. (no scale)

Keyword Arguments:

re (float) – effective earth radius
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. Defaults to 4/3.
dim0 (str) – Name of the first dimension. Defaults to ‘azimuth’.
dim1 (str) – Name of the second dimension. Defaults to ‘range’.

Returns:

dataset (xarray.DataArray) – DataArray

class wradlib.georef.GeorefMethods(obj)[source]#

Bases: XarrayMethods, GeorefMiscMethods, GeorefPolarMethods, GeorefProjectionMethods, GeorefRectMethods, GeorefSatelliteMethods

wradlib xarray SubAccessor methods for Georef.

class wradlib.georef.GeorefMiscMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Misc Methods.

bin_altitude(**kwargs)[source]#

Calculates the height of a radar bin taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the bin altitude is calculated as

\[h = \sqrt{r^2 + (k_e r_e)^2 + 2 r k_e r_e \sin\theta} - k_e r_e\]

Parameters:: obj (xarray.DataArray | xarray.Dataset) – DataArray
Returns:: z (xarray.DataArray) – DataArray

bin_distance(**kwargs)[source]#

Calculates great circle distance from radar site to radar bin over spherical earth, taking the refractivity of the atmosphere into account.

\[s = k_e r_e \arctan\left( \frac{r \cos\theta}{r \cos\theta + k_e r_e + h}\right)\]

where \(h\) would be the radar site altitude amsl.

Parameters:: obj (xarray.DataArray | xarray.Dataset) – DataArray or Dataset
Returns:: bin_distance (xarray.DataArray) – DataArray

site_distance(**kwargs)[source]#

Calculates great circle distance from bin at certain altitude to the radar site over spherical earth, taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the site distance may be calculated as

\[s = k_e r_e \arcsin\left( \frac{r \cos\theta}{k_e r_e + h_n(r, \theta, r_e, k_e)}\right)\]

where \(h_n\) would be provided by bin_altitude.

Parameters:: obj (xarray.DataArray | xarray.Dataset) – DataArray or Dataset
Returns:: z (xarray.DataArray) – DataArray

class wradlib.georef.GeorefPolarMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Polar Methods.

georeference(**kwargs)[source]#

Georeference Dataset/DataArray.

Added in version 1.5.

This function adds georeference data to xarray Dataset/DataArray obj.

Parameters:

obj (xarray.Dataset or xarray.DataArray)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
re (float) – earth’s radius [m]
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

obj (xarray.Dataset or xarray.DataArray)

spherical_to_centroids(**kwargs)[source]#

Generate 3-D centroids of the radar bins from the sperical coordinates (r, phi, theta).

For further information refer to the documentation of spherical_to_xyz.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.
PPI. (Currently only works for)

Returns:

xyz (xarray.DataArray) – Array of shape (…, 3). Contains cartesian coordinates.
crs (pyproj.crs.CoordinateSystem) – Reference System (CRS), AEQD, only returned if crs is None.

Note

Azimuth angles of 360 deg are internally converted to 0 deg.

spherical_to_polyvert(**kwargs)[source]#

Generate 3-D polygon vertices directly from spherical coordinates (r, phi, theta).

This is an alternative to centroid_to_polyvert which does not use centroids, but generates the polygon vertices by simply connecting the corners of the radar bins.

Both azimuth and range arrays are assumed to be equidistant and to contain only unique values. For further information refer to the documentation of spherical_to_xyz.

Currently only works for PPI.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

xyz (xarray.DataArray) – Array of shape (…, 3). Contains cartesian coordinates.
crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD, only returned if crs is None.

spherical_to_proj(**kwargs)[source]#

Transforms spherical coordinates (r, phi, theta) to projected coordinates centered at site in given projection.

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

coords (xarray.DataArray) – Array of shape (…, 3). Contains projected map coordinates.

spherical_to_xyz(**kwargs)[source]#

Transforms spherical coordinates (r, phi, theta) to cartesian coordinates (x, y, z) centered at site (aeqd).

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

re (float) – earth’s radius [m], defaults to None (calculating from given latitude)
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

xyz (xarray.DataArray) – Array of shape (…, 3). Contains cartesian coordinates.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD.

class wradlib.georef.GeorefProjectionMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Projection Methods.

get_earth_radius(*, crs=None)[source]#

Get the radius of the Earth (in km) for a given Spheroid model (sr) at a given position.

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

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –

Coordinate Reference System (CRS) of the coordinates. Can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)

Defaults to EPSG(4326).

Returns:

radius (float) – earth radius in meter

reproject(**kwargs)[source]#

Transform coordinates from current projection to a target projection.

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

src_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to source data CRS.
trg_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to WGS84.
coords (dict, optional) – Mapping of coordinates. Defaults to None
area_of_interest (tuple) – floats (WestLongitudeDeg, SouthLatitudeDeg, EastLongitudeDeg, NorthLatitudeDeg), defaults to None (no area of interest).

Returns:

obj (xarray.DataArray | xarray.Dataset) – reprojected Dataset/DataArray

Examples

See notebooks/georeferencing/georef.

class wradlib.georef.GeorefRectMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Rect Methods.

xyz_to_spherical(**kwargs)[source]#

Returns spherical representation (r, theta, phi) of given cartesian coordinates (x, y, z) with respect to the reference altitude (asl) considering earth’s geometry (crs).

Parameters:

obj (xarray.DataArray | xarray.Dataset)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the source coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None.
ke (float) – Adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

obj (xarray.Dataset) – obj with added spherical coordinates.

class wradlib.georef.GeorefSatelliteMethods[source]#

Bases: object

wradlib xarray SubAccessor methods for Georef Satellite Methods.

correct_parallax(drt, **kwargs)[source]#

Adjust the geolocations of the SR pixels

With SR, we refer to precipitation radars based on space-born platforms such as TRMM or GPM.

Parameters:

obj (xarray.Dataset)
drt (float) – Gate lenght of SR in meter.

Returns:

obj (xarray.Dataset) – obj with added coordinates in ground radar projection and range

dist_from_orbit(bw_sr, freq, re)[source]#

Returns range distances of SR bins (in meters) as seen from the orbit

With SR, we refer to precipitation radars based on space-born platforms such as TRMM or GPM.

Parameters:

obj (xarray.Dataset)
bw_sr (float) – Beam width of SR in degree.
freq (int) – Frequency index of PR.
re (float) – earth radius [m]

Returns:

obj (xarray.Dataset) – obj with added PR bin range distances from SR platform in orbit.

wradlib.georef.add_raster_grid_mapping(ds)[source]#

wradlib.georef.as_xarray_dataarray(data, dims, coords)[source]#

Create Xarray DataArray from NumPy Array

Added in version 1.3.

Parameters:

data (numpy.ndarray) – data array
dims (dict) – dictionary describing xarray dimensions
coords (dict) – dictionary describing xarray coordinates

Returns:

dataset (xarray.DataArray) – DataArray

wradlib.georef.bin_altitude(r, theta, sitealt, *, re=6371000, ke=1.3333333333333333)[source]#

wradlib.georef.bin_altitude(obj: DataArray, **kwargs)

wradlib.georef.bin_altitude(obj: Dataset, **kwargs)

Calculates the height of a radar bin taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the bin altitude is calculated as

\[h = \sqrt{r^2 + (k_e r_e)^2 + 2 r k_e r_e \sin\theta} - k_e r_e\]

Parameters:

r (numpy.ndarray) – Array of ranges [m]
theta (scalar or numpy.ndarray) – Array broadcastable to the shape of r elevation angles in degrees with 0° at horizontal and +90° pointing vertically upwards from the radar
sitealt (float) – Altitude in [m] a.s.l. of the referencing radar site
re (float, optional) – earth’s radius [m], defaults to 6371000.
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

altitude (numpy.ndarray) – Array of heights of the radar bins in [m]

wradlib.georef.bin_distance(r, theta, sitealt, *, re=6371000, ke=1.3333333333333333)[source]#

wradlib.georef.bin_distance(obj: DataArray, **kwargs)

wradlib.georef.bin_distance(obj: Dataset, **kwargs)

Calculates great circle distance from radar site to radar bin over spherical earth, taking the refractivity of the atmosphere into account.

\[s = k_e r_e \arctan\left( \frac{r \cos\theta}{r \cos\theta + k_e r_e + h}\right)\]

where \(h\) would be the radar site altitude amsl.

Parameters:

r (numpy.ndarray) – Array of ranges [m]
theta (scalar or numpy.ndarray) – Array broadcastable to the shape of r elevation angles in degrees with 0° at horizontal and +90° pointing vertically upwards from the radar
sitealt (float) – site altitude [m] amsl.
re (float) – earth’s radius [m]
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

distance (numpy.ndarray) – Array of great circle arc distances [m]

wradlib.georef.centroid_to_polyvert(centroid, delta, /)[source]#

Calculates the 2-D Polygon vertices necessary to form a rectangular polygon around the centroid’s coordinates.

The vertices order will be clockwise, as this is the convention used by ESRI’s shapefile format for a polygon.

Parameters:

centroid (array-like) – List of 2-D coordinates of the center point of the rectangle.
delta (scalar or numpy.ndarray) – Symmetric distances of the vertices from the centroid in each direction. If delta is scalar, it is assumed to apply to both dimensions.

Returns:

vertices (numpy.ndarray) – An array with 5 vertices per centroid.

Note

Examples

>>> centroid_to_polyvert([0., 1.], [0.5, 1.5])
array([[-0.5, -0.5],
       [-0.5,  2.5],
       [ 0.5,  2.5],
       [ 0.5, -0.5],
       [-0.5, -0.5]])
>>> centroid_to_polyvert(np.arange(4).reshape((2,2)), 0.5)
array([[[-0.5,  0.5],
        [-0.5,  1.5],
        [ 0.5,  1.5],
        [ 0.5,  0.5],
        [-0.5,  0.5]],

       [[ 1.5,  2.5],
        [ 1.5,  3.5],
        [ 2.5,  3.5],
        [ 2.5,  2.5],
        [ 1.5,  2.5]]])

wradlib.georef.correct_parallax(*args, **kwargs)[source]#
wradlib.georef.correct_parallax(sr_xy: ndarray, nbin, drt, alpha)
wradlib.georef.correct_parallax(obj: Dataset, drt, **kwargs)

wradlib.georef.create_crs(projname, **kwargs)[source]#

Conveniently supports the construction of pyproj Coordinate Reference System (CRS)

Currently, the following projection names (argument projname) are supported:

“aeqd”: Azimuthal Equidistant

needs the following keyword arguments:

lat_0 (latitude at projection center),

lon_0 (longitude at projection center),

x_0 (false Easting, also known as x-offset),

y_0 (false Northing, also known as y-offset)

“dwd-radolan” : RADOLAN Composite Coordinate System

no additional arguments needed.

Polar stereographic projection used by the German Weather Service (DWD) for all Radar composite products. See the final report on the RADOLAN project [DrWeigl et al., 2004] for details.

Parameters:

projname (str) – “aeqd” or “dwd-radolan”
kwargs (dict) – depends on projname - see above!

Returns:

output (pyproj.crs.CoordinateSystem) – pyproj Coordinate Reference System (CRS)

Examples

See wradlib_workflow:georeferencing-and-projection.

wradlib.georef.create_osr(projname, **kwargs)[source]#

Conveniently supports the construction of osr spatial reference objects

Currently, the following projection names (argument projname) are supported:

“aeqd”: Azimuthal Equidistant

needs the following keyword arguments:

lat_0 (latitude at projection center),

lon_0 (longitude at projection center),

x_0 (false Easting, also known as x-offset),

y_0 (false Northing, also known as y-offset)

“dwd-radolan” : RADOLAN Composite Coordinate System

no additional arguments needed.

Polar stereographic projection used by the German Weather Service (DWD) for all Radar composite products. See the final report on the RADOLAN project [DrWeigl et al., 2004] for details.

Parameters:

projname (str) – “aeqd” or “dwd-radolan”
kwargs (dict) – depends on projname - see above!

Returns:

output (osgeo.osr.SpatialReference) – GDAL/OSR object defining projection

Examples

See wradlib_workflow:georeferencing-and-projection.

Deprecated since version 2.4: This will be removed in 3.0. Use wradlib.georef.projection.create_crs instead.

wradlib.georef.create_raster_dataset(data, coords, *, crs=None, nodata=-9999)[source]#

Create In-Memory Raster Dataset

Parameters:

data (numpy.ndarray) – Array of shape (rows, cols) or (bands, rows, cols) containing the data values.
coords (numpy.ndarray) – Array of shape (nrows, ncols, 2) containing pixel center coordinates or Array of shape (nrows+1, ncols+1, 2) containing pixel edge coordinates
crs – Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None.
nodata (int) – Value of NODATA

Returns:

dataset (osgeo.gdal.Dataset) – In-Memory raster dataset

Note

The origin of the provided data and coordinates is UPPER LEFT.

wradlib.georef.create_raster_geographic(bounds, resolution, resolution_in_meters=False)[source]#

Create an empty geographic raster as xarray dataset following CF conventions. The arc-second is used as unit to avoid precision issues.

Parameters:

bounds (tuple of int) – Bounding box in degrees: (min_lon, min_lat, max_lon, max_lat).
resolution (int or tuple of int) – Resolution value. Interpreted as meters if resolution_in_meters=True, otherwise as arcseconds.
resolution_in_meters (bool, optional) – If True, converts resolution from meters to arcseconds and snaps it to evenly divide the bounds.

Returns:

xarray.Dataset – xarray dataset with WGS84 CRS and coordinates using arcsecond unit, following CF conventions.

wradlib.georef.create_raster_xarray(crs, bounds, resolution, snap_bounds=False, snap_resolution=False)[source]#

Create empty raster as xarray dataset following CF conventions

Parameters:

crs – Coordinate Reference System (CRS) mapping geographic to x,y coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
bounds (tuple of int) – Bounding box as (min_x, min_y, max_x, max_y), as integer.
resolution (int or tuple of int) – Grid resolution in x and y directions. If a single int is provided, it applies to both axes.
snap_bounds (bool, optional) – If True, adjusts bounds to align with resolution grid.
snap_resolution (bool, optional) – If True, adjusts resolution to evenly divide the extent.

Returns:

xarray.Dataset – An xarray Dataset with ‘x’ and ‘y’ pixel center coordinates, ‘spatial_ref’ coordinate with crs_wkt and GeoTransform attributes

wradlib.georef.create_xarray_dataarray(data, *, r=None, phi=None, theta=None, site=None, sweep_mode='azimuth_surveillance', rf=1.0, **kwargs)[source]#

Create Xarray DataArray from Polar Radar Data

Added in version 1.3.

Parameters:

data (numpy.ndarray) – The data array. It is assumed that the first dimension is over the azimuth angles, while the second dimension is over the range bins
r (numpy.ndarray) – The ranges. Units may be chosen arbitrarily, m preferred.
phi (numpy.ndarray) – The azimuth angles in degrees.
theta (numpy.ndarray) – The elevation angles in degrees.
site (tuple) – Coordinates of the radar site.
sweep_mode (str) – Defaults to ‘azimuth_surveillance’.
rf (float) – factor to scale range, defaults to 1. (no scale)

Keyword Arguments:

re (float) – effective earth radius
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. Defaults to 4/3.
dim0 (str) – Name of the first dimension. Defaults to ‘azimuth’.
dim1 (str) – Name of the second dimension. Defaults to ‘range’.

Returns:

dataset (xarray.DataArray) – DataArray

wradlib.georef.dist_from_orbit(*args, **kwargs)[source]#
wradlib.georef.dist_from_orbit(sr_alt: float, alpha, beta, r_sr_inv, *, re=6371000)
wradlib.georef.dist_from_orbit(obj: Dataset, bw_sr, freq, re)

wradlib.georef.ensure_crs(crs, trg='pyproj')[source]#

Return CRS object from given entry. Default to pyproj CRS.

Parameters:

crs – Coordinate Reference System (CRS) of the coordinates. Must be given and can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
None

Keyword Arguments:

trg (str, {"pyproj", "cartopy", "osr"}) – Target Coordinate Reference System type

Returns:

crs (pyproj.crs.CoordinateSystem, cartopy.crs.CRS, osgeo.osr.SpatialReference or None)

wradlib.georef.epsg_to_osr(epsg=None)[source]#

Create osr spatial reference object from EPSG number

Parameters:: epsg (int) – EPSG-Number defining the coordinate system
Returns:: crs (osgeo.osr.SpatialReference) – GDAL/OSR object defining projection

Deprecated since version 2.4: This will be removed in 3.0. Use pyproj.CRS.from_epsg instead.

wradlib.georef.extract_raster_dataset(dataset, *, mode='center', nodata=None)[source]#

Extract data, coordinates and projection information

Parameters:

dataset (osgeo.gdal.Dataset) – raster dataset
mode (str) – either ‘center’ or ‘edge’
nodata (float) – replace nodata values

Returns:

values (numpy.ndarray) – Array of shape (nrows, ncols) or (nbands, nrows, ncols) containing the data values.
coords (numpy.ndarray) – Array of shape (nrows,ncols,2) containing xy coordinates. The array indexing follows image convention with origin at the upper left pixel (northup). The shape is (nrows+1,ncols+1,2) if mode == edge.
crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS) of the used coordinates.

wradlib.georef.georeference(obj, **kwargs)[source]#

Georeference Dataset/DataArray.

Added in version 1.5.

This function adds georeference data to xarray Dataset/DataArray obj.

Parameters:

obj (xarray.Dataset or xarray.DataArray)

Keyword Arguments:

crs –
Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.
re (float) – earth’s radius [m]
ke (float) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

obj (xarray.Dataset or xarray.DataArray)

wradlib.georef.get_centroid(polyg)[source]#

Return centroid of a polygon

Parameters:: polyg (numpy.ndarray) – of shape (num vertices, 2) or ogr.Geometry object
Returns:: out (tuple) – x and y coordinate of the centroid

wradlib.georef.get_default_projection()[source]#: Create a default projection object (wgs84)

wradlib.georef.get_earth_projection(model='ellipsoid', arcsecond=False)[source]#

Get a default earth projection based on WGS

Parameters:

model (str) – earth model used, defaults to ellipsoid:
- ‘ellipsoid’ - WGS84 with ellipsoid heights -> EPSG 4979
- ‘geoid’ - WGS84 with egm96 geoid heights -> EPSG 4326 + 5773
- ‘sphere’ - GRS 1980 authalic sphere -> EPSG 4047
arcsecond (boolean) – true to use arcsecond as unit instead of degree, defaults to False

Returns:

crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS)

wradlib.georef.get_earth_radius(latitude, *, crs=None)[source]#

wradlib.georef.get_earth_radius(obj: DataArray, *, crs=None)

wradlib.georef.get_earth_radius(obj: Dataset, *, crs=None)

Get the radius of the Earth (in km) for a given Spheroid model (sr) at a given position.

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

Parameters:

latitude (float) – geodetic latitude in degrees

Keyword Arguments:

crs –

Coordinate Reference System (CRS) of the coordinates. Must be provided and can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)

Defaults to EPSG(4326).

Returns:

radius (float) – earth radius in meter

wradlib.georef.get_extent(coords)[source]#

Get the extent of 2d coordinates

Parameters:: coords (numpy.ndarray) – coordinates array with shape (…,(x,y))
Returns:: extent (tuple) – (xmin, xmax, ymin, ymax)

wradlib.georef.get_radar_projection(site)[source]#

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

Parameters:: site (sequence) – the WGS84 lon / lat coordinates of the radar location
Returns:: crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS) - radar centric AEQD

wradlib.georef.get_radolan_coordinates(nrows=None, ncols=None, **kwargs)[source]#

Calculates x/y coordinates of radolan projection of the German Weather Service

Returns the 1D x,y coordinates of the radolan projection for the given grid layout.

Parameters:

nrows (int) – number of rows (460, 900 by default, 1100, 1500)
ncols (int) – number of columns (460, 900 by default, 1400)

Keyword Arguments:

wgs84 (bool) – if True, output coordinates are in wgs84 lonlat format (default: False)
mode (str) – ‘radolan’ - lower left pixel coordinates ‘center’ - pixel center coordinates ‘edge’ - pixel edge coordinates
crs –
Coordinate Reference System (CRS) of the DWD grid with spheroid model. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
- string “trig” to use trigonometric formulas for calculation (only for earth model - sphere).
Defaults to None (earth model - sphere).

Returns:

wradlib.georef.get_radolan_coords(lon, lat, **kwargs)[source]#

Calculates x,y coordinates of radolan grid from lon, lat

Parameters:

lon (float or numpy.ndarray) – longitude
lat (float, numpy.ndarray) – latitude

Keyword Arguments:

crs –

Coordinate Reference System (CRS) of the DWD grid with spheroid model. Can be one of:

A pyproj.crs.CoordinateSystem instance
A cartopy.crs.CRS instance
A osgeo.osr.SpatialReference instance
A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
string “trig” to use trigonometric formulas for calculation (only for earth model - sphere).

Defaults to None (earth model - sphere).

wradlib.georef.get_radolan_grid(nrows=None, ncols=None, **kwargs)[source]#

Calculates x/y coordinates of radolan grid of the German Weather Service

Coordinates for 900km x 900km grid#
Coordinate	lon	lat	x	y
LowerLeft	3.5889E	46.9526N	-523.4622	-4658.645
LowerRight	14.6209E	47.0705N	376.5378	-4658.645
UpperRight	15.7208E	54.7405N	376.5378	-3758.645
UpperLeft	2.0715E	54.5877N	-523.4622	-3758.645

Coordinates for 1100km x 900km grid#
Coordinate	lon	lat	x	y
LowerLeft	4.6759E	46.1929N	-443.4622	-4758.645
LowerRight	15.4801E	46.1827N	456.5378	-4758.645
UpperRight	17.1128E	55.5342N	456.5378	-3658.645
UpperLeft	3.0889E	55.5482N	-443.4622	-3658.645

Coordinates for 1500km x 1400km grid#
Coordinate	lon	lat	x	y
LowerLeft	2.3419E	43.9336N	-673.4622	-5008.645

Parameters:

nrows (int) – number of rows (460, 900 by default, 1100, 1500)
ncols (int) – number of columns (460, 900 by default, 1400)

Keyword Arguments:

wgs84 (bool) – if True, output coordinates are in wgs84 lonlat format (default: False)
mode (str) – ‘radolan’ - lower left pixel coordinates ‘center’ - pixel center coordinates ‘edge’ - pixel edge coordinates
crs –
Coordinate Reference System (CRS) of the DWD grid with spheroid model. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
- string “trig” to use trigonometric formulas for calculation (only for earth model - sphere).
Defaults to None (earth model - sphere).

Returns:

Examples

>>> # using osr spatial reference transformation
>>> import wradlib.georef as georef  # noqa
>>> radolan_grid = georef.get_radolan_grid()
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(900, 900, 2), (-523.4622, -4658.6447)

>>> # using pure trigonometric transformations
>>> import wradlib.georef as georef
>>> radolan_grid = georef.get_radolan_grid(crs="trig")
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(900, 900, 2), (-523.4622, -4658.6447)

>>> # using osr spatial reference transformation
>>> import wradlib.georef as georef
>>> radolan_grid = georef.get_radolan_grid(1500, 1400)
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(1500, 1400, 2), (-673.4622, -5008.6447)

>>> # using osr spatial reference transformation
>>> import wradlib.georef as georef
>>> radolan_grid = georef.get_radolan_grid(900, 900, wgs84=True)
>>> print("{0}, ({1:.4f}, {2:.4f})".format(radolan_grid.shape, *radolan_grid[0,0,:]))  # noqa
(900, 900, 2), (3.5889, 46.9526)

See radolan_grid:polar-stereographic-projection.

Raises:: TypeError –

wradlib.georef.get_raster_coordinates(ds, mode='center')[source]#

Generate a 2D array of raster coordinates and return it as a DataArray.

Parameters:

ds (xarray.Dataset) – A dataset containing 1D coordinates ‘x’ and ‘y’.
mode (str, optional) – Determines how coordinates are computed: - “center”: uses the center points of each grid cell (default). - “edge” : computes coordinates at the edges between cells.

Returns:

xarray.DataArray – A 3D DataArray of shape (y, x, 2), where the last dimension contains the spatial coordinates [x, y] for each grid cell.

wradlib.georef.get_raster_elevation(dataset, *, resample=None, **kwargs)[source]#

Return surface elevation corresponding to raster dataset: The resampling algorithm is chosen based on scale ratio

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)
resample (gdal:osgeo.gdalconst.ResampleAlg) – If None the best algorithm is chosen based on scales. GRA_NearestNeighbour = 0, GRA_Bilinear = 1, GRA_Cubic = 2, GRA_CubicSpline = 3, GRA_Lanczos = 4, GRA_Average = 5, GRA_Mode = 6, GRA_Max = 8, GRA_Min = 9, GRA_Med = 10, GRA_Q1 = 11, GRA_Q3 = 12
kwargs (dict) – keyword arguments passed to wradlib.io.dem.get_srtm

Returns:

elevation (numpy.ndarray) – Array of shape (rows, cols, 2) containing elevation

wradlib.georef.get_raster_extent(dataset, *, geo=False, window=True)[source]#

Get the coordinates of the 4 corners of the raster dataset

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)
geo (bool) – True to get geographical coordinates
window (bool) – True to get the window containing the corners

Returns:

extent (numpy.ndarray) – corner coordinates [ul,ll,lr,ur] or window extent [xmin, xmax, ymin, ymax]

wradlib.georef.get_vector_coordinates(layer, **kwargs)[source]#

Function iterates over gdal ogr layer features and packs extracted vector coordinate points into nested ndarray

It transforms coordinates to a given destination osr spatial reference if trg_crs is given and a geotransform is necessary.

Parameters:

layer (osgeo.ogr.Layer)

Keyword Arguments:

src_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Source Projection
trg_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Destination Projection
key (str) – attribute key to extract from layer feature

Returns:

shp (numpy.ndarray) – Dimension of subarrays Nx2 extracted shape coordinate points
attrs (list) – List of attributes extracted from features

wradlib.georef.get_vector_points(geom)[source]#

Extract coordinate points from given ogr geometry as generator object

If geometries are nested, function recurses.

Parameters:: geom (osgeo.ogr.Geometry)
Yields:: result (numpy.ndarray) – expands to Nx2 dimensional nested point arrays

wradlib.georef.grid_to_polyvert(grid, *, ravel=False)[source]#

Get polygonal vertices from rectangular grid coordinates.

Parameters:

grid (numpy.ndarray) – grid edge coordinates
ravel (bool, optional) – option to flatten the grid, defaults to False.

Returns:

polyvert (numpy.ndarray) – A 3-d array of polygon vertices with shape (…, 5, 2).

wradlib.georef.maximum_intensity_projection(data, *, r=None, az=None, angle=None, elev=None, autoext=True)[source]#

Computes the maximum intensity projection along an arbitrary cut through the ppi from polar data.

Parameters:

data (numpy.ndarray) – Array containing polar data (azimuth, range)
r (numpy.ndarray, optional) – Array containing range data
az (numpy.ndarray, optional) – Array containing azimuth data
angle (float, optional) – angle of slice, Defaults to 0. Should be between 0 and 180. 0. means horizontal slice, 90. means vertical slice
elev (float, optional) – elevation angle of scan, Defaults to 0.
autoext (bool, optional) – This routine uses numpy.numpy.digitize to bin the data. As this function needs bounds, we create one set of coordinates more than would usually be provided by r and az. Defaults to True.

Returns:

xs (numpy.ndarray) – meshgrid x array
ys (numpy.ndarray) – meshgrid y array
mip (numpy.ndarray) – Array containing the maximum intensity projection (range, range*2)

wradlib.georef.merge_raster_datasets(datasets, **kwargs)[source]#

Merge rasters.

Parameters:

datasets (list) – list of osgeo.gdal.Dataset raster images with georeferencing
kwargs (dict) – keyword arguments passed to gdal.Warp()

Returns:

dataset (osgeo.gdal.Dataset) – merged raster dataset

wradlib.georef.meters_to_degrees(meters, longitude=0.0, latitude=0.0)[source]#

Converts a distance in meters to degrees of latitude and longitude using the WGS84 ellipsoid. If scalar, assumes equal east/north offset.

Parameters:

meters (float or tuple(float, float)) –
Distance in meters.
- If scalar: interpreted as [meters, meters] (diagonal NE).
- If 2D: interpreted as [east, north] in meters.
latitude (float) – Reference latitude in degrees.
longitude (float) – Reference longitude in degrees.

Returns:

tuple – (delta_latitude, delta_longitude) in degrees

wradlib.georef.numpy_to_ogr(vert, geom_name)[source]#

Convert a vertex array to gdal/ogr geometry.

Using JSON as a vehicle to efficiently deal with numpy arrays.

Parameters:

vert (array-like) – a numpy array of vertices of shape (num vertices, 2)
geom_name (str) – Name of Geometry

Returns:

out (osgeo.ogr.Geometry) – object of type geom_name

wradlib.georef.ogr_add_feature(ds, src, *, name=None)[source]#

Creates OGR.Feature objects in OGR.Layer object.

OGR.Features are built from numpy src points or polygons.

OGR.Features ‘FID’ and ‘index’ corresponds to source data element

Parameters:

ds (osgeo.gdal.Dataset) – object
src (numpy.ndarray) – source data
name (str) – name of wanted Layer

wradlib.georef.ogr_add_geometry(layer, geom, attrs)[source]#

Copies single OGR.Geometry object to an OGR.Layer object.

Given OGR.Geometry is copied to new OGR.Feature and written to given OGR.Layer by given index. Attributes are attached.

Parameters:

layer (osgeo.ogr.Layer) – object
geom (osgeo.ogr.Geometry) – object
attrs (list) – attributes referring to layer fields

wradlib.georef.ogr_copy_layer(src_ds, index, dst_ds, *, reset=True)[source]#

Copy OGR.Layer object.

Copy OGR.Layer object from src_ds gdal.Dataset to dst_ds gdal.Dataset

Parameters:

src_ds (osgeo.gdal.Dataset) – object
index (int) – layer index
dst_ds (osgeo.gdal.Dataset) – object
reset (bool) – if True resets src_layer

wradlib.georef.ogr_copy_layer_by_name(src_ds, name, dst_ds, *, reset=True)[source]#

Copy OGR.Layer object.

Copy OGR.Layer object from src_ds gdal.Dataset to dst_ds gdal.Dataset

Parameters:

src_ds (osgeo.gdal.Dataset) – object
name (str) – layer name
dst_ds (osgeo.gdal.Dataset) – object
reset (bool) – if True resets src_layer

wradlib.georef.ogr_create_layer(ds, name, *, crs=None, geom_type=None, fields=None)[source]#

Creates OGR.Layer objects in gdal.Dataset object.

Creates one OGR.Layer with given name in given gdal.Dataset object using given OGR.GeometryType and FieldDefinitions

Parameters:

ds (osgeo.gdal.Dataset) – object
name (str) – OGRLayer name
crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input object
geom_type (gdal:osgeo.ogr.GeometryType) – (e.g. ogr.wkbPolygon)
fields (list) – list of 2 element tuples (str, gdal:osgeo.ogr.DataType) field name, field type

Returns:

out (osgeo.ogr.Layer) – object

wradlib.georef.ogr_geocol_to_numpy(ogrobj)[source]#

Backconvert a gdal/ogr geometry Collection to a numpy vertex array.

This extracts only Polygon geometries!

Using JSON as a vehicle to efficiently deal with numpy arrays.

Parameters:: ogrobj (osgeo.ogr.Geometry) – Collection object
Returns:: out (numpy.ndarray) – a nested ndarray of vertices of shape (num vertices, 2)

wradlib.georef.ogr_reproject_layer(src_lyr, dst_lyr, trg_crs, *, src_crs=None)[source]#

Reproject src_lyr to dst_lyr.

Creates one OGR.Layer with given name in given gdal.Dataset object using given OGR.GeometryType and FieldDefinitions

Parameters:

src_lyr (osgeo.ogr.Layer) – OGRLayer source layer
dst_lyr (osgeo.ogr.Layer) – OGRLayer destination layer
trg_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Projection Target crs
src_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Projection Source crs

Returns:

dst_lyr (osgeo.ogr.Layer) – OGRLayer destination layer

wradlib.georef.ogr_to_numpy(ogrobj)[source]#

Backconvert a gdal/ogr geometry to a numpy vertex array.

Using JSON as a vehicle to efficiently deal with numpy arrays.

Parameters:: ogrobj (osgeo.ogr.Geometry) – object
Returns:: out (numpy.ndarray) – a nested ndarray of vertices of shape (num vertices, 2)

wradlib.georef.project_bounds(bounds, crs)[source]#

Get geographic bounds in projected coordinate system

Parameters:

bounds (tuple of float) – (lon_min, lon_max, lat_min, lat_max) geographic bounds
crs – Coordinate Reference System (CRS) to be used for projection. Must be provided and can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)

Returns:

bounds (tuople of float) – (xmin, xmax, ymin, ymax) projected bounds

wradlib.georef.projstr_to_osr(projstr)[source]#

Transform a PROJ string to an osr spatial reference object

Parameters:: projstr (str) – PROJ string describing projection
Returns:: crs (osgeo.osr.SpatialReference) – GDAL OSR SRS object defining projection

Examples

See radolan_grid:proj.

Deprecated since version 2.4: This will be removed in 3.0. Use pyproj.CRS.from_proj4 instead.

wradlib.georef.raster_to_polyvert(dataset)[source]#

Get raster polygonal vertices from gdal dataset.

Parameters:: dataset (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)
Returns:: polyvert (numpy.ndarray) – A N-d array of polygon vertices with shape (…, 5, 2).

wradlib.georef.read_gdal_coordinates(dataset, *, mode='center')[source]#

Get the projected coordinates from a GDAL dataset.

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing
mode (str) – either ‘center’ or ‘edge’

Returns:

Examples

See notebooks/classify/clutter_cloud.

wradlib.georef.read_gdal_projection(dataset)[source]#

Get a projection (OSR object) from a GDAL dataset.

Parameters:: dataset (osgeo.gdal.Dataset) – raster image with georeferencing
Returns:: crs (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS)

Examples

See notebooks/classify/clutter_cloud.

wradlib.georef.read_gdal_values(dataset, *, nodata=None)[source]#

Read values from a gdal object.

Parameters:

dataset (osgeo.gdal.Dataset) – raster image with georeferencing
nodata (float) – replace nodata values

Returns:

values (numpy.ndarray) – Array of shape (nrows, ncols) or (nbands, nrows, ncols) containing the data values.

Examples

See notebooks/classify/clutter_cloud.

wradlib.georef.reproject(*args, **kwargs)[source]#

wradlib.georef.reproject(obj: Dataset, **kwargs)

wradlib.georef.reproject(obj: DataArray, **kwargs)

Transform coordinates from a source projection to a target projection.

Call signatures:

reproject(C, **kwargs)
reproject(X, Y, **kwargs)
reproject(X, Y, Z, **kwargs)

C is the np array of source coordinates. X, Y and Z specify arrays of x, y, and z coordinate values

Parameters:

C (numpy.ndarray) – Array of shape (…,2) or (…,3) with coordinates (x,y) or (x,y,z) respectively
X (numpy.ndarray) – Array of x coordinates
Y (numpy.ndarray) – Array of y coordinates
Z (numpy.ndarray) – Array of z coordinates

Keyword Arguments:

src_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
trg_crs –
Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
area_of_interest (tuple) – floats (WestLongitudeDeg, SouthLatitudeDeg, EastLongitudeDeg, NorthLatitudeDeg), defaults to None (no area of interest)

Returns:

trans (numpy.ndarray) – Array of reprojected coordinates x,y (…,2) or x,y,z (…,3) depending on input array.
X, Y (numpy.ndarray) – Arrays of reprojected x,y coordinates, shape depending on input array
X, Y, Z (numpy.ndarray) – Arrays of reprojected x,y,z coordinates, shape depending on input array

Examples

See notebooks/georeferencing/georef.

wradlib.georef.reproject_raster_dataset(src_ds, **kwargs)[source]#

Reproject/Resample given dataset according to keyword arguments

Parameters:

src_ds (osgeo.gdal.Dataset) – raster image with georeferencing (GeoTransform at least)

Keyword Arguments:

spacing (float) – float or tuple of two floats pixel spacing of destination dataset, same unit as pixel coordinates
size (int) – tuple of two ints X/YRasterSize of destination dataset
resample (gdal:osgeo.gdalconst.ResampleAlg) – defaults to GRA_Bilinear GRA_NearestNeighbour = 0, GRA_Bilinear = 1, GRA_Cubic = 2, GRA_CubicSpline = 3, GRA_Lanczos = 4, GRA_Average = 5, GRA_Mode = 6, GRA_Max = 8, GRA_Min = 9, GRA_Med = 10, GRA_Q1 = 11, GRA_Q3 = 12
src_crs –
Coordinate Reference System (CRS) of source dataset. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None (get projection from source dataset)
trg_crs –
Coordinate Reference System (CRS) of source dataset. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to None.
align (bool or tuple) – If False, there is no destination grid aligment. If True, aligns the destination grid to the next integer multiple of destination grid. If tuple (upper-left x,y-coordinate), the destination grid is aligned to this point.

Returns:

dst_ds (osgeo.gdal.Dataset) – reprojected/resampled raster dataset

wradlib.georef.set_coordinate_indexing(coords, *, indexing='xy')[source]#

Sets Coordinates Indexing Scheme

This converts coordinate layout from row-major to column major indexing.

Parameters:

coords (numpy.ndarray) – Array of shape (…, M, N, 2) containing xy-coordinates.
indexing (str) – ‘xy’ or ‘ij’, indexing scheme in which to convert data and coordinates.

Returns:

coords (numpy.ndarray) – Array of shape (…, N, M, 2) containing xy-coordinates.

wradlib.georef.set_raster_indexing(data, coords, *, indexing='xy')[source]#

Sets Data and Coordinates Indexing Scheme

This converts data and coordinate layout from row-major to column major indexing.

Parameters:

data (numpy.ndarray) – Array of shape (…, M, N) containing the data values.
coords (numpy.ndarray) – Array of shape (…, M, N, 2) containing xy-coordinates.
indexing (str) – ‘xy’ or ‘ij’, indexing scheme in which to convert data and coordinates.

Returns:

data (numpy.ndarray) – Array of shape (…, N, M) containing the data values.
coords (numpy.ndarray) – Array of shape (…, N, M, 2) containing xy-coordinates.

wradlib.georef.set_raster_origin(data, coords, direction)[source]#

Converts Data and Coordinates Origin

Parameters:

data (numpy.ndarray) – Array of shape (rows, cols) or (bands, rows, cols) containing the data values.
coords (numpy.ndarray) – Array of shape (rows, cols, 2) containing xy-coordinates.
direction (str) – ‘lower’ or ‘upper’, direction in which to convert data and coordinates.

Returns:

data (numpy.ndarray) – Array of shape (rows, cols) or (bands, rows, cols) containing the data values.
coords (numpy.ndarray) – Array of shape (rows, cols, 2) containing xy-coordinates.

wradlib.georef.site_distance(r, theta, binalt, *, re=6371000, ke=1.3333333333333333)[source]#

wradlib.georef.site_distance(obj: DataArray, **kwargs)

wradlib.georef.site_distance(obj: Dataset, **kwargs)

Calculates great circle distance from bin at certain altitude to the radar site over spherical earth, taking the refractivity of the atmosphere into account.

Based on [Doviak et al., 1993] the site distance may be calculated as

\[s = k_e r_e \arcsin\left( \frac{r \cos\theta}{k_e r_e + h_n(r, \theta, r_e, k_e)}\right)\]

where \(h_n\) would be provided by bin_altitude.

Parameters:

r (numpy.ndarray) – Array of ranges [m]
theta (scalar or numpy.ndarray) – Array broadcastable to the shape of r elevation angles in degrees with 0° at horizontal and +90° pointing vertically upwards from the radar
binalt (numpy.ndarray) – site altitude [m] amsl. same shape as r.
re (float, optional) – earth’s radius [m], defaults to 6371000.
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths

Returns:

distance (numpy.ndarray) – Array of great circle arc distances [m]

wradlib.georef.snap_bounds(bounds, resolution)[source]#

Adjusts integer bounds so they align with the resolution grid.

This ensures that the width and height of the bounds are evenly divisible by the resolution, centering the snapped bounds around the original center.

Parameters:

bounds (tuple of int) – (minx, miny, maxx, maxy) in integer units.
resolution (int or tuple of int) – Desired resolution per axis.

Returns:

tuple of int – Snapped bounds aligned to the resolution grid.

wradlib.georef.snap_resolution(extent, target)[source]#

Snap an integer resolution so it evenly divides the extent, staying close to the target.

This function finds the closest integer resolution that divides the given extent exactly, while minimizing the difference from the desired target resolution.

Parameters:

extent (int) – Total extent along one axis (e.g., width or height in arcseconds or pixels).
target (int) – Desired resolution (size of each cell), must be a positive integer.

Returns:

int – Snapped resolution that evenly divides the extent and is closest to the target.

wradlib.georef.spherical_to_centroids(r, phi, theta, site, *, crs=None)[source]#

wradlib.georef.spherical_to_centroids(obj: Dataset, **kwargs)

Generate 3-D centroids of the radar bins from the sperical coordinates (r, phi, theta).

For further information refer to the documentation of spherical_to_xyz.

Parameters:

r (numpy.ndarray) – Array of ranges [m]; r defines the exterior boundaries of the range bins! (not the centroids). Thus, values must be positive!
phi (numpy.ndarray) – Array of azimuth angles containing values between 0° and 360°. The angles are assumed to describe the pointing direction fo the main beam lobe! The first angle can start at any values, but make sure the array is sorted continuously positively clockwise and the angles are equidistant. An angle if 0 degree is pointing north.
theta (float) – Elevation angle of scan
site (sequence) – the lon/lat/alt coordinates of the radar location
crs – Coordinate Reference System (CRS). Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.

Returns:

output (numpy.ndarray) – A 3-d array of bin centroids with shape(num_rays, num_bins, 3). The last dimension carries the xyz-coordinates either in aeqd or given crs.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD, only returned if crs is None.

Note

Azimuth angles of 360 deg are internally converted to 0 deg.

wradlib.georef.spherical_to_polyvert(r, phi, theta, site, *, crs=None)[source]#

wradlib.georef.spherical_to_polyvert(obj: DataArray, **kwargs)

wradlib.georef.spherical_to_polyvert(obj: Dataset, **kwargs)

Generate 3-D polygon vertices directly from spherical coordinates (r, phi, theta).

This is an alternative to centroid_to_polyvert which does not use centroids, but generates the polygon vertices by simply connecting the corners of the radar bins.

Both azimuth and range arrays are assumed to be equidistant and to contain only unique values. For further information refer to the documentation of spherical_to_xyz.

Parameters:

r (numpy.ndarray) – Array of ranges [m]; r defines the exterior boundaries of the range bins! (not the centroids). Thus, values must be positive!
phi (numpy.ndarray) – Array of azimuth angles containing values between 0° and 360°. The angles are assumed to describe the pointing direction fo the main beam lobe! The first angle can start at any values, but make sure the array is sorted continuously positively clockwise and the angles are equidistant. An angle if 0 degree is pointing north.
theta (float) – Elevation angle of scan
site (sequence) – the lon/lat/alt coordinates of the radar location
crs – Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj:pyproj.CRS instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to AEQD.

Returns:

output (numpy.ndarray) – A 3-d array of polygon vertices with shape(num_vertices, num_vertex_nodes, 2). The last dimension carries the xyz-coordinates either in aeqd or given crs.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD, only returned if crs is None.

Examples

>>> import wradlib.georef as georef  # noqa
>>> import numpy as np
>>> from matplotlib import collections
>>> import matplotlib.pyplot as plt
>>> # define the polar coordinates and the site coordinates in lat/lon
>>> r = np.array([50., 100., 150., 200.]) * 1000
>>> az = np.array([0., 45., 90., 135., 180., 225., 270., 315.])
>>> el = 1.0
>>> site = (9.0, 48.0, 0)
>>> polygons, aeqd = georef.spherical_to_polyvert(r, az, el, site)
>>> # plot the resulting mesh
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> polycoll = collections.PolyCollection(polygons[...,:2], closed=True, facecolors='None')  # noqa
>>> ret = ax.add_collection(polycoll, autolim=True)
>>> plt.autoscale()
>>> plt.show()

wradlib.georef.spherical_to_proj(r, phi, theta, site, *, crs=None, re=None, ke=1.3333333333333333)[source]#

wradlib.georef.spherical_to_proj(obj: DataArray, **kwargs)

wradlib.georef.spherical_to_proj(obj: Dataset, **kwargs)

Transforms spherical coordinates (r, phi, theta) to projected coordinates centered at site in given projection.

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

r (numpy.ndarray) – Contains the radial distances.
phi (numpy.ndarray) – Contains the azimuthal angles.
theta (numpy.ndarray) – Contains the elevation angles.
site (sequence) – the lon / lat coordinates of the radar location and its altitude a.m.s.l. (in meters) if site is of length two, altitude is assumed to be zero
crs – Coordinate Reference System (CRS) of the coordinates. Can be one of:
- A pyproj.crs.CoordinateSystem instance
- A cartopy.crs.CRS instance
- A osgeo.osr.SpatialReference instance
- A type accepted by pyproj.crs.CRS.from_user_input (e.g., EPSG code, PROJ string, dictionary, WKT, or any object with a to_wkt() method)
Defaults to EPSG(4326).
re (float, optional) – earth’s radius [m], defaults to None (calculating from given latitude).
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Returns:

coords (numpy.ndarray) – Array of shape (…, 3). Contains projected map coordinates.

Examples

A few standard directions (North, South, North, East, South, West) with different distances (amounting to roughly 1°) from a site located at 48°N 9°E

>>> r  = np.array([0.,   0., 111., 111., 111., 111.,])*1000
>>> az = np.array([0., 180.,   0.,  90., 180., 270.,])
>>> th = np.array([0.,   0.,   0.,   0.,   0.,  0.5,])
>>> csite = (9.0, 48.0, 0)
>>> coords = spherical_to_proj(r, az, th, csite)
>>> for coord in coords:
...     print( '{0:7.4f}, {1:7.4f}, {2:7.4f}'.format(*coord))
...
 9.0000, 48.0000,  0.0000
 9.0000, 48.0000,  0.0000
 9.0000, 48.9981, 725.7160
10.4872, 47.9904, 725.7160
 9.0000, 47.0017, 725.7160
 7.5131, 47.9904, 1694.2234

Here, the coordinates of the east and west directions won’t come to lie on the latitude of the site because the beam doesn’t travel along the latitude circle but along a great circle.

See wradlib_workflow:georeferencing-and-projection.

wradlib.georef.spherical_to_xyz(r, phi, theta, site, *, re=None, ke=1.3333333333333333, **kwargs)[source]#

wradlib.georef.spherical_to_xyz(obj: DataArray, **kwargs)

wradlib.georef.spherical_to_xyz(obj: Dataset, **kwargs)

Transforms spherical coordinates (r, phi, theta) to cartesian coordinates (x, y, z) centered at site (aeqd).

It takes the shortening of the great circle distance with increasing elevation angle as well as the resulting increase in height into account.

Parameters:

r (numpy.ndarray) – Contains the radial distances in meters.
phi (numpy.ndarray) – Contains the azimuthal angles in degree.
theta (numpy.ndarray) – Contains the elevation angles in degree.
site (sequence) – the lon / lat / alt coordinates of the radar location and its altitude a.m.s.l. (in meters)
re (float, optional) – earth’s radius [m], defaults to None (calculating from given latitude).
ke (float, optional) – adjustment factor to account for the refractivity gradient that affects radar beam propagation. In principle this is wavelength- dependend. The default of 4/3 is a good approximation for most weather radar wavelengths.

Keyword Arguments:

squeeze (bool, optional) – If True, returns squeezed array. Defaults to False.
strict_dims (bool, optional) – If True, generates output of (theta, phi, r, 3) in any case. If False, dimensions with same length are “merged”. Defaults to False.

Returns:

xyz (numpy.ndarray) – Array of shape (…, 3). Contains cartesian coordinates.
aeqd (pyproj.crs.CoordinateSystem) – Coordinate Reference System (CRS), AEQD.

wradlib.georef.sweep_centroids(nrays, rscale, nbins, elangle, /)[source]#

Construct sweep centroids native coordinates.

Parameters:

nrays (int) – number of rays
rscale (float) – length [m] of a range bin
nbins (int) – number of range bins
elangle (float) – elevation angle [deg]

Returns:

coordinates (numpy.ndarray) – array of shape (nrays,nbins,3) containing native centroid radar coordinates (slant range, azimuth, elevation)

wradlib.georef.transform_geometry(geom, trg_crs, **kwargs)[source]#

Perform geotransformation to given destination SpatialReferenceSystem

It transforms coordinates to a given destination osr spatial reference if a geotransform is neccessary.

Parameters:

geom (osgeo.ogr.Geometry)
trg_crs (pyproj:pyproj.CRS, cartopy.crs.CRS, osgeo.osr.SpatialReference, str or int) – Any of the given projection objects, or anything which can be consumed by pyproj.CRS.from_user_input Destination Projection

Keyword Arguments:

Returns:

geom (osgeo.ogr.Geometry) – Transformed Geometry

wradlib.georef.wkt_to_osr(wkt=None)[source]#

Create osr spatial reference object from WKT string

Parameters:: wkt (str) – WTK string defining the coordinate reference system
Returns:: crs (osgeo.osr.SpatialReference) – GDAL/OSR object defining projection

Deprecated since version 2.4: This will be removed in 3.0. Use pyproj.CRS.from_wkt instead.

wradlib.georef.xyz_to_spherical(*args, **kwargs)[source]#
wradlib.georef.xyz_to_spherical(xyz: ndarray, *, altitude=0, crs=None, ke=1.3333333333333333)
wradlib.georef.xyz_to_spherical(obj: DataArray, **kwargs)
wradlib.georef.xyz_to_spherical(obj: Dataset, **kwargs)

Georeferencing

Contents

Georeferencing#

Miscellaneous#

Polar Grid Functions#

Projection Functions#

Raster Functions#

Rectangular Grid Functions#

Satellite Functions#

Vector Functions (GDAL)#

Xarray Functions#