Merge pull request #1738 from lbdreyer/lambert-equal

Add lambert azimuthal equal area coord_system

docs/iris/src/whatsnew/contributions_1.11/newfeature_2015-Aug-06_lambert_azimuthal_equal_area.txt

@@ -0,0 +1 @@
+* The coordinate system :class:`iris.coord_systems.LambertAzimuthalEqualArea` has been added with NetCDF saving support.

lib/iris/coord_systems.py

@@ -494,14 +494,14 @@ def __init__(self, latitude_of_projection_origin,
         * longitude_of_projection_origin:
             True longitude of planar origin in degrees.
 
+        Kwargs:
+
         * false_easting
             X offset from planar origin in metres. Defaults to 0.
 
         * false_northing
             Y offset from planar origin in metres. Defaults to 0.
 
-        Kwargs:
-
         * ellipsoid
             :class:`GeogCS` defining the ellipsoid.
 
@@ -577,14 +577,14 @@ def __init__(self, latitude_of_projection_origin,
             Altitude of satellite in metres above the surface of the
             ellipsoid.
 
+        Kwargs:
+
         * false_easting
             X offset from planar origin in metres. Defaults to 0.
 
         * false_northing
             Y offset from planar origin in metres. Defaults to 0.
 
-        Kwargs:
-
         * ellipsoid
             :class:`GeogCS` defining the ellipsoid.
 
@@ -667,14 +667,14 @@ def __init__(self, central_lat, central_lon,
             * central_lon
                     The central longitude, which aligns with the y axis.
 
+        Kwargs:
+
             * false_easting
                     X offset from planar origin in metres. Defaults to 0.
 
             * false_northing
                     Y offset from planar origin in metres. Defaults to 0.
 
-        Kwargs:
-
             * true_scale_lat
                     Latitude of true scale.
 
@@ -739,7 +739,7 @@ def __init__(self, central_lat=39.0, central_lon=-96.0,
         """
         Constructs a LambertConformal coord system.
 
-        Args:
+        Kwargs:
 
             * central_lat
                     The latitude of "unitary scale".
@@ -753,8 +753,6 @@ def __init__(self, central_lat=39.0, central_lon=-96.0,
             * false_northing
                     Y offset from planar origin in metres.
 
-        Kwargs:
-
             * secant_latitudes
                     Latitudes of secant intersection.
 
@@ -837,13 +835,9 @@ def __init__(self, longitude_of_projection_origin=0, ellipsoid=None):
         """
         Constructs a Mercator coord system.
 
-        Args:
-
+        Kwargs:
             * longitude_of_projection_origin
                     True longitude of planar origin in degrees.
-
-        Kwargs:
-
             * ellipsoid
                     :class:`GeogCS` defining the ellipsoid.
 
@@ -870,3 +864,73 @@ def as_cartopy_crs(self):
 
     def as_cartopy_projection(self):
         return self.as_cartopy_crs()
+
+
+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=0.0,
+                 longitude_of_projection_origin=0.0,
+                 false_easting=0.0, false_northing=0.0,
+                 ellipsoid=None):
+        """
+        Constructs a Lambert Azimuthal Equal Area coord system.
+
+        Kwargs:
+
+            * latitude_of_projection_origin
+                    True latitude of planar origin in degrees. Defaults to 0.
+
+            * longitude_of_projection_origin
+                    True longitude of planar origin 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.
+
+            * 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 origin in degrees.
+        self.longitude_of_projection_origin = longitude_of_projection_origin
+        #: X offset from planar origin in metres.
+        self.false_easting = false_easting
+        #: Y offset from planar origin in metres.
+        self.false_northing = false_northing
+        #: Ellipsoid definition.
+        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)
+
+    def as_cartopy_crs(self):
+        if self.ellipsoid is not None:
+            globe = self.ellipsoid.as_cartopy_globe()
+        else:
+            globe = 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)
+
+    def as_cartopy_projection(self):
+        return self.as_cartopy_crs()

lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb

@@ -157,6 +157,26 @@ fc_provides_grid_mapping_lambert_conformal
         facts_cf.provides(coordinate_system, lambert_conformal)
         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 lambert azimuthal equal area.
+#
+# Purpose:
+#   Creates the lambert azimuthal equal area coordinate system.
+#
+fc_provides_grid_mapping_lambert_azimuthal_equal_area
+    foreach
+        facts_cf.grid_mapping($grid_mapping)
+        check is_grid_mapping(engine, $grid_mapping, CF_GRID_MAPPING_LAMBERT_AZIMUTHAL)
+    assert
+        python cf_grid_var = engine.cf_var.cf_group.grid_mappings[$grid_mapping]
+        python coordinate_system = build_lambert_azimuthal_equal_area_coordinate_system(engine, cf_grid_var)
+        python engine.provides['coordinate_system'] = coordinate_system
+        facts_cf.provides(coordinate_system, lambert_azimuthal_equal_area)
+        python engine.rule_triggered.add(rule.name)
+
+
 #
 # Context:
 #   This rule will trigger iff a coordinate() case specific fact
@@ -715,6 +735,45 @@ fc_build_coordinate_projection_y_stereographic
         python engine.rule_triggered.add(rule.name)
 
 
+#
+# Context:
+#   This rule will trigger iff a projection_x_coordinate coordinate exists and
+#   a lambert azimuthal equal area coordinate system exists.
+#
+# Purpose:
+#   Add the projection_x_coordinate coordinate into the cube.
+#
+fc_build_coordinate_projection_x_lambert_azimuthal_equal_area
+    foreach
+        facts_cf.provides(coordinate, projection_x_coordinate, $coordinate)
+        facts_cf.provides(coordinate_system, lambert_azimuthal_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 lambert azimuthal equal area coordinate system exists.
+#
+# Purpose:
+#   Add the projection_y_coordinate coordinate into the cube.
+#
+fc_build_coordinate_projection_y_lambert_azimuthal_equal_area
+    foreach
+        facts_cf.provides(coordinate, projection_y_coordinate, $coordinate)
+        facts_cf.provides(coordinate_system, lambert_azimuthal_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.
@@ -1298,6 +1357,37 @@ fc_extras
         return cs
 
 
+    ################################################################################
+    def build_lambert_azimuthal_equal_area_coordinate_system(engine, cf_grid_var):
+        """
+        Create a lambert azimuthal 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_projection_origin = getattr(
+            cf_grid_var, CF_ATTR_GRID_LON_OF_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)
+
+        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.LambertAzimuthalEqualArea(
+            latitude_of_projection_origin, longitude_of_projection_origin,
+            false_easting, false_northing, ellipsoid)
+
+        return cs
+
+
     ################################################################################
     def get_attr_units(cf_var, attributes):
         attr_units = getattr(cf_var, CF_ATTR_UNITS, cf_units._UNIT_DIMENSIONLESS)

lib/iris/fileformats/netcdf.py

@@ -1300,7 +1300,8 @@ def _get_dim_names(self, cube):
                 dimension_names.append(dim_name)
         return dimension_names
 
-    def _cf_coord_identity(self, coord):
+    @staticmethod
+    def _cf_coord_identity(coord):
         """
         Determine a suitable units from a given coordinate.
 
@@ -1327,7 +1328,7 @@ def _cf_coord_identity(self, coord):
             elif coord.standard_name == "longitude":
                 units = 'degrees_east'
 
-        return (coord.standard_name, coord.long_name, units)
+        return coord.standard_name, coord.long_name, units
 
     def _ensure_valid_dtype(self, values, src_name, src_object):
         # NetCDF3 does not support int64 or unsigned ints, so we check
@@ -1779,6 +1780,18 @@ def add_ellipsoid(ellipsoid):
                 elif isinstance(cs, iris.coord_systems.OSGB):
                     warnings.warn('OSGB coordinate system not yet handled')
 
+                # lambert azimuthal equal area
+                elif isinstance(cs,
+                                iris.coord_systems.LambertAzimuthalEqualArea):
+                    if cs.ellipsoid:
+                        add_ellipsoid(cs.ellipsoid)
+                    cf_var_grid.longitude_of_projection_origin = (
+                        cs.longitude_of_projection_origin)
+                    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
+
                 # other
                 else:
                     warnings.warn('Unable to represent the horizontal '

lib/iris/tests/integration/test_netcdf.py

@@ -301,5 +301,13 @@ def test_unknown_method(self):
             shutil.rmtree(temp_dirpath)
 
 
+class TestCoordSystem(tests.IrisTest):
+    def test_load_laea_grid(self):
+        cube = iris.load_cube(
+            tests.get_data_path(('NetCDF', 'lambert_azimuthal_equal_area',
+                                 'euro_air_temp.nc')))
+        self.assertCML(cube, ('netcdf', 'netcdf_laea.cml'))
+
+
 if __name__ == "__main__":
     tests.main()

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

@@ -0,0 +1,72 @@
+<?xml version="1.0" ?>
+<cubes xmlns="urn:x-iris:cubeml-0.2">
+  <cube standard_name="air_temperature" units="K" var_name="air_temperature">
+    <attributes>
+      <attribute name="Conventions" value="CF-1.5"/>
+      <attribute name="STASH" value="m01s16i203"/>
+      <attribute name="source" value="Data from Met Office Unified Model"/>
+    </attributes>
+    <coords>
+      <coord>
+        <dimCoord id="1d45e087" points="[6477.0]" shape="(1,)" standard_name="forecast_period" units="Unit('hours')" value_type="float64" var_name="forecast_period"/>
+      </coord>
+      <coord>
+        <dimCoord id="9c8bdf81" points="[246987.0]" shape="(1,)" standard_name="forecast_reference_time" units="Unit('hours since 1970-01-01 00:00:00', calendar='gregorian')" value_type="float64" var_name="forecast_reference_time"/>
+      </coord>
+      <coord>
+        <dimCoord id="6eef7051" long_name="pressure" points="[1000.0]" shape="(1,)" units="Unit('hPa')" value_type="float32" var_name="pressure"/>
+      </coord>
+      <coord datadims="[1]">
+        <dimCoord bounds="[[430357.142857, 869642.857143],
+		[869642.857143, 1308928.57143],
+		[1308928.57143, 1748214.28571],
+		[1748214.28571, 2187500.0],
+		[2187500.0, 2626785.71429],
+		[2626785.71429, 3066071.42857],
+		[3066071.42857, 3505357.14286],
+		[3505357.14286, 3944642.85714],
+		[3944642.85714, 4383928.57143],
+		[4383928.57143, 4823214.28571],
+		[4823214.28571, 5262500.0],
+		[5262500.0, 5701785.71429],
+		[5701785.71429, 6141071.42857],
+		[6141071.42857, 6580357.14286],
+		[6580357.14286, 7019642.85714]]" id="950c6ce8" points="[650000.0, 1089285.71429, 1528571.42857,
+		1967857.14286, 2407142.85714, 2846428.57143,
+		3285714.28571, 3725000.0, 4164285.71429,
+		4603571.42857, 5042857.14286, 5482142.85714,
+		5921428.57143, 6360714.28571, 6800000.0]" shape="(15,)" standard_name="projection_x_coordinate" units="Unit('m')" value_type="float64" var_name="projection_x_coordinate">
+          <lambertAzimuthalEqualArea ellipsoid="None" false_easting="4321000" false_northing="3210000" latitude_of_projection_origin="52" longitude_of_projection_origin="10"/>
+        </dimCoord>
+      </coord>
+      <coord datadims="[0]">
+        <dimCoord bounds="[[414285.714286, 785714.285714],
+		[785714.285714, 1157142.85714],
+		[1157142.85714, 1528571.42857],
+		[1528571.42857, 1900000.0],
+		[1900000.0, 2271428.57143],
+		[2271428.57143, 2642857.14286],
+		[2642857.14286, 3014285.71429],
+		[3014285.71429, 3385714.28571],
+		[3385714.28571, 3757142.85714],
+		[3757142.85714, 4128571.42857],
+		[4128571.42857, 4500000.0],
+		[4500000.0, 4871428.57143],
+		[4871428.57143, 5242857.14286],
+		[5242857.14286, 5614285.71429],
+		[5614285.71429, 5985714.28571]]" id="fbb8fa7a" points="[600000.0, 971428.571429, 1342857.14286,
+		1714285.71429, 2085714.28571, 2457142.85714,
+		2828571.42857, 3200000.0, 3571428.57143,
+		3942857.14286, 4314285.71429, 4685714.28571,
+		5057142.85714, 5428571.42857, 5800000.0]" shape="(15,)" standard_name="projection_y_coordinate" units="Unit('m')" value_type="float64" var_name="projection_y_coordinate">
+          <lambertAzimuthalEqualArea ellipsoid="None" false_easting="4321000" false_northing="3210000" latitude_of_projection_origin="52" longitude_of_projection_origin="10"/>
+        </dimCoord>
+      </coord>
+      <coord>
+        <dimCoord id="cb784457" points="[253464.0]" shape="(1,)" standard_name="time" units="Unit('hours since 1970-01-01 00:00:00', calendar='gregorian')" value_type="float64" var_name="time"/>
+      </coord>
+    </coords>
+    <cellMethods/>
+    <data checksum="0xf8a2bc63" dtype="float32" shape="(15, 15)"/>
+  </cube>
+</cubes>