Fix up the satellite projections (Fixes #1144)

lib/cartopy/crs.py

@@ -2003,20 +2003,14 @@ def __init__(self, projection, satellite_height=35785831,
             proj4_params.append(('sweep', sweep_axis))
         super(_Satellite, self).__init__(proj4_params, globe=globe)
 
-        # TODO: Let the globe return the semimajor axis always.
-        a = np.float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS)
-        b = np.float(self.globe.semiminor_axis or a)
-        h = np.float(satellite_height)
-        max_x = h * np.arcsin(a / (a + h))
-        max_y = h * np.arcsin(b / (a + h))
-
+    def _set_bounds(self, max_x, max_y):
+        false_easting = self.proj4_params['x_0']
+        false_northing = self.proj4_params['y_0']
         coords = _ellipse_boundary(max_x, max_y,
                                    false_easting, false_northing, 61)
         self._boundary = sgeom.LinearRing(coords.T)
-        mins = np.min(coords, axis=1)
-        maxs = np.max(coords, axis=1)
-        self._x_limits = mins[0], maxs[0]
-        self._y_limits = mins[1], maxs[1]
+        self._x_limits = -max_x + false_easting, max_x + false_easting
+        self._y_limits = -max_y + false_northing, max_y + false_northing
         self._threshold = np.diff(self._x_limits)[0] * 0.02
 
     @property
@@ -2038,12 +2032,38 @@ def y_limits(self):
 
 class Geostationary(_Satellite):
     """
+    A view appropriate for satellites in Geostationary Earth orbit.
+
     Perspective view looking directly down from above a point on the equator.
 
+    In this projection, the projected coordinates are scanning angles measured
+    from the satellite looking directly downward, multiplied by the height of
+    the satellite.
+
     """
     def __init__(self, central_longitude=0.0, satellite_height=35785831,
                  false_easting=0, false_northing=0, globe=None,
                  sweep_axis='y'):
+        """
+        Parameters
+        ----------
+        central_longitude: float, optional
+            The central longitude. Defaults to 0.
+        satellite_height: float, optional
+            The height of the satellite. Defaults to 35785831 meters
+            (true geostationary orbit).
+        false_easting:
+            X offset from planar origin in metres. Defaults to 0.
+        false_northing:
+            Y offset from planar origin in metres. Defaults to 0.
+        globe: :class:`cartopy.crs.Globe`, optional
+            If omitted, a default globe is created.
+        sweep_axis: 'x' or 'y', optional. Defaults to 'y'.
+            Controls which axis is scanned first, and thus which angle is
+            applied first. The default is appropriate for Meteosat, while
+            'x' should be used for GOES.
+        """
+
         super(Geostationary, self).__init__(
             projection='geos',
             satellite_height=satellite_height,
@@ -2054,15 +2074,59 @@ def __init__(self, central_longitude=0.0, satellite_height=35785831,
             globe=globe,
             sweep_axis=sweep_axis)
 
+        # TODO: Let the globe return the semimajor axis always.
+        a = np.float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS)
+        b = np.float(self.globe.semiminor_axis or a)
+        h = np.float(satellite_height)
+        max_x = h * np.arcsin(a / (a + h))
+        max_y = h * np.arcsin(b / (a + h))
+        self._set_bounds(max_x, max_y)
+
 
 class NearsidePerspective(_Satellite):
     """
     Perspective view looking directly down from above a point on the globe.
 
+    In this projection, the projected coordinates are x and y measured from
+    the origin of a plane tangent to the Earth directly below the perspective
+    point (e.g. a satellite).
+
     """
     def __init__(self, central_longitude=0.0, central_latitude=0.0,
                  satellite_height=35785831,
                  false_easting=0, false_northing=0, globe=None):
+        """
+        Parameters
+        ----------
+        central_longitude: float, optional
+            The central longitude. Defaults to 0.
+        central_latitude: float, optional
+            The central latitude. Defaults to 0.
+        satellite_height: float, optional
+            The height of the satellite. Defaults to 35785831 meters
+            (true geostationary orbit).
+        false_easting:
+            X offset from planar origin in metres. Defaults to 0.
+        false_northing:
+            Y offset from planar origin in metres. Defaults to 0.
+        globe: :class:`cartopy.crs.Globe`, optional
+            If omitted, a default globe is created.
+
+            .. note::
+                This projection does not handle elliptical globes.
+
+        """
+        if globe is None:
+            globe = Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS, ellipse=None)
+
+        # TODO: Let the globe return the semimajor axis always.
+        a = globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS
+        b = globe.semiminor_axis or a
+
+        if b != a or globe.ellipse is not None:
+            warnings.warn('The proj "nsper" projection does not handle '
+                          'elliptical globes.')
+
         super(NearsidePerspective, self).__init__(
             projection='nsper',
             satellite_height=satellite_height,
@@ -2072,6 +2136,10 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0,
             false_northing=false_northing,
             globe=globe)
 
