From eeb2e19c84c43c7f6361b0855152b8fe43d4b0ad Mon Sep 17 00:00:00 2001 From: Elliott Sales de Andrade Date: Sun, 12 Jul 2020 03:34:45 -0400 Subject: [PATCH 1/3] Fix limits and threshold of TransverseMercator. It formerly did not scale with Globe radius, and so would produce a very very small result within a large boundary if using any smaller Globe radius. --- lib/cartopy/crs.py | 14 +++++-- .../tests/crs/test_transverse_mercator.py | 41 +++++++++++++++++++ 2 files changed, 52 insertions(+), 3 deletions(-) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 2fa1a2d90..ed3936756 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -841,9 +841,17 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, proj4_params += [('approx', None)] super().__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) + + self._x_limits = (-np.pi * a, np.pi * a) + self._y_limits = (-np.pi / 2 * b, np.pi / 2 * b) + self._threshold = min(a, b) * 1.5e-3 # Approximately 1e4 for defaults. + @property def threshold(self): - return 1e4 + return self._threshold @property def boundary(self): @@ -855,11 +863,11 @@ def boundary(self): @property def x_limits(self): - return (-2e7, 2e7) + return self._x_limits @property def y_limits(self): - return (-1e7, 1e7) + return self._y_limits class OSGB(TransverseMercator): diff --git a/lib/cartopy/tests/crs/test_transverse_mercator.py b/lib/cartopy/tests/crs/test_transverse_mercator.py index 3b35d038d..56540714f 100644 --- a/lib/cartopy/tests/crs/test_transverse_mercator.py +++ b/lib/cartopy/tests/crs/test_transverse_mercator.py @@ -12,6 +12,7 @@ import pytest import cartopy.crs as ccrs +from .helpers import check_proj_params @pytest.mark.parametrize('approx', [True, False]) @@ -21,16 +22,56 @@ def setup_class(self): self.point_b = (0.5, 50.5) self.src_crs = ccrs.PlateCarree() + def check_args(self, approx, proj, other_args): + if ccrs.PROJ4_VERSION < (6, 0, 0): + other_args = {*other_args, 'lon_0=0.0', 'lat_0=0.0', 'k=1.0', + 'x_0=0.0', 'y_0=0.0', 'units=m'} + check_proj_params('tmerc' if approx else 'etmerc', proj, + other_args) + else: + other_args = {*other_args, 'lon_0=0.0', 'lat_0=0.0', 'k=1.0', + 'x_0=0.0', 'y_0=0.0', 'units=m'} + if approx: + other_args.add('approx') + check_proj_params('tmerc', proj, other_args) + def test_default(self, approx): proj = ccrs.TransverseMercator(approx=approx) + self.check_args(approx, proj, {'ellps=WGS84'}) + + np.testing.assert_array_almost_equal(proj.x_limits, (-2e7, 2e7), + decimal=-5) + np.testing.assert_array_almost_equal(proj.y_limits, (-1e7, 1e7), + decimal=-5) + res = proj.transform_point(*self.point_a, src_crs=self.src_crs) np.testing.assert_array_almost_equal(res, (-245269.53181, 5627508.74355), decimal=5) + res = proj.transform_point(*self.point_b, src_crs=self.src_crs) np.testing.assert_array_almost_equal(res, (35474.63566645, 5596583.41949901)) + def test_sphere_globe(self, approx): + if not approx: + pytest.skip('Not supported by proj.') + + globe = ccrs.Globe(semimajor_axis=1000, ellipse=None) + proj = ccrs.TransverseMercator(approx=approx, globe=globe) + self.check_args(approx, proj, {'a=1000'}) + + np.testing.assert_array_almost_equal(proj.x_limits, + (-3141.592654, 3141.592654)) + np.testing.assert_array_almost_equal(proj.y_limits, + (-1570.796327, 1570.796327)) + + res = proj.transform_point(*self.point_a, src_crs=self.src_crs) + np.testing.assert_array_almost_equal(res, (-38.377488, 886.259630)) + + res = proj.transform_point(*self.point_b, src_crs=self.src_crs) + np.testing.assert_array_almost_equal(res, (5.550816, 881.409961)) + def test_osgb_vals(self, approx): proj = ccrs.TransverseMercator(central_longitude=-2, central_latitude=49, From 336c09fe5de9cc9a866778e264ee5b6bbdd42944 Mon Sep 17 00:00:00 2001 From: Elliott Sales de Andrade Date: Sun, 12 Jul 2020 04:45:07 -0400 Subject: [PATCH 2/3] Scale threshold to Globe radius in more projections. Fixes #1572, testing the example in the other projections as well. --- lib/cartopy/crs.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index ed3936756..aadc6cd9f 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -1649,9 +1649,15 @@ def __init__(self, central_longitude=0, false_easting=None, false_northing=false_northing, 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) + + self._threshold = min(a, b) / 63.78137 # About 1e5 for defaults. + @property def threshold(self): - return 1e5 + return self._threshold class EckertI(_Eckert): @@ -1821,9 +1827,15 @@ def __init__(self, central_longitude=0, globe=None, false_northing=false_northing, 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) + + self._threshold = min(a, b) / 63.78137 # About 1e5 for defaults. + @property def threshold(self): - return 1e5 + return self._threshold class Robinson(_WarpedRectangularProjection): @@ -2317,6 +2329,7 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, maxs = np.max(coords, axis=1) self._x_limits = mins[0], maxs[0] self._y_limits = mins[1], maxs[1] + self._threshold = min(a, b) / 63.78137 # About 1e5 for defaults. @property def boundary(self): @@ -2324,7 +2337,7 @@ def boundary(self): @property def threshold(self): - return 1e5 + return self._threshold @property def x_limits(self): From 82f9a8dbf4bf86de980ea9904270c98ddea0c8ff Mon Sep 17 00:00:00 2001 From: Elliott Sales de Andrade Date: Sun, 12 Jul 2020 20:53:09 -0400 Subject: [PATCH 3/3] Decrease threshold for AzimuthalEquidistant. Using the example from #1572 shows that the lines are jagged and the threshold is too high. --- lib/cartopy/crs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index aadc6cd9f..f193cbda2 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -2329,7 +2329,7 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, maxs = np.max(coords, axis=1) self._x_limits = mins[0], maxs[0] self._y_limits = mins[1], maxs[1] - self._threshold = min(a, b) / 63.78137 # About 1e5 for defaults. + self._threshold = min(a, b) * 1.5e-3 # Approximately 1e4 for defaults. @property def boundary(self):