Merge pull request #929 from dopplershift/fix_issues_339_455

Fix stereographic scaling issues

Author: QuLogic

lib/cartopy/crs.py

@@ -1310,12 +1310,27 @@ def y_limits(self):
 class Stereographic(Projection):
     def __init__(self, central_latitude=0.0, central_longitude=0.0,
                  false_easting=0.0, false_northing=0.0,
-                 true_scale_latitude=None, globe=None):
+                 true_scale_latitude=None,
+                 scale_factor=None, globe=None):
         proj4_params = [('proj', 'stere'), ('lat_0', central_latitude),
                         ('lon_0', central_longitude),
                         ('x_0', false_easting), ('y_0', false_northing)]
-        if true_scale_latitude:
+
+        if true_scale_latitude is not None:
+            if central_latitude not in (-90., 90.):
+                warnings.warn('"true_scale_latitude" parameter is only used '
+                              'for polar stereographic projections. Consider '
+                              'the use of "scale_factor" instead.')
             proj4_params.append(('lat_ts', true_scale_latitude))
+
+        if scale_factor is not None:
+            if true_scale_latitude is not None:
+                raise ValueError('It does not make sense to provide both '
+                                 '"scale_factor" and "true_scale_latitude". '
+                                 'Ignoring "scale_factor".')
+            else:
+                proj4_params.append(('k_0', scale_factor))
+
         super(Stereographic, self).__init__(proj4_params, globe=globe)
 
         # TODO: Let the globe return the semimajor axis always.
@@ -1358,17 +1373,23 @@ def y_limits(self):
 
 
 class NorthPolarStereo(Stereographic):
-    def __init__(self, central_longitude=0.0, globe=None):
+    def __init__(self, central_longitude=0.0, true_scale_latitude=None,
+                 globe=None):
         super(NorthPolarStereo, self).__init__(
             central_latitude=90,
-            central_longitude=central_longitude, globe=globe)
+            central_longitude=central_longitude,
+            true_scale_latitude=true_scale_latitude,  # None is +90
+            globe=globe)
 
 
 class SouthPolarStereo(Stereographic):
-    def __init__(self, central_longitude=0.0, globe=None):
+    def __init__(self, central_longitude=0.0, true_scale_latitude=None,
+                 globe=None):
         super(SouthPolarStereo, self).__init__(
             central_latitude=-90,
-            central_longitude=central_longitude, globe=globe)
+            central_longitude=central_longitude,
+            true_scale_latitude=true_scale_latitude,  # None is -90
+            globe=globe)
 
 
 class Orthographic(Projection):

lib/cartopy/tests/crs/test_stereographic.py

@@ -1,4 +1,4 @@
-# (C) British Crown Copyright 2013 - 2017, Met Office
+# (C) British Crown Copyright 2013 - 2018, Met Office
 #
 # This file is part of cartopy.
 #
@@ -52,14 +52,42 @@ def test_eccentric_globe(self):
                             [-3932.82587779, 3932.82587779], decimal=4)
 
     def test_true_scale(self):
-        # The "true_scale_latitude" parameter to Stereographic appears
-        # meaningless. This test just ensures that the correct proj4
-        # string is being created. (#339)
-        stereo = ccrs.Stereographic(true_scale_latitude=10)
-        expected = ('+ellps=WGS84 +proj=stere +lat_0=0.0 +lon_0=0.0 '
-                    '+x_0=0.0 +y_0=0.0 +lat_ts=10 +no_defs')
+        # The "true_scale_latitude" parameter only makes sense for
+        # polar stereographic projections (#339 and #455).
+        # For now only the proj4 string creation is tested
+        # See test_scale_factor for test on projection.
+        globe = ccrs.Globe(ellipse='sphere')
+        stereo = ccrs.NorthPolarStereo(true_scale_latitude=30, globe=globe)
+        expected = ('+ellps=sphere +proj=stere +lat_0=90 +lon_0=0.0 '
+                    '+x_0=0.0 +y_0=0.0 +lat_ts=30 +no_defs')
+        assert stereo.proj4_init == expected
+
+    def test_scale_factor(self):
+        # See #455
+        # Use spherical Earth in North Polar Stereographic to check
+        # equivalence between true_scale and scale_factor.
+        # In these conditions a scale factor of 0.75 corresponds exactly to
+        # a standard parallel of 30N.
+        globe = ccrs.Globe(ellipse='sphere')
+        stereo = ccrs.Stereographic(central_latitude=90., scale_factor=0.75,
+                                    globe=globe)
+        expected = ('+ellps=sphere +proj=stere +lat_0=90.0 +lon_0=0.0 '
+                    '+x_0=0.0 +y_0=0.0 +k_0=0.75 +no_defs')
         assert stereo.proj4_init == expected
 
+        # Now test projections
+        lon, lat = 10, 10
+        projected_scale_factor = stereo.transform_point(lon, lat,
+                                                        ccrs.Geodetic())
+
+        # should be equivalent to North Polar Stereo with
+        # true_scale_latitude = 30
+        nstereo = ccrs.NorthPolarStereo(globe=globe, true_scale_latitude=30)
+        projected_true_scale = nstereo.transform_point(lon, lat,
+                                                       ccrs.Geodetic())
+
+        assert projected_true_scale == projected_scale_factor
+
     def test_eastings(self):
         stereo = ccrs.Stereographic()
         stereo_offset = ccrs.Stereographic(false_easting=1234,