Implement full Mercator parameter support

.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.7"
+  IRIS_TEST_DATA_VERSION: "2.8"
   # Base directory for the iris-test-data.
   IRIS_TEST_DATA_DIR: ${HOME}/iris-test-data
 

docs/src/whatsnew/latest.rst

@@ -40,6 +40,10 @@ This document explains the changes made to Iris for this release
 #. `@pp-mo`_ fixed cube arithmetic operation for cubes with meshes.
    (:issue:`4454`, :pull:`4651`)
 
+#. `@wjbenfold`_ added support for CF-compliant treatment of
+   ``standard_parallel`` and ``scale_factor_at_projection_origin`` to
+   :class:`~iris.coord_system.Mercator`. (:issue:`3844`, :pull:`4609`)
+
 
 🐛 Bugs Fixed
 =============

lib/iris/coord_systems.py

@@ -1083,6 +1083,7 @@ def __init__(
         longitude_of_projection_origin=None,
         ellipsoid=None,
         standard_parallel=None,
+        scale_factor_at_projection_origin=None,
         false_easting=None,
         false_northing=None,
     ):
@@ -1100,12 +1101,18 @@ def __init__(
         * standard_parallel:
             The latitude where the scale is 1. Defaults to 0.0 .
 
+        * scale_factor_at_projection_origin:
+            Scale factor at natural origin. Defaults to unused.
+
         * false_easting:
             X offset from the planar origin in metres. Defaults to 0.0.
 
         * false_northing:
             Y offset from the planar origin in metres. Defaults to 0.0.
 
+        Note: 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(
@@ -1115,8 +1122,24 @@ def __init__(
         #: 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 = _arg_default(standard_parallel, 0)
+        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)
@@ -1130,6 +1153,8 @@ def __repr__(self):
             "{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})"
         )
@@ -1142,6 +1167,7 @@ def as_cartopy_crs(self):
             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,
         )

lib/iris/fileformats/_nc_load_rules/helpers.py

@@ -445,8 +445,9 @@ def build_mercator_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 Mercator projections with specific
-    # scale_factor_at_projection_origin. This is checked elsewhere.
+    scale_factor_at_projection_origin = getattr(
+        cf_grid_var, CF_ATTR_GRID_SCALE_FACTOR_AT_PROJ_ORIGIN, None
+    )
 
     ellipsoid = None
     if (
@@ -460,6 +461,7 @@ def build_mercator_coordinate_system(engine, cf_grid_var):
         longitude_of_projection_origin,
         ellipsoid=ellipsoid,
         standard_parallel=standard_parallel,
+        scale_factor_at_projection_origin=scale_factor_at_projection_origin,
         false_easting=false_easting,
         false_northing=false_northing,
     )
@@ -1251,17 +1253,20 @@ def has_supported_mercator_parameters(engine, cf_name):
     is_valid = True
     cf_grid_var = engine.cf_var.cf_group[cf_name]
 
+    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 (
         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 "
-            "Mercator projections"
+            "It does not make sense to provide both "
+            '"scale_factor_at_projection_origin" and "standard_parallel".'
         )
         is_valid = False
 

lib/iris/fileformats/netcdf.py

@@ -2563,7 +2563,13 @@ def add_ellipsoid(ellipsoid):
                     )
                     cf_var_grid.false_easting = cs.false_easting
                     cf_var_grid.false_northing = cs.false_northing
-                    cf_var_grid.scale_factor_at_projection_origin = 1.0
+                    # Only one of these should be set
+                    if cs.standard_parallel is not None:
+                        cf_var_grid.standard_parallel = cs.standard_parallel
+                    elif cs.scale_factor_at_projection_origin is not None:
+                        cf_var_grid.scale_factor_at_projection_origin = (
+                            cs.scale_factor_at_projection_origin
+                        )
 
                 # lcc
                 elif isinstance(cs, iris.coord_systems.LambertConformal):

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

@@ -1,2 +1,2 @@
 <?xml version="1.0" ?>
-<mercator ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="90.0" standard_parallel="0.0"/>
+<mercator ellipsoid="GeogCS(semi_major_axis=6377563.396, semi_minor_axis=6356256.909)" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="90.0" scale_factor_at_projection_origin="None" standard_parallel="0.0"/>