+        h = np.float(satellite_height)
+        max_x = a * np.sqrt(h / (2 * a + h))
+        self._set_bounds(max_x, max_x)
+
 
 class AlbersEqualArea(Projection):
     """

lib/cartopy/tests/crs/test_geostationary.py

@@ -26,14 +26,13 @@
 import cartopy.crs as ccrs
 
 
-# Note: code here is now shared with the NearsidePerspective test.
 def check_proj4_params(name, crs, other_args):
     expected = other_args | {'proj={}'.format(name), 'units=m', 'no_defs'}
     pro4_params = set(crs.proj4_init.lstrip('+').split(' +'))
     assert expected == pro4_params
 
 
-class GeostationaryTestsMixin(object):
+class TestGeostationary(object):
     test_class = ccrs.Geostationary
     expected_proj_name = 'geos'
 
@@ -84,8 +83,6 @@ def test_eastings(self):
                              10434177.81588539, 5309177.81588539),
                             decimal=4)
 
-
-class TestGeostationary(GeostationaryTestsMixin, object):
     def test_sweep(self):
         geos = ccrs.Geostationary(sweep_axis='x')
         other_args = {'ellps=WGS84', 'h=35785831', 'lat_0=0.0', 'lon_0=0.0',

lib/cartopy/tests/crs/test_nearside_perspective.py

@@ -23,28 +23,46 @@
 
 from numpy.testing import assert_almost_equal
 
-from cartopy.tests.crs.test_geostationary import (GeostationaryTestsMixin,
-                                                  check_proj4_params)
+from cartopy.tests.crs.test_geostationary import check_proj4_params
 
-from cartopy.crs import NearsidePerspective
+import cartopy.crs as ccrs
 
 
-class TestEquatorialDefault(GeostationaryTestsMixin, object):
-    # Check that it behaves just like Geostationary, in the absence of a
-    # central_latitude parameter.
-    test_class = NearsidePerspective
-    expected_proj_name = 'nsper'
+def test_default():
+    geos = ccrs.NearsidePerspective()
+    other_args = {'a=6378137.0', 'h=35785831', 'lat_0=0.0', 'lon_0=0.0',
+                  'x_0=0', 'y_0=0'}
 
+    check_proj4_params('nsper', geos, other_args)
 
-class TestOwnSpecifics(object):
-    def test_central_latitude(self):
-        # Check the effect of the added 'central_latitude' key.
-        geos = NearsidePerspective(central_latitude=53.7)
-        other_args = {'ellps=WGS84', 'h=35785831', 'lat_0=53.7', 'lon_0=0.0',
-                      'x_0=0', 'y_0=0'}
-        check_proj4_params('nsper', geos, other_args)
+    assert_almost_equal(geos.boundary.bounds,
+                        (-5476336.098, -5476336.098,
+                         5476336.098, 5476336.098),
+                        decimal=3)
 
-        assert_almost_equal(geos.boundary.bounds,
-                            (-5434177.81588539, -5434177.81588539,
-                             5434177.81588539, 5434177.81588539),
-                            decimal=4)
+
+def test_offset():
+    geos = ccrs.NearsidePerspective(false_easting=5000000,
+                                    false_northing=-123000,)
+    other_args = {'a=6378137.0', 'h=35785831', 'lat_0=0.0', 'lon_0=0.0',
+                  'x_0=5000000', 'y_0=-123000'}
+
+    check_proj4_params('nsper', geos, other_args)
+
+    assert_almost_equal(geos.boundary.bounds,
+                        (-476336.098, -5599336.098,
+                         10476336.098, 5353336.098),
+                        decimal=4)
+
+
+def test_central_latitude():
+    # Check the effect of the added 'central_latitude' key.
+    geos = ccrs.NearsidePerspective(central_latitude=53.7)
+    other_args = {'a=6378137.0', 'h=35785831', 'lat_0=53.7', 'lon_0=0.0',
+                  'x_0=0', 'y_0=0'}
+    check_proj4_params('nsper', geos, other_args)
+
+    assert_almost_equal(geos.boundary.bounds,
+                        (-5476336.098, -5476336.098,
+                         5476336.098, 5476336.098),
+                        decimal=3)