PolarStereographic implementation

.cirrus.yml

@@ -38,7 +38,7 @@ env:
   # Conda packages to be installed.
   CONDA_CACHE_PACKAGES: "nox pip"
   # Git commit hash for iris test data.
-  IRIS_TEST_DATA_VERSION: "2.8"
+  IRIS_TEST_DATA_VERSION: "2.9"
   # Base directory for the iris-test-data.
   IRIS_TEST_DATA_DIR: ${HOME}/iris-test-data
 

docs/src/whatsnew/latest.rst

@@ -64,6 +64,9 @@ This document explains the changes made to Iris for this release
    preserving the time of year. (:issue:`1422`, :issue:`4098`, :issue:`4665`,
    :pull:`4723`)
 
+#. `@wjbenfold`_ implemented the :class:`~iris.coord_systems.PolarStereographic`
+   CRS. (:issue:`4770`, :pull:`4773`)
+
 
 🐛 Bugs Fixed
 =============

lib/iris/coord_systems.py

@@ -1061,32 +1061,39 @@ def __init__(
         false_northing=None,
         true_scale_lat=None,
         ellipsoid=None,
+        scale_factor_at_projection_origin=None,
     ):
         """
         Constructs a Stereographic coord system.
 
-        Args:
+        Parameters
+        ----------
 
-        * central_lat:
+        central_lat : float
             The latitude of the pole.
 
-        * central_lon:
+        central_lon : float
             The central longitude, which aligns with the y axis.
 
-        Kwargs:
-
-        * false_easting:
-            X offset from planar origin in metres. Defaults to 0.0 .
+        false_easting : float, optional
+            X offset from planar origin in metres.
 
-        * false_northing:
-            Y offset from planar origin in metres. Defaults to 0.0 .
+        false_northing : float, optional
+            Y offset from planar origin in metres.
 
-        * true_scale_lat:
+        true_scale_lat : float, optional
             Latitude of true scale.
 
-        * ellipsoid (:class:`GeogCS`):
+        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 one of true_scale_lat and scale_factor_at_projection_origin
+
         """
 
         #: True latitude of planar origin in degrees.
@@ -1105,27 +1112,42 @@ def __init__(
         self.true_scale_lat = _arg_default(
             true_scale_lat, None, cast_as=_float_or_None
         )
-        # N.B. the way we use this parameter, we need it to default to 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__(self):
-        return (
-            "Stereographic(central_lat={!r}, central_lon={!r}, "
-            "false_easting={!r}, false_northing={!r}, "
-            "true_scale_lat={!r}, "
-            "ellipsoid={!r})".format(
-                self.central_lat,
-                self.central_lon,
-                self.false_easting,
-                self.false_northing,
-                self.true_scale_lat,
-                self.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()
+
     def as_cartopy_crs(self):
         globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
 
@@ -1135,13 +1157,81 @@ def as_cartopy_crs(self):
             self.false_easting,
             self.false_northing,
             self.true_scale_lat,
+            self.scale_factor_at_projection_origin,
             globe=globe,
         )
 
     def as_cartopy_projection(self):
         return self.as_cartopy_crs()
 
 
+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 one of `true_scale_lat` and `scale_factor_at_projection_origin`.
+        TODO: What if you give neither?
+
+
+        """
+
+        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()
+
+
 class LambertConformal(CoordSystem):
     """
     A coordinate system in the Lambert Conformal conic projection.

lib/iris/fileformats/_nc_load_rules/actions.py

@@ -110,9 +110,13 @@ def action_default(engine):
         hh.build_transverse_mercator_coordinate_system,
     ),
     hh.CF_GRID_MAPPING_STEREO: (
-        hh.has_supported_stereographic_parameters,
+        None,
         hh.build_stereographic_coordinate_system,
     ),
+    hh.CF_GRID_MAPPING_POLAR: (
+        hh.has_supported_polar_stereographic_parameters,
+        hh.build_polar_stereographic_coordinate_system,
+    ),
     hh.CF_GRID_MAPPING_LAMBERT_CONFORMAL: (
         None,
         hh.build_lambert_conformal_coordinate_system,

lib/iris/fileformats/_nc_load_rules/helpers.py

@@ -145,6 +145,7 @@
 CF_ATTR_GRID_SEMI_MINOR_AXIS = "semi_minor_axis"
 CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN = "latitude_of_projection_origin"
 CF_ATTR_GRID_LON_OF_PROJ_ORIGIN = "longitude_of_projection_origin"
+CF_ATTR_GRID_STRAIGHT_VERT_LON = "straight_vertical_longitude_from_pole"
 CF_ATTR_GRID_STANDARD_PARALLEL = "standard_parallel"
 CF_ATTR_GRID_FALSE_EASTING = "false_easting"
 CF_ATTR_GRID_FALSE_NORTHING = "false_northing"
@@ -418,8 +419,6 @@ def build_stereographic_coordinate_system(engine, cf_grid_var):
     )
     false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None)
     false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None)
