Add Albers Equal Area Projection and some test as per #2496 (#2943)

Add Albers Equal Area Projection and NetCDF rules to handle it

Author: cvelascof

lib/iris/coord_systems.py

@@ -934,3 +934,84 @@ def as_cartopy_crs(self):
 
     def as_cartopy_projection(self):
         return self.as_cartopy_crs()
+
+
+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=0.0,
+                 longitude_of_central_meridian=0.0,
+                 false_easting=0.0, false_northing=0.0,
+                 standard_parallels=(20.0, 50.0),
+                 ellipsoid=None):
+        """
+        Constructs a Albers Conical Equal Area coord system.
+
+        Kwargs:
+
+            * latitude_of_projection_origin
+                    True latitude of planar origin in degrees.
+                    Defaults to 0.
+
+            * longitude_of_central_meridian
+                    True longitude of planar central meridian in degrees.
+                    Defaults to 0.
+
+            * false_easting
+                    X offset from planar origin in metres. Defaults to 0.
+
+            * false_northing
+                    Y offset from planar origin in metres. Defaults to 0.
+
+            * standard_parallels
+                    The one or two latitudes of correct scale.
+                    Defaults to (20,50).
+            * ellipsoid
+                    :class:`GeogCS` defining the ellipsoid.
+
+        """
+        #: True latitude of planar origin in degrees.
+        self.latitude_of_projection_origin = latitude_of_projection_origin
+        #: True longitude of planar central meridian in degrees.
+        self.longitude_of_central_meridian = longitude_of_central_meridian
+        #: X offset from planar origin in metres.
+        self.false_easting = false_easting
+        #: Y offset from planar origin in metres.
+        self.false_northing = false_northing
+        #: The one or two latitudes of correct scale.
+        self.standard_parallels = standard_parallels
+        #: Ellipsoid definition.
+        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)
+
+    def as_cartopy_crs(self):
+        if self.ellipsoid is not None:
+            globe = self.ellipsoid.as_cartopy_globe()
+        else:
+            globe = 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)
+
+    def as_cartopy_projection(self):
+        return self.as_cartopy_crs()

lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb

@@ -176,6 +176,25 @@ fc_provides_grid_mapping_lambert_azimuthal_equal_area
         facts_cf.provides(coordinate_system, lambert_azimuthal_equal_area)
         python engine.rule_triggered.add(rule.name)
 
+#
+# Context:
+#   This rule will trigger iff a grid_mapping() case specific fact
+#   has been asserted that refers to a albers conical equal area.
+#
+# Purpose:
+#   Creates the albers conical equal area coordinate system.
+#
+fc_provides_grid_mapping_albers_equal_area
+    foreach
+        facts_cf.grid_mapping($grid_mapping)
+        check is_grid_mapping(engine, $grid_mapping, CF_GRID_MAPPING_ALBERS)
+    assert
+        python cf_grid_var = engine.cf_var.cf_group.grid_mappings[$grid_mapping]
+        python coordinate_system = build_albers_equal_area_coordinate_system(engine, cf_grid_var)
+        python engine.provides['coordinate_system'] = coordinate_system
+        facts_cf.provides(coordinate_system, albers_equal_area)
+        python engine.rule_triggered.add(rule.name)
+
 
 #
 # Context:
@@ -774,6 +793,45 @@ fc_build_coordinate_projection_y_lambert_azimuthal_equal_area
                                           coord_system=engine.provides['coordinate_system'])
         python engine.rule_triggered.add(rule.name)
 
