diff --git a/pyresample/future/spherical/__init__.py b/pyresample/future/spherical/__init__.py new file mode 100644 index 000000000..b21f05de9 --- /dev/null +++ b/pyresample/future/spherical/__init__.py @@ -0,0 +1,27 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (c) 2021 Pyresample developers +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License along +# with this program. If not, see . +"""Future features that are backwards incompatible with current functionality.""" + +from .extent import SExtent # noqa +from .point import SMultiPoint, SPoint # noqa + +__all__ = [ + "SExtent", + "SPoint", + "SMultiPoint", +] diff --git a/pyresample/future/spherical/extent.py b/pyresample/future/spherical/extent.py new file mode 100644 index 000000000..873f3c656 --- /dev/null +++ b/pyresample/future/spherical/extent.py @@ -0,0 +1,297 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Copyright (c) 2013-2022 Pyresample Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +"""Definition of SExtent class.""" +import numbers + +import numpy as np +from shapely.geometry import MultiPolygon, Polygon +from shapely.ops import unary_union + + +def bounds_from_extent(extent): + """Get shapely bounds from a matplotlib/cartopy extent. + + Shapely bounds: [x_min, y_min, x_max, y_max] + Matplotlib/Cartopy extent: [x_min, x_max, y_min, y_max] + """ + bounds = [extent[0], extent[2], extent[1], extent[3]] + return bounds + + +def _get_list_extents_from_args(args): + """Get list_extent from SExtent input arguments.""" + args = list(args) + # Check something is passed to SExtent + if len(args) == 0: + raise ValueError("No extent passed to SExtent.") + # Check that the first args element is a list, tuple of np.array + if np.any([not isinstance(arg, (list, tuple, np.ndarray)) for arg in args]): + raise TypeError("Specify the extent(s) using list, tuple or np.array.") + # Get list extent + # - If SExtent(extent) or SExtent(list_extent) + if len(args) == 1: + list_extents = list(*args) + if len(list_extents) == 0: # input case: [], or () + raise ValueError("No extent passed to SExtent.") + else: + list_extents = args + # If a single list, check correct length and contain a number + if len(list_extents) == 4 and isinstance(list_extents[0], numbers.Number): + list_extents = [list_extents] + # Ensure that here after is a list of extents + if not isinstance(list_extents[0], (tuple, list, np.ndarray)): + raise ValueError("SExtent expects a single extent or a list of extents.") + return list_extents + + +def _check_extent_order(extent): + """Check extent order validity.""" + if extent[0] > extent[1]: + raise ValueError('extent[0] (lon_min) must be smaller than extent[1] (lon_max).') + if extent[2] > extent[3]: + raise ValueError('extent[2] (lat_min) must be smaller than extent[3] (lat_max).') + + +def _check_extent_values(extent, use_radians=False): + """Check extent value validity.""" + lat_valid_range = [-90.0, 90.0] + lon_valid_range = [-180.0, 180.0] + extent = np.array(extent) + if use_radians: + lat_valid_range = np.deg2rad(lat_valid_range) + lon_valid_range = np.deg2rad(lon_valid_range) + if np.any(np.logical_or(extent[0:2] < lon_valid_range[0], extent[0:2] > lon_valid_range[1])): + raise ValueError(f'extent longitude values must be within {lon_valid_range}.') + if np.any(np.logical_or(extent[2:] < lat_valid_range[0], extent[2:] > lat_valid_range[1])): + raise ValueError(f'extent latitude values must be within {lat_valid_range}.') + + +def _check_is_not_point_extent(extent): + """Check the extent does not represent a point.""" + if extent[0] == extent[1] and extent[2] == extent[3]: + raise ValueError("An extent can not be defined by a point.") + + +def _check_is_not_line_extent(extent): + """Check the extent does not represent a line.""" + if extent[0] == extent[1] or extent[2] == extent[3]: + raise ValueError("An extent can not be defined by a line.") + + +def _check_valid_extent(extent, use_radians=False): + """Check lat/lon extent validity.""" + # Check extent length + if len(extent) != 4: + raise ValueError("'extent' must have length 4: [lon_min, lon_max, lat_min, lat_max].") + # Convert to np.array + extent = np.array(extent) + # Check order and values + _check_extent_order(extent) + _check_extent_values(extent, use_radians=use_radians) + # Check is not point or line + _check_is_not_point_extent(extent) + _check_is_not_line_extent(extent) + # Ensure is a list + extent = extent.tolist() + return extent + + +def _check_non_overlapping_extents(polygons): + """Check that the extents polygons are non-overlapping. + + Note: Touching extents are considered valids. + """ + for poly in polygons.geoms: + # Intersects includes within/contain/equals/touches ! + n_intersects = np.sum([p.intersects(poly) and not p.touches(poly) for p in polygons.geoms]) + if n_intersects >= 2: + raise ValueError("The extents composing SExtent can not be duplicates or overlapping.") + + +def _is_global_extent(polygons): + """Check if a list of extents composes a global extent.""" + from shapely.geometry import Polygon + from shapely.ops import unary_union + + # Try union the polygons + unioned_polygon = unary_union(polygons) + # If still a MultiPolygon, not a global extent + if not isinstance(unioned_polygon, Polygon): + return False + # Define global polygon + global_polygon = Polygon.from_bounds(-180, -90, 180, 90) + # Assess if global + is_global = global_polygon.equals(unioned_polygon) + return is_global + + +class SExtent: + """Spherical Extent. + + SExtent longitudes are defined between -180 and 180 degree. + Geometrical operations are performed on a planar coordinate system. + + A spherical geometry crossing the anti-meridian will have an SExtent + composed of [lon_start, 180, ...,...] and [-180, lon_end, ..., ...] + + Important notes: + - SExtents touching at the anti-meridian are considered to not touching ! + - SExtents.intersects predicate includes also contains/within cases. + - In comparison to shapely polygon behaviour, SExtent.intersects is + False if SExtents are just touching. + + There is not an upper limit on the number of extents composing SExtent !!! + The only conditions is that the extents do not intersect/overlap each other. + + More specifically, extents composing an SExtent: + - must represent a polygon (not a point or a line), + - can touch each other, + - but can not intersect/contain/overlap each other. + + Examples of valid SExtent definitions: + extent = [x_min, x_max, y_min, ymax] + sext = SExtent(extent) + sext = SExtent([extent, extent]) + sext = SExtent(extent, extent, ...) + + """ + + def __init__(self, *args): + # Get list of extents + list_extents = _get_list_extents_from_args(args) + # Check extents format and value validity + self._extents = [_check_valid_extent(ext) for ext in list_extents] + # Pre-compute shapely polygon + list_polygons = [Polygon.from_bounds(*bounds_from_extent(ext)) for ext in self._extents] + self._polygons = MultiPolygon(list_polygons) + # Check topological validity of extents polygons + _check_non_overlapping_extents(self._polygons) + + def to_shapely(self): + """Return the shapely extent rectangle(s) polygon(s).""" + return self._polygons + + def __str__(self): + """Get simplified representation of SExtent.""" + return str(self._extents) + + def __repr__(self): + """Get simplified representation of SExtent.""" + return str(self._extents) + + def __iter__(self): + """Get extents iterator.""" + return self._extents.__iter__() + + def _repr_svg_(self): + """Display the SExtent in the Ipython terminal.""" + return self.to_shapely()._repr_svg_() + + @property + def is_global(self): + """Check if the extent is global.""" + # Try union the polygons + unioned_polygon = unary_union(self._polygons) + # If still a MultiPolygon, not a global extent + if not isinstance(unioned_polygon, Polygon): + return False + # Define global polygon + global_polygon = Polygon.from_bounds(-180.0, -90.0, 180.0, 90.0) + # Assess if global + is_global = global_polygon.equals(unioned_polygon) + return is_global + + def intersects(self, other): + """Check if SExtent is intersecting the other SExtent. + + Touching SExtent are considered to not intersect ! + SExtents intersection also includes contains/within occurence. + """ + if not isinstance(other, SExtent): + raise TypeError("SExtent.intersects() expects a SExtent class instance.") + # Important caveat in comparison to shapely ! + # In shapely, intersects includes touching geometries ! + bl = (self.to_shapely().intersects(other.to_shapely()) and + not self.to_shapely().touches(other.to_shapely())) + return bl + + def disjoint(self, other): + """Check if SExtent does not intersect (and do not touch) the other SExtent.""" + if not isinstance(other, SExtent): + raise TypeError("SExtent.disjoint() expects a SExtent class instance.") + return self.to_shapely().disjoint(other.to_shapely()) + + def within(self, other): + """Check if SExtent is within another SExtent.""" + if not isinstance(other, SExtent): + raise TypeError("SExtent.within() expects a SExtent class instance.") + return self.to_shapely().within(other.to_shapely()) + + def contains(self, other): + """Check if the SExtent contains another SExtent.""" + if not isinstance(other, SExtent): + raise TypeError("SExtent.contains() expects a SExtent class instance.") + return self.to_shapely().contains(other.to_shapely()) + + def touches(self, other): + """Check if SExtent external touches another SExtent. + + Important notes: + - Touching SExtents are not disjoint ! + - Touching SExtents are considered to not intersect ! + - SExtents which are contained/within each other AND are "touching in the interior" + are not considered touching ! + """ + if not isinstance(other, SExtent): + raise TypeError("SExtent.touches() expects a SExtent class instance.") + return self.to_shapely().touches(other.to_shapely()) + + def equals(self, other): + """Check if SExtent is equals to SExtent.""" + if not isinstance(other, SExtent): + raise TypeError("SExtent.equals() expects a SExtent class instance.") + return self.to_shapely().equals(other.to_shapely()) + + def plot(self, ax=None, **plot_kwargs): + """Plot the SLine using Cartopy.""" + import matplotlib.pyplot as plt + try: + import cartopy.crs as ccrs + except ModuleNotFoundError: + raise ModuleNotFoundError("Install cartopy to plot spherical geometries.") + + # Retrieve shapely polygon + geom = self.to_shapely() + + # Create figure if ax not provided + ax_not_provided = False + if ax is None: + ax_not_provided = True + proj_crs = ccrs.PlateCarree() + fig, ax = plt.subplots(subplot_kw=dict(projection=proj_crs)) + + # Add extent polygon + ax.add_geometries([geom], crs=ccrs.PlateCarree(), **plot_kwargs) + # Beautify plot by default + if ax_not_provided: + ax.stock_img() + ax.coastlines() + gl = ax.gridlines(draw_labels=True, linestyle='--') + gl.xlabels_top = False + gl.ylabels_right = False + return ax diff --git a/pyresample/test/test_sgeom/__init__.py b/pyresample/test/test_sgeom/__init__.py new file mode 100644 index 000000000..b496d5604 --- /dev/null +++ b/pyresample/test/test_sgeom/__init__.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (c) 2021 Pyresample developers +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License along +# with this program. If not, see . +"""Test for spherical geometries.""" diff --git a/pyresample/test/test_sgeom/test_extent.py b/pyresample/test/test_sgeom/test_extent.py new file mode 100644 index 000000000..cd535009f --- /dev/null +++ b/pyresample/test/test_sgeom/test_extent.py @@ -0,0 +1,589 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Copyright (c) 2013-2022 Pyresample Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +"""Test cases for SExtent class.""" +import numpy as np +import pytest + +from pyresample.future.spherical import SExtent + + +class TestSExtent: + """Test SExtent.""" + + def test_sextent_with_bad_value(self): + """Test it raise error when providing bad values to SExtent.""" + # lon_min > lon_max + extent = [10, 0, -40, 40] + with pytest.raises(ValueError): + SExtent(extent) + + # lat_min > lat_max + extent = [-10, 0, 41, 40] + with pytest.raises(ValueError): + SExtent(extent) + + # lon_min < -180 + extent = [-181, 0, -40, 40] + with pytest.raises(ValueError): + SExtent(extent) + + # lon_max > 180 + extent = [-180, 181, -40, 40] + with pytest.raises(ValueError): + SExtent(extent) + + # lat_min < -90 + extent = [50, 60, -91, 40] + with pytest.raises(ValueError): + SExtent(extent) + + # lat_max > 90 + extent = [50, 60, -40, 91] + with pytest.raises(ValueError): + SExtent(extent) + + # Passing a list of extent with bad values + list_extent = [[-180, 181, -40, 91], [-180, 181, -40, 91]] + with pytest.raises(ValueError): + SExtent(list_extent) + + def test_sextent_with_bad_input(self): + """Test it raise error when providing bad inputs to SExtent.""" + # Length 5 + extent = [-180, 181, -40, 91, 50] + with pytest.raises(ValueError): + SExtent(extent) + # Length 1 + extent = [-180] + with pytest.raises(ValueError): + SExtent(extent) + # 4 single digits + with pytest.raises(TypeError): + SExtent(1, 2, 3, 4) + # Length 0 + extent = [] + with pytest.raises(ValueError): + SExtent(extent) + # No argument + with pytest.raises(ValueError): + SExtent() + # None + with pytest.raises(TypeError): + SExtent(None) + # Passing a list of extent with bad length + list_extent = [[-180, 180, -40, 90, 50], [-180, 180, -40, 90]] + with pytest.raises(ValueError): + SExtent(list_extent) + + def test_sextent_with_correct_format(self): + """Test SExtent when providing correct extent(s) format and values.""" + # Accept list + extent = [-180, -175, -40, 40] + assert list(SExtent(extent))[0] == extent + # Accept tuple + extent = (-180, -175, -40, 40) + assert list(SExtent(extent))[0] == list(extent) + # Accept numpy array + extent = np.array([-180, -175, -40, 40]) + assert list(SExtent(extent))[0] == extent.tolist() + + # Accept list of extents (list) + extent1 = [50, 60, -40, 40] + extent2 = [175, 180, -40, 40] + list_extent = [extent1, extent2] + assert np.allclose(list(SExtent(list_extent)), list_extent) + # Accept list of extents (tuple) + extent1 = (50, 60, -40, 40) + extent2 = (175, 180, -40, 40) + list_extent = [extent1, extent2] + assert np.allclose(list(SExtent(list_extent)), list_extent) + # Accept list of extents (np.array) + extent1 = np.array([0, 60, -40, 40]) + extent2 = np.array([175, 180, -40, 40]) + list_extent = [extent1, extent2] + assert np.allclose(list(SExtent(list_extent)), list_extent) + # Accept multiple extents (list) + extent1 = [50, 60, -40, 40] + extent2 = [175, 180, -40, 40] + list_extent = [extent1, extent2] + assert np.allclose(list(SExtent(extent1, extent2)), list_extent) + # Accept multiple extents (tuple) + extent1 = (50, 60, -40, 40) + extent2 = (175, 180, -40, 40) + list_extent = [extent1, extent2] + assert np.allclose(list(SExtent(extent1, extent2)), list_extent) + # Accept multiple extents (np.array) + extent1 = np.array([0, 60, -40, 40]) + extent2 = np.array([175, 180, -40, 40]) + list_extent = [extent1, extent2] + assert np.allclose(list(SExtent(extent1, extent2)), list_extent) + + def test_single_sextent_bad_topology(self): + """Test that raise error when the extents is a point or a line.""" + # - Point extent + extent = [0, 0, 40, 40] + with pytest.raises(ValueError): + SExtent(extent) + # - Line extent + extent = [0, 10, 40, 40] + with pytest.raises(ValueError): + SExtent(extent) + extent = [0, 0, -40, 50] + with pytest.raises(ValueError): + SExtent(extent) + + def test_multiple_touching_extents(self): + """Test that touching extents composing SExtent do not raise error.""" + extent1 = [0, 40, 0, 40] + extent2 = [0, 40, -40, 0] + _ = SExtent(extent1, extent2) + + def test_multiple_overlapping_extents(self): + """Test that raise error when the extents composing SExtent overlaps.""" + # Intersecting raise error + extent1 = [0, 40, 0, 40] + extent2 = [20, 60, 20, 60] + with pytest.raises(ValueError): + SExtent(extent1, extent2) + + # Duplicate extent raise errors + extent1 = [0, 40, 0, 40] + with pytest.raises(ValueError): + SExtent(extent1, extent1) + + # Within extent raise errors + extent1 = [0, 40, 0, 40] + extent2 = [10, 20, 10, 20] + with pytest.raises(ValueError): + SExtent(extent1, extent2) + + def test_to_shapely(self): + """Test shapely conversion.""" + from shapely.geometry import MultiPolygon, Polygon + extent = [0, 20, 10, 30] + bounds = [0, 10, 20, 30] + shapely_sext = SExtent(extent).to_shapely() + shapely_polygon = MultiPolygon([Polygon.from_bounds(*bounds)]) + assert shapely_sext.equals(shapely_polygon) + + def test_iter(self): + """Test extent list extraction.""" + # Check single extent + extent = [0, 20, 10, 30] + sext = SExtent(extent) + list_extents = list(sext) + assert list_extents == [extent] + + # Check multiple extents + extent1 = [0, 20, 10, 30] + extent2 = [0, 10, 40, 60] + sext = SExtent(extent1, extent2) + list_extents = list(sext) + assert list_extents == [extent1, extent2] + + def test_str(self): + """Check the string representation.""" + extent = [0, 20, 10, 30] + sext = SExtent(extent) + assert str(sext) == '[[0, 20, 10, 30]]' + + def test_repr(self): + """Check the representation.""" + extent = [0, 20, 10, 30] + sext = SExtent(extent) + assert repr(sext) == '[[0, 20, 10, 30]]' + + def test_repr_svg(self): + """Check the sketch Ipython representation.""" + extent = [-180, 180, -90, 90] + sext = SExtent(extent) + repr_svg = '' + \ + '' + \ + '' + assert sext._repr_svg_() == repr_svg + + def test_is_global(self): + """Check is_global property.""" + # Is global + extent = [-180, 180, -90, 90] + sext = SExtent(extent) + assert sext.is_global + + # Is clearly not global (single extent) + extent = [-175, 180, -90, 90] + sext = SExtent(extent) + assert not sext.is_global + + # Is clearly not global (multiple extent) + extent1 = [-170, -160, -90, 90] + extent2 = [0, 10, -90, 90] + sext = SExtent(extent1, extent2) + assert not sext.is_global + + # Is global, but with multiple extents + extent1 = [-180, 0, -90, 90] + extent2 = [0, 180, -90, 90] + sext = SExtent(extent1, extent2) + assert sext.is_global + + # Is not global, but polgon bounds ... + extent1 = [-180, 0, -90, 90] + extent2 = [0, 180, 0, 90] + sext = SExtent(extent1, extent2) + assert not sext.is_global + + def test_SExtent_single_not_intersect(self): + """Check disjoint extents.""" + # - Not intersecting across longitude + extent1 = [50, 60, -40, 40] + extent2 = [175, 180, -40, 40] + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.disjoint(sext2) + assert sext2.disjoint(sext1) + + assert not sext1.intersects(sext2) + assert not sext2.intersects(sext1) + + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + # - Not intersecting across longitude (at the antimeridian) + extent1 = [175, 180, -40, 40] + extent2 = [-180, -175, -40, 40] + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.disjoint(sext2) + assert sext2.disjoint(sext1) + + assert not sext1.intersects(sext2) + assert not sext2.intersects(sext1) + + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + # - Not intersecting across latitude + extent1 = [50, 60, -50, -40] + extent2 = [50, 60, 0, 40] + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.disjoint(sext2) + assert sext2.disjoint(sext1) + + assert not sext1.intersects(sext2) + assert not sext2.intersects(sext1) + + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + def test_SExtent_single_touches(self): + """Check touching extents.""" + extent1 = [0, 10, 40, 60] + extent2 = [0, 10, 30, 40] + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.touches(sext2) + assert sext2.touches(sext1) + + # Touching extents are not disjoint ! + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + # Touching extents do no intersect ! + assert not sext1.intersects(sext2) + assert not sext2.intersects(sext1) + + # Touching extents does not contain or are within + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + def test_SExtent_single_intersect(self): + """Check intersecting extents.""" + extent1 = [0, 10, 40, 60] + extent2 = [0, 20, 30, 50] + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.intersects(sext2) + assert sext2.intersects(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + def test_SExtent_single_untouching_within_contains(self): + """Check untouching extents which contains/are within each other.""" + extent1 = [0, 40, 0, 40] # contains extent2 + extent2 = [10, 20, 10, 20] # is within extent1 + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.intersects(sext2) + assert sext2.intersects(sext1) + + assert sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert sext2.within(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + def test_SExtent_single_touching_within_contains(self): + """Check touching extents which contains/are within each other.""" + extent1 = [0, 40, 0, 40] # contains extent2 + extent2 = [10, 40, 10, 40] # is within extent1 + sext1 = SExtent(extent1) + sext2 = SExtent(extent2) + + assert sext1.intersects(sext2) + assert sext2.intersects(sext1) + + assert sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert sext2.within(sext1) + + # Although they touch interiorly, touches is False ! + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + def test_SExtent_multiple_not_intersect(self): + """Check non-intersecting extents.""" + extent11 = [50, 60, -40, 40] + extent21 = [-180, -170, -40, 40] + extent22 = [90, 100, -40, 40] + sext1 = SExtent(extent11) + sext2 = SExtent(extent21, extent22) + + assert sext1.disjoint(sext2) + assert sext2.disjoint(sext1) + + assert not sext1.intersects(sext2) + assert not sext2.intersects(sext1) + + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + def test_SExtent_multiple_touches(self): + """Check touching extents.""" + # - Touches (not crossing the interior) + extent11 = [50, 60, -40, 40] + extent21 = [60, 70, -40, 40] + extent22 = [-175, -170, -40, 40] + sext1 = SExtent(extent11) + sext2 = SExtent(extent21, extent22) + + assert sext1.touches(sext2) + assert sext2.touches(sext1) + + # Touching extents do no intersect ! + assert not sext1.intersects(sext2) + assert not sext2.intersects(sext1) + + # Touching extents are not disjoint ! + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + # Touching extents does not contain or are within + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + def test_SExtent_multiple_intersect(self): + """Check intersecting extents.""" + # - Intersect 1 + extent11 = [50, 60, -40, 40] + extent21 = [55, 70, -40, 40] + extent22 = [-175, -170, -40, 40] + sext1 = SExtent(extent11) + sext2 = SExtent(extent21, extent22) + + assert sext1.intersects(sext2) + assert sext2.intersects(sext1) + + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + assert not sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert not sext2.within(sext1) + + def test_SExtent_multiple_within(self): + """Check touching extents which contains/are within each other.""" + extent11 = [-180, 180, -40, 40] # containts extent 2 + extent21 = [55, 70, -40, 40] # within extent 1 + extent22 = [-180, -170, -40, 40] # within extent 1 + sext1 = SExtent(extent11) + sext2 = SExtent(extent21, extent22) + + assert sext2.intersects(sext1) + assert sext1.intersects(sext2) + + assert sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert sext2.within(sext1) + + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + # Touches is False because touches in the interior ! + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + def test_SExtent_multiple_within_at_antimeridian(self): + """Check touching extents which contains/are within each other.""" + # 2 within 1 + extent11 = [160, 180, -40, 40] # containts extent 2 + extent12 = [-180, -160, -40, 40] # containts extent 2 + extent21 = [170, 180, -20, 20] # within extent 1 + extent22 = [-180, -170, -20, 20] # within extent 1 + sext1 = SExtent(extent11, extent12) + sext2 = SExtent(extent21, extent22) + + assert sext1.intersects(sext2) + assert sext2.intersects(sext1) + + assert sext1.contains(sext2) + assert not sext2.contains(sext1) + + assert not sext1.within(sext2) + assert sext2.within(sext1) + + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + # Touches is False because touches in the interior ! + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + def test_SExtent_multiple_equals(self): + """Check extent equality.""" + extent11 = [160, 180, -40, 40] + extent12 = [-180, -160, -40, 40] + extent21 = [160, 170, -40, 40] + extent22 = [170, 180, -40, 40] + extent23 = [-180, -160, -40, 40] + sext1 = SExtent(extent11, extent12) + sext2 = SExtent(extent21, extent22, extent23) + + assert sext1.equals(sext2) + assert sext2.equals(sext1) + + # Equal extents intersects, are within/contains each other + assert sext1.intersects(sext2) + assert sext2.intersects(sext1) + + assert sext1.within(sext2) + assert sext2.within(sext1) + + assert sext1.contains(sext2) + assert sext2.contains(sext1) + + # Equal extents are not disjoint ! + assert not sext1.disjoint(sext2) + assert not sext2.disjoint(sext1) + + # Equal extents do not touches ! + assert not sext1.touches(sext2) + assert not sext2.touches(sext1) + + def test_SExtent_binary_predicates_bad_input(self): + extent = [160, 180, -40, 40] + sext = SExtent(extent) + # Disjoint + with pytest.raises(TypeError): + sext.disjoint("bad_dtype") + # Intersects + with pytest.raises(TypeError): + sext.intersects("bad_dtype") + # Within + with pytest.raises(TypeError): + sext.within("bad_dtype") + # Contains + with pytest.raises(TypeError): + sext.contains("bad_dtype") + # Touches + with pytest.raises(TypeError): + sext.touches("bad_dtype") + # Equals + with pytest.raises(TypeError): + sext.equals("bad_dtype") diff --git a/pyresample/test/test_sgeom/test_shapely.py b/pyresample/test/test_sgeom/test_shapely.py new file mode 100644 index 000000000..49075a68f --- /dev/null +++ b/pyresample/test/test_sgeom/test_shapely.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Copyright (c) 2013-2022 Pyresample Developers +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +"""Test cases for shapely.""" +import numpy as np +from shapely.geometry import MultiPolygon, Polygon +from shapely.ops import unary_union + + +class TestShapely: + """Test Shapely methods.""" + + def test_polygon_unary_union(self): + """Test MultiPolygon unary union.""" + bounds1 = (-180.0, -90.0, 0.0, 90.0) # extent1 = [-180, 0, -90, 90] + bounds2 = (0.0, 0.0, 180.0, 90.0) # extent2 = [0, 180, -90, 90] + polygon1 = Polygon.from_bounds(*bounds1) + polygon2 = Polygon.from_bounds(*bounds2) + multipolygon = MultiPolygon([polygon1, polygon2]) + unified_polygon = unary_union(multipolygon) + assert isinstance(unified_polygon, Polygon) + vertices = np.array(unified_polygon.exterior.coords) + expected_vertices = np.array([[-180., 90.], + [0., 90.], + [180., 90.], + [180., 0.], + [0., 0.], + [0., -90.], + [-180., -90.], + [-180., 90.]]) + np.testing.assert_allclose(expected_vertices, vertices) + + def test_polygon_from_bounds(self): + """Test Polygon definition from bounds.""" + global_polygon = Polygon.from_bounds(-180, -90, 180, 90) + assert global_polygon.bounds == (-180.0, -90.0, 180.0, 90.0) + + def test_polygon_equals(self): + # First polygon goes from -180 to 0 longitude + vertices1 = np.array([[-180., -90.], + [-180., 90.], + [0., 90.], + [0., -90.], + [-180., -90.]]) + polygon1 = Polygon(vertices1) + # Second polygon goes from 0 to 180 longitude + vertices2 = np.array([[0., -90.], + [0., 90.], + [180., 90.], + [180., -90.], + [0., -90.]]) + polygon2 = Polygon(vertices2) + + # Global separate polygon from -180 to 180 + global_separate = MultiPolygon([polygon1, polygon2]) + + # Global single polygon from -180 to 180 + global_polygon = Polygon.from_bounds(-180, -90, 180, 90) + + # Unioned polygon + global_union = unary_union(global_separate) + + # Checks topological equality + assert global_polygon.equals(global_union) + assert global_polygon.equals(global_separate)