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

Add ability to specify custom bounds to Projections #1888

Open
wants to merge 5 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
45 changes: 35 additions & 10 deletions lib/cartopy/crs.py
Original file line number Diff line number Diff line change
Expand Up @@ -644,27 +644,52 @@ class Projection(CRS, metaclass=ABCMeta):
# Whether or not this projection can handle wrapped coordinates
_wrappable = False

def __init__(self, *args, **kwargs):
def __init__(self, *args, bounds=None, transform_bounds=False, **kwargs):
"""
Parameters
----------
bounds: list or tuple, optional
Four element projection bounds (xmin, xmax, ymin, ymax).
If not provided defaults to "area_of_use" of the Coordinate
Reference System if it exists. This can be used to define an area
of use when one doesn't exist or to define a subset of the full
projection space.
transform_bounds: bool, optional
Transform bounds provided as lon(x)/lat(y) degrees to projection
space. This is a convenience when bounds are only known in degrees
for a non-geographic CRS. Defaults to False.

Additional arguments are passed directly to :class:`cartopy.crs.CRS`.

"""
super().__init__(*args, **kwargs)
self.bounds = None
if self.area_of_use:
if bounds is None and self.area_of_use is not None:
bounds = (
self.area_of_use.west,
self.area_of_use.east,
self.area_of_use.south,
self.area_of_use.north,
)
transform_bounds = True

if bounds is not None and transform_bounds:
# Convert lat/lon bounds to projected bounds.
# Geographic area of the entire dataset referenced to WGS 84
# NB. We can't use a polygon transform at this stage because
# that relies on the existence of the map boundary... the very
# thing we're trying to work out! ;-)
x0 = self.area_of_use.west
x1 = self.area_of_use.east
y0 = self.area_of_use.south
y1 = self.area_of_use.north
x0, x1, y0, y1 = bounds
lons = np.array([x0, x0, x1, x1])
lats = np.array([y0, y1, y1, y0])
points = self.transform_points(self.as_geodetic(), lons, lats)
x = points[:, 0]
y = points[:, 1]
self.bounds = (x.min(), x.max(), y.min(), y.max())
x0, x1, y0, y1 = self.bounds
self.threshold = min(x1 - x0, y1 - y0) / 100.
bounds = (x.min(), x.max(), y.min(), y.max())

self.bounds = None
if bounds is not None:
x0, x1, y0, y1 = self.bounds = tuple(bounds)
self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.

@property
def boundary(self):
Expand Down
35 changes: 35 additions & 0 deletions lib/cartopy/tests/test_crs.py
Original file line number Diff line number Diff line change
Expand Up @@ -332,3 +332,38 @@ def test_projection__from_string():

def test_crs__from_pyproj_crs():
assert ccrs.CRS(pyproj.CRS("EPSG:4326")) == "EPSG:4326"


@pytest.mark.parametrize(
("crs_input", "bounds", "transform_bounds", "exp_bounds", "exp_threshold"),
[
("EPSG:4326", None, False,
(-180.0, 180.0, -90.0, 90.0), 1.8),
("EPSG:4326", [-100.0, -90.0, 35.0, 45.0], False,
(-100.0, -90.0, 35.0, 45.0), 0.1),
("EPSG:4326", [-100.0, -90.0, 35.0, 45.0], True,
(-100.0, -90.0, 35.0, 45.0), 0.1),
("EPSG:6932", [-40000., 40000., -40000., 40000.], False,
(-40000., 40000., -40000., 40000.), 800.0),
("EPSG:6932", [20.0, 25.0, -85.0, -78.0], True,
(190942.4609561831, 565329.6096711607,
505972.06807129766, 1257011.612161974), 3743.871487149776),
# Meteosat Second Generation (MSG) - SEVIRI - Flipped Geostationary
({'proj': 'geos', 'lon_0': 0.0, 'a': 6378169.00, 'b': 6356583.80,
'h': 35785831.00, 'units': 'm'},
[5500000, -5500000, -5500000, 5500000], False,
(5500000, -5500000, -5500000, 5500000), 110000.0),
({'proj': 'geos', 'lon_0': 0.0, 'a': 6378169.00, 'b': 6356583.80,
'h': 35785831.00, 'units': 'm'},
[-10.0, 10.0, -10.0, 10.0], True,
(-1084697.7494802547, 1084697.7494802547,
-1093480.6233566382, 1093480.6233566382), 21693.954989605096),
]
)
def test_projection_with_bounds(crs_input, bounds, transform_bounds,
exp_bounds, exp_threshold):
pcrs = pyproj.CRS.from_user_input(crs_input)
crs = ccrs.Projection(pcrs, bounds=bounds,
transform_bounds=transform_bounds)
assert crs.bounds == exp_bounds
assert crs.threshold == exp_threshold