lib/iris/tests/results/netcdf/netcdf_merc.cml

@@ -53,15 +53,15 @@
 		 45.5158, 45.9993]]" shape="(192, 192)" standard_name="longitude" units="Unit('degrees')" value_type="float32" var_name="lon"/>
       </coord>
       <coord datadims="[1]">
-        <dimCoord id="970d5aa6" points="[-5.16102e+06, -5.10719e+06, -5.05336e+06, ...,
+        <dimCoord id="732bd7eb" points="[-5.16102e+06, -5.10719e+06, -5.05336e+06, ...,
 		5.01299e+06, 5.06682e+06, 5.12065e+06]" shape="(192,)" standard_name="projection_x_coordinate" units="Unit('m')" value_type="float32" var_name="x">
-          <mercator ellipsoid="GeogCS(6378169.0)" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="0.0" standard_parallel="0.0"/>
+          <mercator ellipsoid="GeogCS(6378169.0)" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="0.0" scale_factor_at_projection_origin="1.0" standard_parallel="None"/>
         </dimCoord>
       </coord>
       <coord datadims="[0]">
-        <dimCoord id="072b8f17" points="[5.16101e+06, 5.10718e+06, 5.05335e+06, ...,
+        <dimCoord id="2a475a6a" points="[5.16101e+06, 5.10718e+06, 5.05335e+06, ...,
 		-5.01294e+06, -5.06678e+06, -5.12061e+06]" shape="(192,)" standard_name="projection_y_coordinate" units="Unit('m')" value_type="float32" var_name="y">
-          <mercator ellipsoid="GeogCS(6378169.0)" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="0.0" standard_parallel="0.0"/>
+          <mercator ellipsoid="GeogCS(6378169.0)" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="0.0" scale_factor_at_projection_origin="1.0" standard_parallel="None"/>
         </dimCoord>
       </coord>
       <coord>

lib/iris/tests/results/netcdf/netcdf_merc_false.cml

@@ -6,17 +6,17 @@
     </attributes>
     <coords>
       <coord datadims="[2]">
-        <dimCoord id="7979c65f" long_name="x-coordinate in Cartesian system" points="[-5950000.0, -5925000.0, -5900000.0, -5875000.0,
+        <dimCoord id="94e59ce8" long_name="x-coordinate in Cartesian system" points="[-5950000.0, -5925000.0, -5900000.0, -5875000.0,
 		-5850000.0, -5825000.0, -5800000.0, -5775000.0,
 		-5750000.0, -5725000.0]" shape="(10,)" standard_name="projection_x_coordinate" units="Unit('m')" value_type="float64" var_name="x">
-          <mercator ellipsoid="GeogCS(6371229.0)" false_easting="-12500.0" false_northing="-12500.0" longitude_of_projection_origin="12.0" standard_parallel="-2.0"/>
+          <mercator ellipsoid="GeogCS(6371229.0)" false_easting="-12500.0" false_northing="-12500.0" longitude_of_projection_origin="12.0" scale_factor_at_projection_origin="None" standard_parallel="-2.0"/>
         </dimCoord>
       </coord>
       <coord datadims="[1]">
-        <dimCoord id="f5c6807e" long_name="y-coordinate in Cartesian system" points="[-6200000.0, -6175000.0, -6150000.0, -6125000.0,
+        <dimCoord id="a72902b5" long_name="y-coordinate in Cartesian system" points="[-6200000.0, -6175000.0, -6150000.0, -6125000.0,
 		-6100000.0, -6075000.0, -6050000.0, -6025000.0,
 		-6000000.0, -5975000.0]" shape="(10,)" standard_name="projection_y_coordinate" units="Unit('m')" value_type="float64" var_name="y">
-          <mercator ellipsoid="GeogCS(6371229.0)" false_easting="-12500.0" false_northing="-12500.0" longitude_of_projection_origin="12.0" standard_parallel="-2.0"/>
+          <mercator ellipsoid="GeogCS(6371229.0)" false_easting="-12500.0" false_northing="-12500.0" longitude_of_projection_origin="12.0" scale_factor_at_projection_origin="None" standard_parallel="-2.0"/>
         </dimCoord>
       </coord>
       <coord datadims="[0]">

