# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Definitions of coordinate systems."""
from abc import ABCMeta, abstractmethod
from functools import cached_property
import re
from typing import ClassVar
import warnings
import cartopy.crs as ccrs
import numpy as np
from iris._deprecation import warn_deprecated
import iris.warnings
def _arg_default(value, default, cast_as=float):
"""Apply a default value and type for an optional kwarg."""
if value is None:
value = default
value = cast_as(value)
return value
def _1or2_parallels(arg):
"""Accept 1 or 2 inputs as a tuple of 1 or 2 floats."""
try:
values_tuple = tuple(arg)
except TypeError:
values_tuple = (arg,)
values_tuple = tuple([float(x) for x in values_tuple])
nvals = len(values_tuple)
if nvals not in (1, 2):
emsg = "Allows only 1 or 2 parallels or secant latitudes : got {!r}"
raise ValueError(emsg.format(arg))
return values_tuple
def _float_or_None(arg):
"""Cast as float, except for allowing None as a distinct valid value."""
if arg is not None:
arg = float(arg)
return arg
[docs]
class CoordSystem(metaclass=ABCMeta):
"""Abstract base class for coordinate systems."""
grid_mapping_name: ClassVar[str | None] = None
[docs]
def __eq__(self, other):
"""Override equality.
The `_globe` and `_crs` attributes are not compared because they are
cached properties and completely derived from other attributes. The
nature of caching means that they can appear on one object and not on
another despite the objects being identical, and them being completely
derived from other attributes means they will only differ if other
attributes that are being tested for equality differ.
"""
if self.__class__ != other.__class__:
return False
self_keys = set(self.__dict__.keys())
other_keys = set(other.__dict__.keys())
check_keys = (self_keys | other_keys) - {"_globe", "_crs"}
for key in check_keys:
try:
if self.__dict__[key] != other.__dict__[key]:
return False
except KeyError:
return False
return True
def __ne__(self, other):
# Must supply __ne__, Python does not defer to __eq__ for
# negative equality.
return not (self == other)
[docs]
def xml_element(self, doc, attrs=None):
"""Perform default behaviour for coord systems."""
# attrs - optional list of (k,v) items, used for alternate output
xml_element_name = type(self).__name__
# lower case the first char
first_char = xml_element_name[0]
xml_element_name = xml_element_name.replace(first_char, first_char.lower(), 1)
coord_system_xml_element = doc.createElement(xml_element_name)
if attrs is None:
attrs = self.__dict__.items()
attrs = sorted(attrs, key=lambda attr: attr[0])
for name, value in attrs:
if isinstance(value, float):
value_str = "{:.16}".format(value)
elif isinstance(value, np.float32):
value_str = "{:.8}".format(value)
else:
value_str = "{}".format(value)
coord_system_xml_element.setAttribute(name, value_str)
return coord_system_xml_element
@staticmethod
def _ellipsoid_to_globe(ellipsoid, globe_default):
if ellipsoid is not None:
globe = ellipsoid.as_cartopy_globe()
else:
globe = globe_default
return globe
[docs]
@abstractmethod
def as_cartopy_crs(self):
"""Return a cartopy CRS representing our native coordinate system."""
pass
[docs]
@abstractmethod
def as_cartopy_projection(self):
"""Return a cartopy projection representing our native map.
This will be the same as the :func:`~CoordSystem.as_cartopy_crs` for
map projections but for spherical coord systems (which are not map
projections) we use a map projection, such as PlateCarree.
"""
pass
_short_datum_names = {
"OSGB 1936": "OSGB36",
"OSGB_1936": "OSGB36",
"WGS 84": "WGS84",
}
[docs]
class GeogCS(CoordSystem):
"""A geographic (ellipsoidal) coordinate system.
A geographic (ellipsoidal) coordinate system, defined by the shape of
the Earth and a prime meridian.
"""
grid_mapping_name = "latitude_longitude"
def __init__(
self,
semi_major_axis=None,
semi_minor_axis=None,
inverse_flattening=None,
longitude_of_prime_meridian=None,
):
"""Create a new GeogCS.
Parameters
----------
semi_major_axis, semi_minor_axis : optional
Axes of ellipsoid, in metres. At least one must be given (see note
below).
inverse_flattening : optional
Can be omitted if both axes given (see note below). Default 0.0.
longitude_of_prime_meridian : optional
Specifies the prime meridian on the ellipsoid, in degrees. Default 0.0.
Notes
-----
If just semi_major_axis is set, with no semi_minor_axis or
inverse_flattening, then a perfect sphere is created from the given
radius.
If just two of semi_major_axis, semi_minor_axis, and inverse_flattening
are given the missing element is calculated from the formula:
:math:`flattening = (major - minor) / major`
Currently, Iris will not allow over-specification (all three ellipsoid
parameters).
After object creation, altering any of these properties will not update
the others. semi_major_axis and semi_minor_axis are used when creating
Cartopy objects.
Examples::
cs = GeogCS(6371229)
pp_cs = GeogCS(iris.fileformats.pp.EARTH_RADIUS)
airy1830 = GeogCS(semi_major_axis=6377563.396,
semi_minor_axis=6356256.909)
airy1830 = GeogCS(semi_major_axis=6377563.396,
inverse_flattening=299.3249646)
custom_cs = GeogCS(6400000, 6300000)
"""
# No ellipsoid specified? (0 0 0)
if (
(semi_major_axis is None)
and (semi_minor_axis is None)
and (inverse_flattening is None)
):
raise ValueError("No ellipsoid specified")
# Ellipsoid over-specified? (1 1 1)
if (
(semi_major_axis is not None)
and (semi_minor_axis is not None)
and (inverse_flattening is not None)
):
raise ValueError("Ellipsoid is overspecified")
# Perfect sphere (semi_major_axis only)? (1 0 0)
elif semi_major_axis is not None and (
semi_minor_axis is None and not inverse_flattening
):
semi_minor_axis = semi_major_axis
inverse_flattening = 0.0
# Calculate semi_major_axis? (0 1 1)
elif semi_major_axis is None and (
semi_minor_axis is not None and inverse_flattening is not None
):
semi_major_axis = -semi_minor_axis / (
(1.0 - inverse_flattening) / inverse_flattening
)
# Calculate semi_minor_axis? (1 0 1)
elif semi_minor_axis is None and (
semi_major_axis is not None and inverse_flattening is not None
):
semi_minor_axis = semi_major_axis - (
(1.0 / inverse_flattening) * semi_major_axis
)
# Calculate inverse_flattening? (1 1 0)
elif inverse_flattening is None and (
semi_major_axis is not None and semi_minor_axis is not None
):
if semi_major_axis == semi_minor_axis:
inverse_flattening = 0.0
else:
inverse_flattening = 1.0 / (
(semi_major_axis - semi_minor_axis) / semi_major_axis
)
# We didn't get enough to specify an ellipse.
else:
raise ValueError("Insufficient ellipsoid specification")
#: Major radius of the ellipsoid in metres.
self._semi_major_axis = float(semi_major_axis)
#: Minor radius of the ellipsoid in metres.
self._semi_minor_axis = float(semi_minor_axis)
#: :math:`1/f` where :math:`f = (a-b)/a`.
self._inverse_flattening = float(inverse_flattening)
self._datum = None
#: Describes 'zero' on the ellipsoid in degrees.
self.longitude_of_prime_meridian = _arg_default(longitude_of_prime_meridian, 0)
def _pretty_attrs(self):
attrs = [("semi_major_axis", self.semi_major_axis)]
if self.semi_major_axis != self.semi_minor_axis:
attrs.append(("semi_minor_axis", self.semi_minor_axis))
if self.longitude_of_prime_meridian != 0.0:
attrs.append(
(
"longitude_of_prime_meridian",
self.longitude_of_prime_meridian,
)
)
# An unknown crs datum will be treated as None
if self.datum is not None and self.datum != "unknown":
attrs.append(
(
"datum",
self.datum,
)
)
return attrs
def __repr__(self):
attrs = self._pretty_attrs()
# Special case for 1 pretty attr
if len(attrs) == 1 and attrs[0][0] == "semi_major_axis":
return "GeogCS(%r)" % self.semi_major_axis
else:
return "GeogCS(%s)" % ", ".join(["%s=%r" % (k, v) for k, v in attrs])
def __str__(self):
attrs = self._pretty_attrs()
# Special case for 1 pretty attr
if len(attrs) == 1 and attrs[0][0] == "semi_major_axis":
return "GeogCS({:.16})".format(self.semi_major_axis)
else:
text_attrs = []
for k, v in attrs:
if isinstance(v, float):
text_attrs.append("{}={:.16}".format(k, v))
elif isinstance(v, np.float32):
text_attrs.append("{}={:.8}".format(k, v))
else:
text_attrs.append("{}={}".format(k, v))
return "GeogCS({})".format(", ".join(text_attrs))
[docs]
def xml_element(self, doc):
# Special output for spheres
attrs = self._pretty_attrs()
if len(attrs) == 1 and attrs[0][0] == "semi_major_axis":
attrs = [("earth_radius", self.semi_major_axis)]
return CoordSystem.xml_element(self, doc, attrs)
[docs]
def as_cartopy_crs(self):
return self._crs
[docs]
def as_cartopy_projection(self):
return ccrs.PlateCarree(
central_longitude=self.longitude_of_prime_meridian,
globe=self.as_cartopy_globe(),
)
[docs]
def as_cartopy_globe(self):
return self._globe
@cached_property
def _globe(self):
"""A representation of this CRS as a Cartopy Globe.
Note
----
This property is created when required and then cached for speed. That
cached value is cleared when an assignment is made to a property of the
class that invalidates the cache.
"""
if self._datum is not None:
short_datum = _short_datum_names.get(self._datum, self._datum)
# Cartopy doesn't actually enact datums unless they're provided without
# ellipsoid axes, so only provide the datum
return ccrs.Globe(short_datum, ellipse=None)
else:
return ccrs.Globe(
ellipse=None,
semimajor_axis=self._semi_major_axis,
semiminor_axis=self._semi_minor_axis,
)
@cached_property
def _crs(self):
"""A representation of this CRS as a Cartopy CRS.
Note
----
This property is created when required and then cached for speed. That
cached value is cleared when an assignment is made to a property of the
class that invalidates the cache.
"""
return ccrs.Geodetic(self._globe)
def _wipe_cached_properties(self):
"""Wipe the cached properties on the object.
Wipe the cached properties on the object as part of any update to a
value that invalidates the cache.
"""
try:
delattr(self, "_crs")
except AttributeError:
pass
try:
delattr(self, "_globe")
except AttributeError:
pass
@property
def semi_major_axis(self):
if self._semi_major_axis is not None:
return self._semi_major_axis
else:
return self._crs.ellipsoid.semi_major_metre
@semi_major_axis.setter
def semi_major_axis(self, value):
"""Assign semi_major_axis.
Setting this property to a different value invalidates the current datum
(if any) because a datum encodes a specific semi-major axis. This also
invalidates the cached `cartopy.Globe` and `cartopy.CRS`.
"""
value = float(value)
if not np.isclose(self.semi_major_axis, value):
self._datum = None
self._wipe_cached_properties()
self._semi_major_axis = value
@property
def semi_minor_axis(self):
if self._semi_minor_axis is not None:
return self._semi_minor_axis
else:
return self._crs.ellipsoid.semi_minor_metre
@semi_minor_axis.setter
def semi_minor_axis(self, value):
"""Assign semi_minor_axis.
Setting this property to a different value invalidates the current datum
(if any) because a datum encodes a specific semi-minor axis. This also
invalidates the cached `cartopy.Globe` and `cartopy.CRS`.
"""
value = float(value)
if not np.isclose(self.semi_minor_axis, value):
self._datum = None
self._wipe_cached_properties()
self._semi_minor_axis = value
@property
def inverse_flattening(self):
if self._inverse_flattening is not None:
return self._inverse_flattening
else:
self._crs.ellipsoid.inverse_flattening
@inverse_flattening.setter
def inverse_flattening(self, value):
"""Assign inverse_flattening.
Setting this property to a different value does not affect the behaviour
of this object any further than the value of this property.
"""
wmsg = (
"Setting inverse_flattening does not affect other properties of "
"the GeogCS object. To change other properties set them explicitly"
" or create a new GeogCS instance."
)
warnings.warn(wmsg, category=iris.warnings.IrisUserWarning)
value = float(value)
self._inverse_flattening = value
@property
def datum(self):
if self._datum is None:
return None
else:
datum = self._datum
return datum
@datum.setter
def datum(self, value):
"""Assign datum.
Setting this property to a different value invalidates the current
values of the ellipsoid measurements because a datum encodes its own
ellipse. This also invalidates the cached `cartopy.Globe` and
`cartopy.CRS`.
"""
if self._datum != value:
self._semi_major_axis = None
self._semi_minor_axis = None
self._inverse_flattening = None
self._wipe_cached_properties()
self._datum = value
[docs]
@classmethod
def from_datum(cls, datum, longitude_of_prime_meridian=None):
crs = super().__new__(cls)
crs._semi_major_axis = None
crs._semi_minor_axis = None
crs._inverse_flattening = None
#: Describes 'zero' on the ellipsoid in degrees.
crs.longitude_of_prime_meridian = _arg_default(longitude_of_prime_meridian, 0)
crs._datum = datum
return crs
[docs]
class RotatedGeogCS(CoordSystem):
"""A coordinate system with rotated pole, on an optional :class:`GeogCS`."""
grid_mapping_name = "rotated_latitude_longitude"
def __init__(
self,
grid_north_pole_latitude,
grid_north_pole_longitude,
north_pole_grid_longitude=None,
ellipsoid=None,
):
"""Construct a coordinate system with rotated pole, on an optional :class:`GeogCS`.
Parameters
----------
grid_north_pole_latitude :
The true latitude of the rotated pole in degrees.
grid_north_pole_longitude :
The true longitude of the rotated pole in degrees.
north_pole_grid_longitude : optional
Longitude of true north pole in rotated grid, in degrees.
Defaults to 0.0.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
Examples
--------
::
rotated_cs = RotatedGeogCS(30, 30)
another_cs = RotatedGeogCS(30, 30,
ellipsoid=GeogCS(6400000, 6300000))
"""
#: The true latitude of the rotated pole in degrees.
self.grid_north_pole_latitude = float(grid_north_pole_latitude)
#: The true longitude of the rotated pole in degrees.
self.grid_north_pole_longitude = float(grid_north_pole_longitude)
#: Longitude of true north pole in rotated grid in degrees.
self.north_pole_grid_longitude = _arg_default(north_pole_grid_longitude, 0)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def _pretty_attrs(self):
attrs = [
("grid_north_pole_latitude", self.grid_north_pole_latitude),
("grid_north_pole_longitude", self.grid_north_pole_longitude),
]
if self.north_pole_grid_longitude != 0.0:
attrs.append(("north_pole_grid_longitude", self.north_pole_grid_longitude))
if self.ellipsoid is not None:
attrs.append(("ellipsoid", self.ellipsoid))
return attrs
def __repr__(self):
attrs = self._pretty_attrs()
result = "RotatedGeogCS(%s)" % ", ".join(["%s=%r" % (k, v) for k, v in attrs])
# Extra prettiness
result = result.replace("grid_north_pole_latitude=", "")
result = result.replace("grid_north_pole_longitude=", "")
return result
def __str__(self):
attrs = self._pretty_attrs()
text_attrs = []
for k, v in attrs:
if isinstance(v, float):
text_attrs.append("{}={:.16}".format(k, v))
elif isinstance(v, np.float32):
text_attrs.append("{}={:.8}".format(k, v))
else:
text_attrs.append("{}={}".format(k, v))
result = "RotatedGeogCS({})".format(", ".join(text_attrs))
# Extra prettiness
result = result.replace("grid_north_pole_latitude=", "")
result = result.replace("grid_north_pole_longitude=", "")
return result
[docs]
def xml_element(self, doc):
return CoordSystem.xml_element(self, doc, self._pretty_attrs())
def _ccrs_kwargs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, None)
cartopy_kwargs = {
"central_rotated_longitude": self.north_pole_grid_longitude,
"pole_longitude": self.grid_north_pole_longitude,
"pole_latitude": self.grid_north_pole_latitude,
"globe": globe,
}
return cartopy_kwargs
[docs]
def as_cartopy_crs(self):
return ccrs.RotatedGeodetic(**self._ccrs_kwargs())
[docs]
def as_cartopy_projection(self):
return ccrs.RotatedPole(**self._ccrs_kwargs())
[docs]
class TransverseMercator(CoordSystem):
"""A cylindrical map projection, with XY coordinates measured in metres."""
grid_mapping_name = "transverse_mercator"
def __init__(
self,
latitude_of_projection_origin,
longitude_of_central_meridian,
false_easting=None,
false_northing=None,
scale_factor_at_central_meridian=None,
ellipsoid=None,
):
"""Construct a TransverseMercator object.
Parameters
----------
latitude_of_projection_origin :
True latitude of planar origin in degrees.
longitude_of_central_meridian :
True longitude of planar origin in degrees.
false_easting : optional
X offset from planar origin in metres.
Defaults to 0.0.
false_northing : optional
Y offset from planar origin in metres.
Defaults to 0.0.
scale_factor_at_central_meridian : optional
Reduces the cylinder to slice through the ellipsoid
(secant form). Used to provide TWO longitudes of zero
distortion in the area of interest.
Defaults to 1.0 .
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
Examples
--------
::
airy1830 = GeogCS(6377563.396, 6356256.909)
osgb = TransverseMercator(49, -2, 400000, -100000, 0.9996012717,
ellipsoid=airy1830)
"""
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = float(latitude_of_projection_origin)
#: True longitude of planar origin in degrees.
self.longitude_of_central_meridian = float(longitude_of_central_meridian)
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: Scale factor at the centre longitude.
self.scale_factor_at_central_meridian = _arg_default(
scale_factor_at_central_meridian, 1.0
)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"TransverseMercator(latitude_of_projection_origin={!r}, "
"longitude_of_central_meridian={!r}, false_easting={!r}, "
"false_northing={!r}, scale_factor_at_central_meridian={!r}, "
"ellipsoid={!r})".format(
self.latitude_of_projection_origin,
self.longitude_of_central_meridian,
self.false_easting,
self.false_northing,
self.scale_factor_at_central_meridian,
self.ellipsoid,
)
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, None)
return ccrs.TransverseMercator(
central_longitude=self.longitude_of_central_meridian,
central_latitude=self.latitude_of_projection_origin,
false_easting=self.false_easting,
false_northing=self.false_northing,
scale_factor=self.scale_factor_at_central_meridian,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class OSGB(TransverseMercator):
"""A Specific transverse mercator projection on a specific ellipsoid."""
def __init__(self):
TransverseMercator.__init__(
self,
49,
-2,
400000,
-100000,
0.9996012717,
GeogCS(6377563.396, 6356256.909),
)
[docs]
def as_cartopy_crs(self):
return ccrs.OSGB()
[docs]
def as_cartopy_projection(self):
return ccrs.OSGB()
[docs]
class Orthographic(CoordSystem):
"""An orthographic map projection."""
grid_mapping_name = "orthographic"
def __init__(
self,
latitude_of_projection_origin,
longitude_of_projection_origin,
false_easting=None,
false_northing=None,
ellipsoid=None,
):
"""Construct an Orthographic coord system.
Parameters
----------
latitude_of_projection_origin :
True latitude of planar origin in degrees.
longitude_of_projection_origin :
True longitude of planar origin in degrees.
false_easting : optional
X offset from planar origin in metres. Defaults to 0.0.
false_northing : optional
Y offset from planar origin in metres. Defaults to 0.0.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
"""
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = float(latitude_of_projection_origin)
#: True longitude of planar origin in degrees.
self.longitude_of_projection_origin = float(longitude_of_projection_origin)
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"Orthographic(latitude_of_projection_origin={!r}, "
"longitude_of_projection_origin={!r}, "
"false_easting={!r}, false_northing={!r}, "
"ellipsoid={!r})".format(
self.latitude_of_projection_origin,
self.longitude_of_projection_origin,
self.false_easting,
self.false_northing,
self.ellipsoid,
)
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
warnings.warn(
"Discarding false_easting and false_northing that are not used by Cartopy.",
category=iris.warnings.IrisDefaultingWarning,
)
return ccrs.Orthographic(
central_longitude=self.longitude_of_projection_origin,
central_latitude=self.latitude_of_projection_origin,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class VerticalPerspective(CoordSystem):
"""A vertical/near-side perspective satellite image map projection."""
grid_mapping_name = "vertical_perspective"
def __init__(
self,
latitude_of_projection_origin,
longitude_of_projection_origin,
perspective_point_height,
false_easting=None,
false_northing=None,
ellipsoid=None,
):
"""Construct a Vertical Perspective coord system.
Parameters
----------
latitude_of_projection_origin :
True latitude of planar origin in degrees.
longitude_of_projection_origin :
True longitude of planar origin in degrees.
perspective_point_height :
Altitude of satellite in metres above the surface of the
ellipsoid.
false_easting : optional
X offset from planar origin in metres. Defaults to 0.0.
false_northing : optional
Y offset from planar origin in metres. Defaults to 0.0.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
"""
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = float(latitude_of_projection_origin)
#: True longitude of planar origin in degrees.
self.longitude_of_projection_origin = float(longitude_of_projection_origin)
#: Altitude of satellite in metres.
self.perspective_point_height = float(perspective_point_height)
# TODO: test if may be cast to float for proj.4
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"Vertical Perspective(latitude_of_projection_origin={!r}, "
"longitude_of_projection_origin={!r}, "
"perspective_point_height={!r}, "
"false_easting={!r}, false_northing={!r}, "
"ellipsoid={!r})".format(
self.latitude_of_projection_origin,
self.longitude_of_projection_origin,
self.perspective_point_height,
self.false_easting,
self.false_northing,
self.ellipsoid,
)
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
return ccrs.NearsidePerspective(
central_latitude=self.latitude_of_projection_origin,
central_longitude=self.longitude_of_projection_origin,
satellite_height=self.perspective_point_height,
false_easting=self.false_easting,
false_northing=self.false_northing,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class Geostationary(CoordSystem):
"""A geostationary satellite image map projection."""
grid_mapping_name = "geostationary"
def __init__(
self,
latitude_of_projection_origin,
longitude_of_projection_origin,
perspective_point_height,
sweep_angle_axis,
false_easting=None,
false_northing=None,
ellipsoid=None,
):
"""Construct a Geostationary coord system.
Parameters
----------
latitude_of_projection_origin :
True latitude of planar origin in degrees.
longitude_of_projection_origin :
True longitude of planar origin in degrees.
perspective_point_height :
Altitude of satellite in metres above the surface of the ellipsoid.
sweep_angle_axis : str
The axis along which the satellite instrument sweeps - 'x' or 'y'.
false_easting : optional
X offset from planar origin in metres. Defaults to 0.0.
false_northing : optional
Y offset from planar origin in metres. Defaults to 0.0.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
"""
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = float(latitude_of_projection_origin)
if self.latitude_of_projection_origin != 0.0:
raise ValueError(
"Non-zero latitude of projection currently not supported by Cartopy."
)
#: True longitude of planar origin in degrees.
self.longitude_of_projection_origin = float(longitude_of_projection_origin)
#: Altitude of satellite in metres.
self.perspective_point_height = float(perspective_point_height)
# TODO: test if may be cast to float for proj.4
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: The sweep angle axis (string 'x' or 'y').
self.sweep_angle_axis = sweep_angle_axis
if self.sweep_angle_axis not in ("x", "y"):
raise ValueError('Invalid sweep_angle_axis - must be "x" or "y"')
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"Geostationary(latitude_of_projection_origin={!r}, "
"longitude_of_projection_origin={!r}, "
"perspective_point_height={!r}, false_easting={!r}, "
"false_northing={!r}, sweep_angle_axis={!r}, "
"ellipsoid={!r}".format(
self.latitude_of_projection_origin,
self.longitude_of_projection_origin,
self.perspective_point_height,
self.false_easting,
self.false_northing,
self.sweep_angle_axis,
self.ellipsoid,
)
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
return ccrs.Geostationary(
central_longitude=self.longitude_of_projection_origin,
satellite_height=self.perspective_point_height,
false_easting=self.false_easting,
false_northing=self.false_northing,
globe=globe,
sweep_axis=self.sweep_angle_axis,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class Stereographic(CoordSystem):
"""A stereographic map projection."""
grid_mapping_name = "stereographic"
def __init__(
self,
central_lat,
central_lon,
false_easting=None,
false_northing=None,
true_scale_lat=None,
ellipsoid=None,
scale_factor_at_projection_origin=None,
):
"""Construct a Stereographic coord system.
Parameters
----------
central_lat : float
The latitude of the pole.
central_lon : float
The central longitude, which aligns with the y axis.
false_easting : float, optional
X offset from planar origin in metres.
false_northing : float, optional
Y offset from planar origin in metres.
true_scale_lat : float, optional
Latitude of true scale.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
scale_factor_at_projection_origin : float, optional
Scale factor at the origin of the projection.
Notes
-----
It is only valid to provide one of true_scale_lat and
scale_factor_at_projection_origin
"""
#: True latitude of planar origin in degrees.
self.central_lat = float(central_lat)
#: True longitude of planar origin in degrees.
self.central_lon = float(central_lon)
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: Latitude of true scale.
self.true_scale_lat = _arg_default(true_scale_lat, None, cast_as=_float_or_None)
#: Scale factor at projection origin.
self.scale_factor_at_projection_origin = _arg_default(
scale_factor_at_projection_origin, None, cast_as=_float_or_None
)
# N.B. the way we use these parameters, we need them to default to None,
# and *not* to 0.0.
if (
self.true_scale_lat is not None
and self.scale_factor_at_projection_origin is not None
):
raise ValueError(
"It does not make sense to provide both "
'"scale_factor_at_projection_origin" and "true_scale_latitude". '
)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def _repr_attributes(self):
if self.scale_factor_at_projection_origin is None:
scale_info = "true_scale_lat={!r}, ".format(self.true_scale_lat)
else:
scale_info = "scale_factor_at_projection_origin={!r}, ".format(
self.scale_factor_at_projection_origin
)
return (
f"(central_lat={self.central_lat}, central_lon={self.central_lon}, "
f"false_easting={self.false_easting}, false_northing={self.false_northing}, "
f"{scale_info}"
f"ellipsoid={self.ellipsoid})"
)
def __repr__(self):
return "Stereographic" + self._repr_attributes()
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
return ccrs.Stereographic(
self.central_lat,
self.central_lon,
self.false_easting,
self.false_northing,
self.true_scale_lat,
self.scale_factor_at_projection_origin,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class PolarStereographic(Stereographic):
"""A subclass of the stereographic map projection centred on a pole."""
grid_mapping_name = "polar_stereographic"
def __init__(
self,
central_lat,
central_lon,
false_easting=None,
false_northing=None,
true_scale_lat=None,
scale_factor_at_projection_origin=None,
ellipsoid=None,
):
"""Construct a Polar Stereographic coord system.
Parameters
----------
central_lat : {90, -90}
The latitude of the pole.
central_lon : float
The central longitude, which aligns with the y axis.
false_easting : float, optional
X offset from planar origin in metres.
false_northing : float, optional
Y offset from planar origin in metres.
true_scale_lat : float, optional
Latitude of true scale.
scale_factor_at_projection_origin : float, optional
Scale factor at the origin of the projection.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
Notes
-----
It is only valid to provide at most one of `true_scale_lat` and
`scale_factor_at_projection_origin`.
"""
super().__init__(
central_lat=central_lat,
central_lon=central_lon,
false_easting=false_easting,
false_northing=false_northing,
true_scale_lat=true_scale_lat,
scale_factor_at_projection_origin=scale_factor_at_projection_origin,
ellipsoid=ellipsoid,
)
def __repr__(self):
return "PolarStereographic" + self._repr_attributes()
[docs]
class Mercator(CoordSystem):
"""A coordinate system in the Mercator projection."""
grid_mapping_name = "mercator"
def __init__(
self,
longitude_of_projection_origin=None,
ellipsoid=None,
standard_parallel=None,
scale_factor_at_projection_origin=None,
false_easting=None,
false_northing=None,
):
"""Construct a Mercator coord system.
Parameters
----------
longitude_of_projection_origin : optional
True longitude of planar origin in degrees. Defaults to 0.0.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
standard_parallel : optional
The latitude where the scale is 1. Defaults to 0.0.
scale_factor_at_projection_origin : optional
Scale factor at natural origin. Defaults to unused.
false_easting : optional
X offset from the planar origin in metres. Defaults to 0.0.
false_northing : optional
Y offset from the planar origin in metres. Defaults to 0.0.
datum : optional
If given, specifies the datumof the coordinate system. Only
respected if iris.Future.daum_support is set.
Notes
-----
Only one of ``standard_parallel`` and
``scale_factor_at_projection_origin`` should be included.
"""
#: True longitude of planar origin in degrees.
self.longitude_of_projection_origin = _arg_default(
longitude_of_projection_origin, 0
)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
# Initialise to None, then set based on arguments
#: The latitude where the scale is 1.
self.standard_parallel = None
# The scale factor at the origin of the projection
self.scale_factor_at_projection_origin = None
if scale_factor_at_projection_origin is None:
self.standard_parallel = _arg_default(standard_parallel, 0)
else:
if standard_parallel is None:
self.scale_factor_at_projection_origin = _arg_default(
scale_factor_at_projection_origin, 0
)
else:
raise ValueError(
"It does not make sense to provide both "
'"scale_factor_at_projection_origin" and '
'"standard_parallel".'
)
#: X offset from the planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from the planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
def __repr__(self):
res = (
"Mercator(longitude_of_projection_origin="
"{self.longitude_of_projection_origin!r}, "
"ellipsoid={self.ellipsoid!r}, "
"standard_parallel={self.standard_parallel!r}, "
"scale_factor_at_projection_origin="
"{self.scale_factor_at_projection_origin!r}, "
"false_easting={self.false_easting!r}, "
"false_northing={self.false_northing!r})"
)
return res.format(self=self)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
return ccrs.Mercator(
central_longitude=self.longitude_of_projection_origin,
globe=globe,
latitude_true_scale=self.standard_parallel,
scale_factor=self.scale_factor_at_projection_origin,
false_easting=self.false_easting,
false_northing=self.false_northing,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class LambertAzimuthalEqualArea(CoordSystem):
"""A coordinate system in the Lambert Azimuthal Equal Area projection."""
grid_mapping_name = "lambert_azimuthal_equal_area"
def __init__(
self,
latitude_of_projection_origin=None,
longitude_of_projection_origin=None,
false_easting=None,
false_northing=None,
ellipsoid=None,
):
"""Construct a Lambert Azimuthal Equal Area coord system.
Parameters
----------
latitude_of_projection_origin : optional
True latitude of planar origin in degrees. Defaults to 0.0.
longitude_of_projection_origin : optional
True longitude of planar origin in degrees. Defaults to 0.0.
false_easting : optional
X offset from planar origin in metres. Defaults to 0.0.
false_northing : optional
Y offset from planar origin in metres. Defaults to 0.0.
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
"""
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = _arg_default(
latitude_of_projection_origin, 0
)
#: True longitude of planar origin in degrees.
self.longitude_of_projection_origin = _arg_default(
longitude_of_projection_origin, 0
)
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"LambertAzimuthalEqualArea(latitude_of_projection_origin={!r},"
" longitude_of_projection_origin={!r}, false_easting={!r},"
" false_northing={!r}, ellipsoid={!r})"
).format(
self.latitude_of_projection_origin,
self.longitude_of_projection_origin,
self.false_easting,
self.false_northing,
self.ellipsoid,
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
return ccrs.LambertAzimuthalEqualArea(
central_longitude=self.longitude_of_projection_origin,
central_latitude=self.latitude_of_projection_origin,
false_easting=self.false_easting,
false_northing=self.false_northing,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class AlbersEqualArea(CoordSystem):
"""A coordinate system in the Albers Conical Equal Area projection."""
grid_mapping_name = "albers_conical_equal_area"
def __init__(
self,
latitude_of_projection_origin=None,
longitude_of_central_meridian=None,
false_easting=None,
false_northing=None,
standard_parallels=None,
ellipsoid=None,
):
"""Construct a Albers Conical Equal Area coord system.
Parameters
----------
latitude_of_projection_origin : optional
True latitude of planar origin in degrees. Defaults to 0.0.
longitude_of_central_meridian : optional
True longitude of planar central meridian in degrees.
Defaults to 0.0.
false_easting : optional
X offset from planar origin in metres. Defaults to 0.0.
false_northing : optional
Y offset from planar origin in metres. Defaults to 0.0.
standard_parallels : number or iterable of 1 or 2 numbers, optional
The one or two latitudes of correct scale.
Defaults to (20.0, 50.0).
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
"""
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = _arg_default(
latitude_of_projection_origin, 0
)
#: True longitude of planar central meridian in degrees.
self.longitude_of_central_meridian = _arg_default(
longitude_of_central_meridian, 0
)
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: The one or two latitudes of correct scale (tuple of 1 or 2 floats).
self.standard_parallels = _arg_default(
standard_parallels, (20, 50), cast_as=_1or2_parallels
)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"AlbersEqualArea(latitude_of_projection_origin={!r},"
" longitude_of_central_meridian={!r}, false_easting={!r},"
" false_northing={!r}, standard_parallels={!r},"
" ellipsoid={!r})"
).format(
self.latitude_of_projection_origin,
self.longitude_of_central_meridian,
self.false_easting,
self.false_northing,
self.standard_parallels,
self.ellipsoid,
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
return ccrs.AlbersEqualArea(
central_longitude=self.longitude_of_central_meridian,
central_latitude=self.latitude_of_projection_origin,
false_easting=self.false_easting,
false_northing=self.false_northing,
standard_parallels=self.standard_parallels,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class ObliqueMercator(CoordSystem):
"""A cylindrical map projection, with XY coordinates measured in metres.
Designed for regions not well suited to :class:`Mercator` or
:class:`TransverseMercator`, as the positioning of the cylinder is more
customisable.
See Also
--------
RotatedMercator :
:class:`ObliqueMercator` with ``azimuth_of_central_line=90``.
"""
grid_mapping_name = "oblique_mercator"
def __init__(
self,
azimuth_of_central_line,
latitude_of_projection_origin,
longitude_of_projection_origin,
false_easting=None,
false_northing=None,
scale_factor_at_projection_origin=None,
ellipsoid=None,
):
"""Construct an ObliqueMercator object.
Parameters
----------
azimuth_of_central_line : float
Azimuth of centerline clockwise from north at the center point of
the centre line.
latitude_of_projection_origin : float
The true longitude of the central meridian in degrees.
longitude_of_projection_origin : float
The true latitude of the planar origin in degrees.
false_easting : float, optional
X offset from the planar origin in metres.
Defaults to 0.0.
false_northing : float, optional
Y offset from the planar origin in metres.
Defaults to 0.0.
scale_factor_at_projection_origin : float, optional
Scale factor at the central meridian.
Defaults to 1.0 .
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
Examples
--------
>>> from iris.coord_systems import GeogCS, ObliqueMercator
>>> my_ellipsoid = GeogCS(6371229.0, None, 0.0)
>>> ObliqueMercator(90.0, -22.0, -59.0, -25000.0, -25000.0, 1., my_ellipsoid)
ObliqueMercator(azimuth_of_central_line=90.0, latitude_of_projection_origin=-22.0, longitude_of_projection_origin=-59.0, false_easting=-25000.0, false_northing=-25000.0, scale_factor_at_projection_origin=1.0, ellipsoid=GeogCS(6371229.0))
"""
#: Azimuth of centerline clockwise from north.
self.azimuth_of_central_line = float(azimuth_of_central_line)
#: True latitude of planar origin in degrees.
self.latitude_of_projection_origin = float(latitude_of_projection_origin)
#: True longitude of planar origin in degrees.
self.longitude_of_projection_origin = float(longitude_of_projection_origin)
#: X offset from planar origin in metres.
self.false_easting = _arg_default(false_easting, 0)
#: Y offset from planar origin in metres.
self.false_northing = _arg_default(false_northing, 0)
#: Scale factor at the central meridian.
self.scale_factor_at_projection_origin = _arg_default(
scale_factor_at_projection_origin, 1.0
)
#: Ellipsoid definition (:class:`GeogCS` or None).
self.ellipsoid = ellipsoid
def __repr__(self):
return (
"{!s}(azimuth_of_central_line={!r}, "
"latitude_of_projection_origin={!r}, "
"longitude_of_projection_origin={!r}, false_easting={!r}, "
"false_northing={!r}, scale_factor_at_projection_origin={!r}, "
"ellipsoid={!r})".format(
self.__class__.__name__,
self.azimuth_of_central_line,
self.latitude_of_projection_origin,
self.longitude_of_projection_origin,
self.false_easting,
self.false_northing,
self.scale_factor_at_projection_origin,
self.ellipsoid,
)
)
[docs]
def as_cartopy_crs(self):
globe = self._ellipsoid_to_globe(self.ellipsoid, None)
return ccrs.ObliqueMercator(
central_longitude=self.longitude_of_projection_origin,
central_latitude=self.latitude_of_projection_origin,
false_easting=self.false_easting,
false_northing=self.false_northing,
scale_factor=self.scale_factor_at_projection_origin,
azimuth=self.azimuth_of_central_line,
globe=globe,
)
[docs]
def as_cartopy_projection(self):
return self.as_cartopy_crs()
[docs]
class RotatedMercator(ObliqueMercator):
""":class:`ObliqueMercator` with ``azimuth_of_central_line=90``.
As noted in CF versions 1.10 and earlier:
The Rotated Mercator projection is an Oblique Mercator projection
with azimuth = +90.
Notes
-----
.. deprecated:: 3.8.0
This coordinate system was introduced as already scheduled for removal
in a future release, since CF version 1.11 onwards now requires use of
:class:`ObliqueMercator` with ``azimuth_of_central_line=90.`` .
Any :class:`RotatedMercator` instances will always be saved to NetCDF
as the ``oblique_mercator`` grid mapping.
"""
def __init__(
self,
latitude_of_projection_origin,
longitude_of_projection_origin,
false_easting=None,
false_northing=None,
scale_factor_at_projection_origin=None,
ellipsoid=None,
):
"""Construct a RotatedMercator object.
Parameters
----------
latitude_of_projection_origin : float
The true longitude of the central meridian in degrees.
longitude_of_projection_origin : float
The true latitude of the planar origin in degrees.
false_easting : float, optional
X offset from the planar origin in metres.
Defaults to 0.0.
false_northing : float, optional
Y offset from the planar origin in metres.
Defaults to 0.0.
scale_factor_at_projection_origin : float, optional
Scale factor at the central meridian.
Defaults to 1.0 .
ellipsoid : :class:`GeogCS`, optional
If given, defines the ellipsoid.
"""
message = (
"iris.coord_systems.RotatedMercator is deprecated, and will be "
"removed in a future release. Instead please use "
"iris.coord_systems.ObliqueMercator with "
"azimuth_of_central_line=90 ."
)
warn_deprecated(message)
super().__init__(
90.0,
latitude_of_projection_origin,
longitude_of_projection_origin,
false_easting,
false_northing,
scale_factor_at_projection_origin,
ellipsoid,
)
def __repr__(self):
# Remove the azimuth argument from the parent repr.
result = super().__repr__()
result = re.sub(r"azimuth_of_central_line=\d*\.?\d*, ", "", result)
return result