-    # Iris currently only supports Stereographic projections with a scale
-    # factor of 1.0. This is checked elsewhere.
 
     cs = iris.coord_systems.Stereographic(
         latitude_of_projection_origin,
@@ -433,6 +432,42 @@ def build_stereographic_coordinate_system(engine, cf_grid_var):
     return cs
 
 
+################################################################################
+def build_polar_stereographic_coordinate_system(engine, cf_grid_var):
+    """
+    Create a polar stereographic coordinate system from the CF-netCDF
+    grid mapping variable.
+
+    """
+    ellipsoid = _get_ellipsoid(cf_grid_var)
+
+    latitude_of_projection_origin = getattr(
+        cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None
+    )
+    longitude_of_projection_origin = getattr(
+        cf_grid_var, CF_ATTR_GRID_STRAIGHT_VERT_LON, None
+    )
+    true_scale_lat = getattr(cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None)
+    scale_factor_at_projection_origin = getattr(
+        cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None
+    )
+
+    false_easting = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_EASTING, None)
+    false_northing = getattr(cf_grid_var, CF_ATTR_GRID_FALSE_NORTHING, None)
+
+    cs = iris.coord_systems.PolarStereographic(
+        latitude_of_projection_origin,
+        longitude_of_projection_origin,
+        false_easting,
+        false_northing,
+        true_scale_lat,
+        scale_factor_at_projection_origin,
+        ellipsoid=ellipsoid,
+    )
+
+    return cs
+
+
 ################################################################################
 def build_mercator_coordinate_system(engine, cf_grid_var):
     """
@@ -1239,24 +1274,45 @@ def has_supported_mercator_parameters(engine, cf_name):
 
 
 ################################################################################
-def has_supported_stereographic_parameters(engine, cf_name):
-    """Determine whether the CF grid mapping variable has a value of 1.0
-    for the scale_factor_at_projection_origin attribute."""
+def has_supported_polar_stereographic_parameters(engine, cf_name):
+    """Determine whether the CF grid mapping variable has the supported
+    values for the parameters of the Polar Stereographic projection."""
 
     is_valid = True
     cf_grid_var = engine.cf_var.cf_group[cf_name]
 
+    latitude_of_projection_origin = getattr(
+        cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None
+    )
+
+    standard_parallel = getattr(
+        cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None
+    )
     scale_factor_at_projection_origin = getattr(
         cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None
     )
 
+    if (
+        latitude_of_projection_origin != 90
+        and latitude_of_projection_origin != -90
+    ):
+        warnings.warn('"latitude_of_projection_origin" must be +90 or -90.')
+        is_valid = False
+
     if (
         scale_factor_at_projection_origin is not None
-        and scale_factor_at_projection_origin != 1
+        and standard_parallel is not None
     ):
         warnings.warn(
-            "Scale factors other than 1.0 not yet supported for "
-            "stereographic projections"
+            "It does not make sense to provide both "
+            '"scale_factor_at_projection_origin" and "standard_parallel".'
+        )
+        is_valid = False
+
+    if scale_factor_at_projection_origin is None and standard_parallel is None:
+        warnings.warn(
+            'One of "scale_factor_at_projection_origin" and '
+            '"standard_parallel" is required.'
         )
         is_valid = False
 

lib/iris/fileformats/netcdf.py

@@ -2685,6 +2685,28 @@ def add_ellipsoid(ellipsoid):
                     cf_var_grid.false_easting = cs.false_easting
                     cf_var_grid.false_northing = cs.false_northing
 
+                # polar stereo (have to do this before Stereographic because it subclasses it)
+                elif isinstance(cs, iris.coord_systems.PolarStereographic):
+                    if cs.ellipsoid:
+                        add_ellipsoid(cs.ellipsoid)
+                    cf_var_grid.latitude_of_projection_origin = cs.central_lat
+                    cf_var_grid.straight_vertical_longitude_from_pole = (
+                        cs.central_lon
+                    )
+                    cf_var_grid.false_easting = cs.false_easting
+                    cf_var_grid.false_northing = cs.false_northing
+                    # Only one of these should be set
+                    if cs.true_scale_lat is not None:
+                        cf_var_grid.true_scale_lat = cs.true_scale_lat
+                    elif cs.scale_factor_at_projection_origin is not None:
+                        cf_var_grid.scale_factor_at_projection_origin = (
+                            cs.scale_factor_at_projection_origin
+                        )
+                    else:
+                        cf_var_grid.scale_factor_at_projection_origin = (
+                            1.0  # TODO: Is this right?
+                        )
+
                 # stereo
                 elif isinstance(cs, iris.coord_systems.Stereographic):
                     if cs.true_scale_lat is not None:

lib/iris/tests/results/coord_systems/PolarStereographic.xml

@@ -0,0 +1,2 @@
+<?xml version="1.0" ?>
+<polarStereographic central_lat="90.0" central_lon="0.0" ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="0.0" false_northing="0.0" scale_factor_at_projection_origin="None" true_scale_lat="None"/>

lib/iris/tests/results/coord_systems/PolarStereographicScaleFactor.xml

@@ -0,0 +1,2 @@
+<?xml version="1.0" ?>
+<polarStereographic central_lat="90.0" central_lon="0.0" ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="0.0" false_northing="0.0" scale_factor_at_projection_origin="1.1" true_scale_lat="None"/>

lib/iris/tests/results/coord_systems/PolarStereographicStandardParallel.xml

@@ -0,0 +1,2 @@
+<?xml version="1.0" ?>
+<polarStereographic central_lat="90.0" central_lon="0.0" ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="0.0" false_northing="0.0" scale_factor_at_projection_origin="None" true_scale_lat="30.0"/>

lib/iris/tests/results/coord_systems/Stereographic.xml

@@ -1,2 +1,2 @@
 <?xml version="1.0" ?>
-<stereographic central_lat="-90.0" central_lon="-45.0" ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="100.0" false_northing="200.0" true_scale_lat="None"/>
+<stereographic central_lat="-90.0" central_lon="-45.0" ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="100.0" false_northing="200.0" scale_factor_at_projection_origin="None" true_scale_lat="None"/>