lib/iris/tests/results/netcdf/netcdf_merc_scale_factor.cml

@@ -0,0 +1,22 @@
+<?xml version="1.0" ?>
+<cubes xmlns="urn:x-iris:cubeml-0.2">
+  <cube dtype="float64" units="unknown" var_name="wibble">
+    <attributes>
+      <attribute name="Conventions" value="CF-1.7"/>
+    </attributes>
+    <coords>
+      <coord datadims="[1]">
+        <dimCoord id="5ae9804c" points="[0.0, 12.5, 25.0, 37.5, 50.0]" shape="(5,)" standard_name="projection_x_coordinate" units="Unit('m')" value_type="float64" var_name="projection_x_coordinate">
+          <mercator ellipsoid="None" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="0.0" scale_factor_at_projection_origin="1.2" standard_parallel="None"/>
+        </dimCoord>
+      </coord>
+      <coord datadims="[0]">
+        <dimCoord id="502d5aa1" points="[0.0, 12.5, 25.0, 37.5, 50.0]" shape="(5,)" standard_name="projection_y_coordinate" units="Unit('m')" value_type="float64" var_name="projection_y_coordinate">
+          <mercator ellipsoid="None" false_easting="0.0" false_northing="0.0" longitude_of_projection_origin="0.0" scale_factor_at_projection_origin="1.2" standard_parallel="None"/>
+        </dimCoord>
+      </coord>
+    </coords>
+    <cellMethods/>
+    <data checksum="0x1936d0c2" dtype="float64" mask_checksum="no-masked-elements" shape="(5, 5)"/>
+  </cube>
+</cubes>

lib/iris/tests/results/unit/fileformats/netcdf/Saver/write/mercator.cdl

@@ -13,7 +13,7 @@ variables:
 		mercator:longitude_of_projection_origin = 49. ;
 		mercator:false_easting = 0. ;
 		mercator:false_northing = 0. ;
-		mercator:scale_factor_at_projection_origin = 1. ;
+		mercator:standard_parallel = 0. ;
 	int64 projection_y_coordinate(projection_y_coordinate) ;
 		projection_y_coordinate:axis = "Y" ;
 		projection_y_coordinate:units = "m" ;

lib/iris/tests/results/unit/fileformats/netcdf/Saver/write/mercator_no_ellipsoid.cdl

@@ -10,7 +10,7 @@ variables:
 		mercator:longitude_of_projection_origin = 49. ;
 		mercator:false_easting = 0. ;
 		mercator:false_northing = 0. ;
-		mercator:scale_factor_at_projection_origin = 1. ;
+		mercator:standard_parallel = 0. ;
 	int64 projection_y_coordinate(projection_y_coordinate) ;
 		projection_y_coordinate:axis = "Y" ;
 		projection_y_coordinate:units = "m" ;

lib/iris/tests/test_netcdf.py

@@ -218,16 +218,26 @@ def test_load_merc_grid(self):
         )
         self.assertCML(cube, ("netcdf", "netcdf_merc.cml"))
 
-    def test_load_merc_false_en_grid(self):
+    def test_load_complex_merc_grid(self):
         # Test loading a single CF-netCDF file with a Mercator grid_mapping that
-        # includes false easting and northing
+        # includes false easting and northing and a standard parallel
         cube = iris.load_cube(
             tests.get_data_path(
                 ("NetCDF", "mercator", "false_east_north_merc.nc")
             )
         )
         self.assertCML(cube, ("netcdf", "netcdf_merc_false.cml"))
 
+    def test_load_merc_grid_non_unit_scale_factor(self):
+        # Test loading a single CF-netCDF file with a Mercator grid_mapping that
+        # includes a non-unit scale factor at projection origin
+        cube = iris.load_cube(
+            tests.get_data_path(
+                ("NetCDF", "mercator", "non_unit_scale_factor_merc.nc")
+            )
+        )
+        self.assertCML(cube, ("netcdf", "netcdf_merc_scale_factor.cml"))
+
     def test_load_stereographic_grid(self):
         # Test loading a single CF-netCDF file with a stereographic
         # grid_mapping.

