Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix limits/threshold of projections with non-default Globe #1607

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 27 additions & 6 deletions lib/cartopy/crs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree this would be nice.

a = np.float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS)
b = np.float(self.globe.semiminor_axis or a)
Comment on lines +845 to +846
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason for these to be np.float() instead of just leaving them as the values that are returned?


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):
Expand All @@ -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):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -2309,14 +2329,15 @@ 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):
return self._boundary

@property
def threshold(self):
return 1e5
return self._threshold

@property
def x_limits(self):
Expand Down
41 changes: 41 additions & 0 deletions lib/cartopy/tests/crs/test_transverse_mercator.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
import pytest

import cartopy.crs as ccrs
from .helpers import check_proj_params


@pytest.mark.parametrize('approx', [True, False])
Expand All @@ -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):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can get rid of this check now since we are dependent on pyproj's PROJ version and looks like we are at 8+
https://pyproj4.github.io/pyproj/3.3.1/

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,
Expand Down