Correct handling of NetCDF VerticalPerspective and Geostationary coordinate systems

docs/iris/src/whatsnew/contributions_2.3.0/bugfix_2019-Sep-20_vp_coordsystem_definition.txt

@@ -0,0 +1,4 @@
+* :class:`iris.coord_systems.VerticalPerspective` coordinate system now uses
+the `CF Vertical perspective definition <http://cfconventions
+.org/cf-conventions/cf-conventions.html#vertical-perspective>`_; had been
+erroneously using Geostationary.

docs/iris/src/whatsnew/contributions_2.3.0/newfeature_2019-Sep-20_geostationary_coordsystem.txt

@@ -0,0 +1,5 @@
+* New coordinate system: :class:`iris.coord_systems.Geostationary`, based on
+the `CF Geostationary projection definition <http://cfconventions
+.org/cf-conventions/cf-conventions.html#_geostationary_projection>`_.
+* :class:`iris.coord_systems.VerticalPerspective` can now be saved to and
+loaded from NetCDF files.

lib/iris/coord_systems.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2010 - 2018, Met Office
+# (C) British Crown Copyright 2010 - 2019, Met Office
 #
 # This file is part of Iris.
 #
@@ -76,6 +76,15 @@ def xml_element(self, doc, attrs=None):
 
         return coord_system_xml_element
 
+    @staticmethod
+    def _ellipsoid_to_globe(ellipsoid, globe_default):
+        if ellipsoid is not None:
+            globe = ellipsoid.as_cartopy_globe()
+        else:
+            globe = globe_default
+
+        return globe
+
     @abstractmethod
     def as_cartopy_crs(self):
         """
@@ -338,9 +347,8 @@ def xml_element(self, doc):
         return CoordSystem.xml_element(self, doc, self._pretty_attrs())
 
     def _ccrs_kwargs(self):
-        globe = None
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, None)
+
         # Cartopy v0.12 provided the new arg north_pole_grid_longitude
         cartopy_kwargs = {'pole_longitude': self.grid_north_pole_longitude,
                           'pole_latitude': self.grid_north_pole_latitude,
@@ -441,10 +449,7 @@ def __repr__(self):
                                         self.ellipsoid)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = None
+        globe = self._ellipsoid_to_globe(self.ellipsoid, None)
 
         return ccrs.TransverseMercator(
             central_longitude=self.longitude_of_central_meridian,
@@ -534,10 +539,7 @@ def __repr__(self):
                                         self.ellipsoid)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
 
         warnings.warn('Discarding false_easting and false_northing that are '
                       'not used by Cartopy.')
@@ -553,7 +555,7 @@ def as_cartopy_projection(self):
 
 class VerticalPerspective(CoordSystem):
     """
-    An geostationary satellite image map projection.
+    A vertical/near-side perspective satellite image map projection.
 
     """
 
@@ -563,7 +565,7 @@ def __init__(self, latitude_of_projection_origin,
                  longitude_of_projection_origin, perspective_point_height,
                  false_easting=0, false_northing=0, ellipsoid=None):
         """
-        Constructs an Vertical Perspective Geostationary coord system.
+        Constructs a Vertical Perspective coord system.
 
         Args:
 