+#
+# Context:
+#   This rule will trigger iff a projection_x_coordinate coordinate exists and
+#   a albers conical equal area coordinate system exists.
+#
+# Purpose:
+#   Add the projection_x_coordinate coordinate into the cube.
+#
+fc_build_coordinate_projection_x_albers_equal_area
+    foreach
+        facts_cf.provides(coordinate, projection_x_coordinate, $coordinate)
+        facts_cf.provides(coordinate_system, albers_equal_area)
+    assert
+        python cf_coord_var = engine.cf_var.cf_group.coordinates[$coordinate]
+        python build_dimension_coordinate(engine, cf_coord_var,
+                                          coord_name=CF_VALUE_STD_NAME_PROJ_X,
+                                          coord_system=engine.provides['coordinate_system'])
+        python engine.rule_triggered.add(rule.name)
+
+
+#
+# Context:
+#   This rule will trigger iff a projection_y_coordinate coordinate exists and
+#   a albers conical equal area coordinate system exists.
+#
+# Purpose:
+#   Add the projection_y_coordinate coordinate into the cube.
+#
+fc_build_coordinate_projection_y_albers_equal_area
+    foreach
+        facts_cf.provides(coordinate, projection_y_coordinate, $coordinate)
+        facts_cf.provides(coordinate_system, albers_equal_area)
+    assert
+        python cf_coord_var = engine.cf_var.cf_group.coordinates[$coordinate]
+        python build_dimension_coordinate(engine, cf_coord_var,
+                                          coord_name=CF_VALUE_STD_NAME_PROJ_Y,
+                                          coord_system=engine.provides['coordinate_system'])
+        python engine.rule_triggered.add(rule.name)
+
 #
 # Context:
 #   This rule will trigger iff a CF time coordinate exists.
@@ -1387,6 +1445,37 @@ fc_extras
 
         return cs
 
+    ################################################################################
+    def build_albers_equal_area_coordinate_system(engine, cf_grid_var):
+        """
+        Create a albers conical equal area coordinate system from the CF-netCDF
+        grid mapping variable.
+
+        """
+        major, minor, inverse_flattening = _get_ellipsoid(cf_grid_var)
+
+        latitude_of_projection_origin = getattr(
+            cf_grid_var, CF_ATTR_GRID_LAT_OF_PROJ_ORIGIN, None)
+        longitude_of_central_meridian = getattr(
+            cf_grid_var, CF_ATTR_GRID_LON_OF_CENT_MERIDIAN, 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)
+        standard_parallels = getattr(
+            cf_grid_var, CF_ATTR_GRID_STANDARD_PARALLEL, None)
+
+        ellipsoid = None
+        if major is not None or minor is not None or \
+                inverse_flattening is not None:
+            ellipsoid = iris.coord_systems.GeogCS(major, minor,
+                                                  inverse_flattening)
+
+        cs = iris.coord_systems.AlbersEqualArea(
+            latitude_of_projection_origin, longitude_of_central_meridian,
+            false_easting, false_northing, standard_parallels, ellipsoid)
+
+        return cs
 
     ################################################################################
     def get_attr_units(cf_var, attributes):

lib/iris/fileformats/netcdf.py

@@ -1866,6 +1866,19 @@ def add_ellipsoid(ellipsoid):
                     cf_var_grid.false_easting = cs.false_easting
                     cf_var_grid.false_northing = cs.false_northing
 