lib/iris/tests/unit/coord_systems/test_Mercator.py

@@ -30,6 +30,7 @@ def test_repr(self):
             "ellipsoid=GeogCS(semi_major_axis=6377563.396, "
             "semi_minor_axis=6356256.909), "
             "standard_parallel=0.0, "
+            "scale_factor_at_projection_origin=None, "
             "false_easting=0.0, false_northing=0.0)"
         )
         self.assertEqual(expected, repr(self.tm))
@@ -49,13 +50,21 @@ def test_set_optional_args(self):
         self.assertEqualAndKind(crs.false_easting, 13.0)
         self.assertEqualAndKind(crs.false_northing, 12.0)
 
+    def test_set_optional_scale_factor_alternative(self):
+        # Check that setting the optional (non-ellipse) args works.
+        crs = Mercator(
+            scale_factor_at_projection_origin=1.3,
+        )
+        self.assertEqualAndKind(crs.scale_factor_at_projection_origin, 1.3)
+
     def _check_crs_defaults(self, crs):
         # Check for property defaults when no kwargs options were set.
         # NOTE: except ellipsoid, which is done elsewhere.
         self.assertEqualAndKind(crs.longitude_of_projection_origin, 0.0)
         self.assertEqualAndKind(crs.standard_parallel, 0.0)
         self.assertEqualAndKind(crs.false_easting, 0.0)
         self.assertEqualAndKind(crs.false_northing, 0.0)
+        self.assertEqualAndKind(crs.scale_factor_at_projection_origin, None)
 
     def test_no_optional_args(self):
         # Check expected defaults with no optional args.
@@ -67,6 +76,7 @@ def test_optional_args_None(self):
         crs = Mercator(
             longitude_of_projection_origin=None,
             standard_parallel=None,
+            scale_factor_at_projection_origin=None,
             false_easting=None,
             false_northing=None,
         )
@@ -117,6 +127,31 @@ def test_extra_kwargs(self):
         res = merc_cs.as_cartopy_crs()
         self.assertEqual(res, expected)
 
+    def test_extra_kwargs_scale_factor_alternative(self):
+        # Check that a projection with non-default values is correctly
+        # converted to a cartopy CRS.
+        scale_factor_at_projection_origin = 1.3
+        ellipsoid = GeogCS(
+            semi_major_axis=6377563.396, semi_minor_axis=6356256.909
+        )
+
+        merc_cs = Mercator(
+            ellipsoid=ellipsoid,
+            scale_factor_at_projection_origin=scale_factor_at_projection_origin,
+        )
+
+        expected = ccrs.Mercator(
+            globe=ccrs.Globe(
+                semimajor_axis=6377563.396,
+                semiminor_axis=6356256.909,
+                ellipse=None,
+            ),
+            scale_factor=scale_factor_at_projection_origin,
+        )
+
+        res = merc_cs.as_cartopy_crs()
+        self.assertEqual(res, expected)
+
 
 class Test_as_cartopy_projection(tests.IrisTest):
     def test_simple(self):
@@ -159,6 +194,29 @@ def test_extra_kwargs(self):
         res = merc_cs.as_cartopy_projection()
         self.assertEqual(res, expected)
 
+    def test_extra_kwargs_scale_factor_alternative(self):
+        ellipsoid = GeogCS(
+            semi_major_axis=6377563.396, semi_minor_axis=6356256.909
+        )
+        scale_factor_at_projection_origin = 1.3
+
+        merc_cs = Mercator(
+            ellipsoid=ellipsoid,
+            scale_factor_at_projection_origin=scale_factor_at_projection_origin,
+        )
+
+        expected = ccrs.Mercator(
+            globe=ccrs.Globe(
+                semimajor_axis=6377563.396,
+                semiminor_axis=6356256.909,
+                ellipse=None,
+            ),
+            scale_factor=scale_factor_at_projection_origin,
+        )
+
+        res = merc_cs.as_cartopy_projection()
+        self.assertEqual(res, expected)
+
 
 if __name__ == "__main__":
     tests.main()