diff --git a/lib/iris/coord_systems.py b/lib/iris/coord_systems.py index 651067dd46..0bc770f9b8 100644 --- a/lib/iris/coord_systems.py +++ b/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() diff --git a/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb b/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb index 6274bef70b..942450847b 100644 --- a/lib/iris/fileformats/_pyke_rules/fc_rules_cf.krb +++ b/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): diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py index 5aa3e54a11..6ac9c26811 100644 --- a/lib/iris/fileformats/netcdf.py +++ b/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 ' diff --git a/lib/iris/tests/unit/coord_systems/test_AlbersEqualArea.py b/lib/iris/tests/unit/coord_systems/test_AlbersEqualArea.py new file mode 100644 index 0000000000..dd428fe1c8 --- /dev/null +++ b/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 . +""" +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() diff --git a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py b/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py index f4b77df15c..76a5c13861 100644 --- a/lib/iris/tests/unit/fileformats/netcdf/test_Saver.py +++ b/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.