@@ -592,34 +594,28 @@ def __init__(self, latitude_of_projection_origin,
         #: True latitude of planar origin in degrees.
         self.latitude_of_projection_origin = float(
             latitude_of_projection_origin)
-        if self.latitude_of_projection_origin != 0.0:
-            raise ValueError('Non-zero latitude of projection currently not'
-                             ' supported by Cartopy.')
 
         #: True longitude of planar origin in degrees.
         self.longitude_of_projection_origin = float(
             longitude_of_projection_origin)
 
         #: Altitude of satellite in metres.
         # test if perspective_point_height may be cast to float for proj.4
-        test_pph = float(perspective_point_height)
-        self.perspective_point_height = perspective_point_height
+        self.perspective_point_height = float(perspective_point_height)
 
         #: X offset from planar origin in metres.
-        test_fe = float(false_easting)
-        self.false_easting = false_easting
+        self.false_easting = float(false_easting)
 
         #: Y offset from planar origin in metres.
-        test_fn = float(false_northing)
-        self.false_northing = false_northing
+        self.false_northing = float(false_northing)
 
         #: Ellipsoid definition.
         self.ellipsoid = ellipsoid
 
     def __repr__(self):
         return "Vertical Perspective(latitude_of_projection_origin={!r}, "\
                "longitude_of_projection_origin={!r}, "\
-               "perspective_point_height = {!r}, "\
+               "perspective_point_height={!r}, "\
                "false_easting={!r}, false_northing={!r}, "\
                "ellipsoid={!r})".format(self.latitude_of_projection_origin,
                                         self.longitude_of_projection_origin,
@@ -629,12 +625,10 @@ def __repr__(self):
                                         self.ellipsoid)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
 
-        return ccrs.Geostationary(
+        return ccrs.NearsidePerspective(
+            central_latitude=self.latitude_of_projection_origin,
             central_longitude=self.longitude_of_projection_origin,
             satellite_height=self.perspective_point_height,
             false_easting=self.false_easting,
@@ -645,6 +639,104 @@ def as_cartopy_projection(self):
         return self.as_cartopy_crs()
 
 
+class Geostationary(CoordSystem):
+    """
+    A geostationary satellite image map projection.
+
+    """
+
+    grid_mapping_name = 'geostationary'
+
+    def __init__(self, latitude_of_projection_origin,
+                 longitude_of_projection_origin,
+                 perspective_point_height, sweep_angle_axis, false_easting=0,
+                 false_northing=0, ellipsoid=None):
+
+        """
+        Constructs a Geostationary coord system.
+
+        Args:
+
+        * latitude_of_projection_origin (float):
+            True latitude of planar origin in degrees.
+
+        * longitude_of_projection_origin (float):
+            True longitude of planar origin in degrees.
+
+        * perspective_point_height (float):
+            Altitude of satellite in metres above the surface of the ellipsoid.
+
+        * sweep_angle_axis (string):
+            The axis along which the satellite instrument sweeps - 'x' or 'y'.
+
+        Kwargs:
+
+        * false_easting (float):
+            X offset from planar origin in metres. Defaults to 0.
+
+        * false_northing (float):
+            Y offset from planar origin in metres. Defaults to 0.
+
+        * ellipsoid (iris.coord_systems.GeogCS):
+            :class:`GeogCS` defining the ellipsoid.
+
+        """
+        #: True latitude of planar origin in degrees.
+        self.latitude_of_projection_origin = float(
+            latitude_of_projection_origin)
+        if self.latitude_of_projection_origin != 0.0:
+            raise ValueError('Non-zero latitude of projection currently not'
+                             ' supported by Cartopy.')
+
+        #: True longitude of planar origin in degrees.
+        self.longitude_of_projection_origin = float(
+            longitude_of_projection_origin)
+
+        #: Altitude of satellite in metres.
+        # test if perspective_point_height may be cast to float for proj.4
+        self.perspective_point_height = float(perspective_point_height)
+
+        #: X offset from planar origin in metres.
+        self.false_easting = float(false_easting)
+
+        #: Y offset from planar origin in metres.
+        self.false_northing = float(false_northing)
+
+        #: The axis along which the satellite instrument sweeps - 'x' or 'y'.
+        self.sweep_angle_axis = sweep_angle_axis
+        if self.sweep_angle_axis not in ('x', 'y'):
+            raise ValueError('Invalid sweep_angle_axis - must be "x" or "y"')
+
+        #: Ellipsoid definition.
+        self.ellipsoid = ellipsoid
+
+    def __repr__(self):
+        return "Geostationary(latitude_of_projection_origin={!r}, " \
+               "longitude_of_projection_origin={!r}, " \
+               "perspective_point_height={!r}, false_easting={!r}, " \
+               "false_northing={!r}, sweep_angle_axis={!r}, " \
+               "ellipsoid={!r}".format(self.latitude_of_projection_origin,
+                                       self.longitude_of_projection_origin,
+                                       self.perspective_point_height,
+                                       self.false_easting,
+                                       self.false_northing,
+                                       self.sweep_angle_axis, self.ellipsoid)
+
+    def as_cartopy_crs(self):
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
+
+        return ccrs.Geostationary(
+            central_longitude=self.longitude_of_projection_origin,
+            satellite_height=self.perspective_point_height,
+            false_easting=self.false_easting,
+            false_northing=self.false_northing,
+            globe=globe,
+            sweep_axis=self.sweep_angle_axis)
+
+    def as_cartopy_projection(self):
+        return self.as_cartopy_crs()
+
+
 class Stereographic(CoordSystem):
     """
     A stereographic map projection.
@@ -712,10 +804,8 @@ def __repr__(self):
                                         self.ellipsoid)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
+
         return ccrs.Stereographic(
             self.central_lat, self.central_lon,
             self.false_easting, self.false_northing,
@@ -807,10 +897,7 @@ def as_cartopy_crs(self):
         else:
             cutoff = None
 
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
 
         # Cartopy v0.12 deprecated the use of secant_latitudes.
         if cartopy.__version__ < '0.12':
@@ -866,10 +953,7 @@ def __repr__(self):
         return res.format(self=self)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
 
         return ccrs.Mercator(
             central_longitude=self.longitude_of_projection_origin,
@@ -935,10 +1019,8 @@ def __repr__(self):
                     self.ellipsoid)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
+
         return ccrs.LambertAzimuthalEqualArea(
             central_longitude=self.longitude_of_projection_origin,
             central_latitude=self.latitude_of_projection_origin,
@@ -1015,10 +1097,8 @@ def __repr__(self):
                     self.ellipsoid)
 
     def as_cartopy_crs(self):
-        if self.ellipsoid is not None:
-            globe = self.ellipsoid.as_cartopy_globe()
-        else:
-            globe = ccrs.Globe()
+        globe = self._ellipsoid_to_globe(self.ellipsoid, ccrs.Globe())
+
         return ccrs.AlbersEqualArea(
             central_longitude=self.longitude_of_central_meridian,
             central_latitude=self.latitude_of_projection_origin,