diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 2fa1a2d90..f193cbda2 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): @@ -1641,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): @@ -1813,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): @@ -2309,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) * 1.5e-3 # Approximately 1e4 for defaults. @property def boundary(self): @@ -2316,7 +2337,7 @@ def boundary(self): @property def threshold(self): - return 1e5 + return self._threshold @property def x_limits(self): 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,