Merge pull request #1043 from dopplershift/mercator-options

Add more options to Mercator projection

lib/cartopy/crs.py

@@ -945,7 +945,8 @@ class Mercator(Projection):
 
     def __init__(self, central_longitude=0.0,
                  min_latitude=-80.0, max_latitude=84.0,
-                 globe=None, latitude_true_scale=0.0):
+                 globe=None, latitude_true_scale=None,
+                 false_easting=0.0, false_northing=0.0, scale_factor=None):
         """
         Parameters
         ----------
@@ -961,12 +962,34 @@ def __init__(self, central_longitude=0.0,
             If omitted, a default globe is created.
         latitude_true_scale: optional
             The latitude where the scale is 1. Defaults to 0 degrees.
+        false_easting: optional
+            X offset from the planar origin in metres. Defaults to 0.
+        false_northing: optional
+            Y offset from the planar origin in metres. Defaults to 0.
+        scale_factor: optional
+            Scale factor at natural origin. Defaults to unused.
 
+        Only one of ``latitude_true_scale`` and ``scale_factor`` should
+        be included.
         """
         proj4_params = [('proj', 'merc'),
                         ('lon_0', central_longitude),
-                        ('lat_ts', latitude_true_scale),
+                        ('x_0', false_easting),
+                        ('y_0', false_northing),
                         ('units', 'm')]
+
+        # If it's None, we don't pass it to Proj4, in which case its default
+        # of 0.0 will be used.
+        if latitude_true_scale is not None:
+            proj4_params.append(('lat_ts', latitude_true_scale))
+
+        if scale_factor is not None:
+            if latitude_true_scale is not None:
+                raise ValueError('It does not make sense to provide both '
+                                 '"scale_factor" and "latitude_true_scale". ')
+            else:
+                proj4_params.append(('k_0', scale_factor))
+
         super(Mercator, self).__init__(proj4_params, globe=globe)
 
         # Calculate limits.

lib/cartopy/tests/crs/test_mercator.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.
 #
@@ -25,8 +25,8 @@
 def test_default():
     crs = ccrs.Mercator()
 
-    assert crs.proj4_init == ('+ellps=WGS84 +proj=merc +lon_0=0.0 +lat_ts=0.0 '
-                              '+units=m +no_defs')
+    assert crs.proj4_init == ('+ellps=WGS84 +proj=merc +lon_0=0.0 +x_0=0.0 '
+                              '+y_0=0.0 +units=m +no_defs')
     assert_almost_equal(crs.boundary.bounds,
                         [-20037508, -15496571, 20037508, 18764656], decimal=0)
 
@@ -36,7 +36,7 @@ def test_eccentric_globe():
                        ellipse=None)
     crs = ccrs.Mercator(globe=globe, min_latitude=-40, max_latitude=40)
     assert crs.proj4_init == ('+a=10000 +b=5000 +proj=merc +lon_0=0.0 '
-                              '+lat_ts=0.0 +units=m +no_defs')
+                              '+x_0=0.0 +y_0=0.0 +units=m +no_defs')
 
     assert_almost_equal(crs.boundary.bounds,
                         [-31415.93, -2190.5, 31415.93, 2190.5], decimal=2)
@@ -60,7 +60,7 @@ def test_equality():
 def test_central_longitude():
     cl = 10.0
     crs = ccrs.Mercator(central_longitude=cl)
-    proj4_str = ('+ellps=WGS84 +proj=merc +lon_0={} +lat_ts=0.0 '
+    proj4_str = ('+ellps=WGS84 +proj=merc +lon_0={} +x_0=0.0 +y_0=0.0 '
                  '+units=m +no_defs'.format(cl))
     assert crs.proj4_init == proj4_str
 
@@ -71,9 +71,35 @@ def test_central_longitude():
 def test_latitude_true_scale():
     lat_ts = 20.0
     crs = ccrs.Mercator(latitude_true_scale=lat_ts)
-    proj4_str = ('+ellps=WGS84 +proj=merc +lon_0=0.0 +lat_ts={} '
-                 '+units=m +no_defs'.format(lat_ts))
+    proj4_str = ('+ellps=WGS84 +proj=merc +lon_0=0.0 +x_0=0.0 +y_0=0.0 '
+                 '+units=m +lat_ts={} +no_defs'.format(lat_ts))
     assert crs.proj4_init == proj4_str
 
     assert_almost_equal(crs.boundary.bounds,
                         [-18836475, -14567718, 18836475, 17639917], decimal=0)
+
+
+def test_easting_northing():
+    false_easting = 1000000
+    false_northing = -2000000
+    crs = ccrs.Mercator(false_easting=false_easting,
+                        false_northing=false_northing)
+    proj4_str = ('+ellps=WGS84 +proj=merc +lon_0=0.0 +x_0={} +y_0={} '
+                 '+units=m +no_defs'.format(false_easting, false_northing))
+    assert crs.proj4_init == proj4_str
+
+    assert_almost_equal(crs.boundary.bounds,
+                        [-19037508, -17496571, 21037508, 16764656], decimal=0)
+
+
+def test_scale_factor():
+    # Should be same as lat_ts=20 for a sphere
+    scale_factor = 0.939692620786
+    crs = ccrs.Mercator(scale_factor=scale_factor,
+                        globe=ccrs.Globe(ellipse='sphere'))
+    proj4_str = ('+ellps=sphere +proj=merc +lon_0=0.0 +x_0=0.0 +y_0=0.0 '
+                 '+units=m +k_0={:.12f} +no_defs'.format(scale_factor))
+    assert crs.proj4_init == proj4_str
+
+    assert_almost_equal(crs.boundary.bounds,
+                        [-18808021, -14585266, 18808021, 17653216], decimal=0)