+                # albers conical equal area
+                elif isinstance(cs,
+                                iris.coord_systems.AlbersEqualArea):
+                    if cs.ellipsoid:
+                        add_ellipsoid(cs.ellipsoid)
+                    cf_var_grid.longitude_of_central_meridian = (
+                        cs.longitude_of_central_meridian)
+                    cf_var_grid.latitude_of_projection_origin = (
+                        cs.latitude_of_projection_origin)
+                    cf_var_grid.false_easting = cs.false_easting
+                    cf_var_grid.false_northing = cs.false_northing
+                    cf_var_grid.standard_parallel = (cs.standard_parallels)
+
                 # other
                 else:
                     warnings.warn('Unable to represent the horizontal '

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

@@ -0,0 +1,100 @@
+# (C) British Crown Copyright 2015 - 2018, Met Office
+#
+# This file is part of Iris.
+#
+# Iris is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Iris is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with Iris.  If not, see <http://www.gnu.org/licenses/>.
+"""
+Unit tests for the :class:`iris.coord_systems.AlbersEqualArea` class.
+
+"""
+
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip)  # noqa
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+import cartopy.crs as ccrs
+from iris.coord_systems import GeogCS, AlbersEqualArea
+
+
+class Test_as_cartopy_crs(tests.IrisTest):
+    def setUp(self):
+        self.latitude_of_projection_origin = 0.0
+        self.longitude_of_central_meridian = 0.0
+        self.semi_major_axis = 6377563.396
+        self.semi_minor_axis = 6356256.909
+        self.false_easting = 0.0
+        self.false_northing = 0.0
+        self.standard_parallels = (-18., -36.)
+        self.ellipsoid = GeogCS(self.semi_major_axis, self.semi_minor_axis)
+        self.aea_cs = AlbersEqualArea(
+            self.latitude_of_projection_origin,
+            self.longitude_of_central_meridian,
+            self.false_easting,
+            self.false_northing,
+            self.standard_parallels,
+            ellipsoid=self.ellipsoid)
+
+    def test_crs_creation(self):
+        res = self.aea_cs.as_cartopy_crs()
+        globe = ccrs.Globe(semimajor_axis=self.semi_major_axis,
+                           semiminor_axis=self.semi_minor_axis,
+                           ellipse=None)
+        expected = ccrs.AlbersEqualArea(
+            self.longitude_of_central_meridian,
+            self.latitude_of_projection_origin,
+            self.false_easting,
+            self.false_northing,
+            self.standard_parallels,
+            globe=globe)
+        self.assertEqual(res, expected)
+
+
+class Test_as_cartopy_projection(tests.IrisTest):
+    def setUp(self):
+        self.latitude_of_projection_origin = 0.0
+        self.longitude_of_central_meridian = 0.0
+        self.semi_major_axis = 6377563.396
+        self.semi_minor_axis = 6356256.909
+        self.false_easting = 0.0
+        self.false_northing = 0.0
+        self.standard_parallels = (-18., -36.)
+        self.ellipsoid = GeogCS(self.semi_major_axis, self.semi_minor_axis)
+        self.aea_cs = AlbersEqualArea(
+            self.latitude_of_projection_origin,
+            self.longitude_of_central_meridian,
+            self.false_easting,
+            self.false_northing,
+            self.standard_parallels,
+            ellipsoid=self.ellipsoid)
+
+    def test_projection_creation(self):
+        res = self.aea_cs.as_cartopy_projection()
+        globe = ccrs.Globe(semimajor_axis=self.semi_major_axis,
+                           semiminor_axis=self.semi_minor_axis,
+                           ellipse=None)
+        expected = ccrs.AlbersEqualArea(
+            self.latitude_of_projection_origin,
+            self.longitude_of_central_meridian,
+            self.false_easting,
+            self.false_northing,
+            self.standard_parallels,
+            globe=globe)
+        self.assertEqual(res, expected)
+
+
+if __name__ == '__main__':
+    tests.main()

lib/iris/tests/unit/fileformats/netcdf/test_Saver.py

@@ -34,7 +34,8 @@
 from iris._lazy_data import as_lazy_data
 from iris.coord_systems import (GeogCS, TransverseMercator, RotatedGeogCS,
                                 LambertConformal, Mercator, Stereographic,
-                                LambertAzimuthalEqualArea)
+                                LambertAzimuthalEqualArea,
+                                AlbersEqualArea)
 from iris.coords import DimCoord
 from iris.cube import Cube
 from iris.fileformats.netcdf import Saver
@@ -785,6 +786,26 @@ def test_laea_cs(self):
                     }
         self._test(coord_system, expected)
 
+    def test_aea_cs(self):
+        coord_system = AlbersEqualArea(
+            latitude_of_projection_origin=52,
+            longitude_of_central_meridian=10,
+            false_easting=100,
+            false_northing=200,
+            standard_parallels=(38, 50),
+            ellipsoid=GeogCS(6377563.396, 6356256.909))
+        expected = {'grid_mapping_name': b'albers_conical_equal_area',
+                    'latitude_of_projection_origin': 52,
+                    'longitude_of_central_meridian': 10,
+                    'false_easting': 100,
+                    'false_northing': 200,
+                    'standard_parallel': (38, 50),
+                    'semi_major_axis': 6377563.396,
+                    'semi_minor_axis': 6356256.909,
+                    'longitude_of_prime_meridian': 0,
+                    }
+        self._test(coord_system, expected)
+
 
 class Test__create_cf_cell_measure_variable(tests.IrisTest):
     # Saving of masked data is disallowed.