From 5c7e3326ba7feb6310a9b017aaa8dac082578e41 Mon Sep 17 00:00:00 2001 From: Ryan May Date: Wed, 28 Feb 2018 12:28:16 -0700 Subject: [PATCH] Add more options to Mercator projection This adds scale_factor as an additional way to control the scaling (in addition to latitude_true_scale). It also adds support for false_easting and false_northing. These options are all required to be fully capable of representing Mercator projection coming from a CF-compliant netCDF file. --- lib/cartopy/crs.py | 27 +++++++++++++++-- lib/cartopy/tests/crs/test_mercator.py | 40 +++++++++++++++++++++----- 2 files changed, 58 insertions(+), 9 deletions(-) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 57fb8c061..9b8d16441 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -943,7 +943,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 ---------- @@ -959,12 +960,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. diff --git a/lib/cartopy/tests/crs/test_mercator.py b/lib/cartopy/tests/crs/test_mercator.py index 522d5e9c2..fa920915c 100644 --- a/lib/cartopy/tests/crs/test_mercator.py +++ b/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)