diff --git a/.github/workflows/deploy.yaml b/.github/workflows/deploy.yaml index 1bd1f7f4b..cf90d7b87 100644 --- a/.github/workflows/deploy.yaml +++ b/.github/workflows/deploy.yaml @@ -46,10 +46,6 @@ jobs: os: ubuntu-latest python-version: '3.8' docker-image: manylinux2014_x86_64 - - name: manylinux 32-bit - os: ubuntu-latest - python-version: '3.8' - docker-image: manylinux2014_i686 steps: - uses: actions/checkout@v3 @@ -142,4 +138,3 @@ jobs: user: __token__ password: ${{ secrets.pypi_password }} skip_existing: true - diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 9e642420a..de07947f6 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -36,10 +36,9 @@ from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP from pyresample.area_config import create_area_def from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary -from pyresample.utils import ( - check_slice_orientation, +from pyresample.utils import check_slice_orientation, load_cf_area +from pyresample.utils.proj4 import ( get_geostationary_height, - load_cf_area, proj4_dict_to_str, proj4_radius_parameters, ) @@ -1556,7 +1555,7 @@ def _get_geo_boundary_sides(self, frequency=None): # Ensure an even number of vertices for side creation if (frequency % 2) != 0: frequency = frequency + 1 - lons, lats = get_geostationary_bounding_box(self, nb_points=frequency) + lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=frequency) # Retrieve dummy sides for GEO (side1 and side3 always of length 2) side02_step = int(frequency / 2) - 1 lon_sides = [lons[slice(0, side02_step + 1)], @@ -2772,40 +2771,49 @@ def get_geostationary_angle_extent(geos_area): aeq = 1 - req ** 2 / (h ** 2) ap_ = 1 - rp ** 2 / (h ** 2) - # generate points around the north hemisphere in satellite projection - # make it a bit smaller so that we stay inside the valid area - xmax = np.arccos(np.sqrt(aeq)) - ymax = np.arccos(np.sqrt(ap_)) - return xmax, ymax + x_angle = np.arccos(np.sqrt(aeq)) + y_angle = np.arccos(np.sqrt(ap_)) + return x_angle, y_angle def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50): """Get the bbox in geos projection coordinates of the valid pixels inside `geos_area`. - Notes: - - The first and last element of the output vectors are equal. - - If nb_points is even, it will return x and y vectors of length nb_points + 1. - Parameters ---------- nb_points : Number of points on the polygon. """ - xmax, ymax = get_geostationary_angle_extent(geos_area) + x, y = get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points) + ll_x, ll_y, ur_x, ur_y = geos_area.area_extent + + from shapely.geometry import Polygon + geo_bbox = Polygon(np.vstack((x, y)).T) + area_bbox = Polygon(((ll_x, ll_y), (ll_x, ur_y), (ur_x, ur_y), (ur_x, ll_y))) + intersection = area_bbox.intersection(geo_bbox) + try: + x, y = intersection.boundary.xy + except NotImplementedError: + return [], [] + return np.asanyarray(x[:-1]), np.asanyarray(y[:-1]) + + +def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50): + """Get the bbox in geos projection coordinates of the full disk in `geos_area` projection. + + Args: + nb_points: Number of points on the polygon + """ + x_max_angle, y_max_angle = get_geostationary_angle_extent(geos_area) h = get_geostationary_height(geos_area.crs) # generate points around the north hemisphere in satellite projection # make it a bit smaller so that we stay inside the valid area - x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (xmax - 0.0001) - y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (ymax - 0.0001) - - ll_x, ll_y, ur_x, ur_y = geos_area.area_extent - + points_around = np.linspace(-np.pi, np.pi, nb_points, endpoint=False) + x = np.cos(points_around) * (x_max_angle - 0.0001) + y = -np.sin(points_around) * (y_max_angle - 0.0001) x *= h y *= h - # We remove one element with [:-1] to avoid duplicate values at the equator - x = np.clip(np.concatenate([x[:-1], x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x)) - y = np.clip(np.concatenate([y[:-1], -y]), min(ll_y, ur_y), max(ll_y, ur_y)) return x, y diff --git a/pyresample/gradient/__init__.py b/pyresample/gradient/__init__.py index 68d5bf214..57c67c23c 100644 --- a/pyresample/gradient/__init__.py +++ b/pyresample/gradient/__init__.py @@ -39,7 +39,7 @@ from pyresample.geometry import ( AreaDefinition, SwathDefinition, - get_geostationary_bounding_box, + get_geostationary_bounding_box_in_lonlats, ) from pyresample.gradient._gradient_search import ( one_step_gradient_indices, @@ -375,7 +375,7 @@ def _check_input_coordinates(dst_x, dst_y, def get_border_lonlats(geo_def): """Get the border x- and y-coordinates.""" if geo_def.proj_dict['proj'] == 'geos': - lon_b, lat_b = get_geostationary_bounding_box(geo_def, 3600) + lon_b, lat_b = get_geostationary_bounding_box_in_lonlats(geo_def, 3600) else: lons, lats = geo_def.get_boundary_lonlats() lon_b = np.concatenate((lons.side1, lons.side2, lons.side3, lons.side4)) diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index 791cd48d7..883b64ee5 100644 --- a/pyresample/test/test_geometry.py +++ b/pyresample/test/test_geometry.py @@ -2591,8 +2591,100 @@ def _get_fake_antimeridian_lonlats(use_dask: bool) -> tuple: return lons, lats -class TestCrop(unittest.TestCase): - """Test the area helpers.""" +@pytest.fixture +def truncated_geos_area(): + """Create a truncated geostationary area.""" + projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos', + 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'} + area_extent = (5567248.0742, 5570248.4773, -5570248.4773, 1393687.2705) + width = 3712 + height = 1392 + geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent) + return geos_area + + +@pytest.fixture +def truncated_geos_area_in_space(): + """Create a truncated geostationary area.""" + projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos', + 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'} + area_extent = (5575000, 5575000, 5570000, 5570000) + width = 10 + height = 10 + geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent) + return geos_area + + +class TestGeostationaryTools: + """Test the geostationary bbox tools.""" + + def test_get_full_geostationary_bbox(self, truncated_geos_area): + nb_points = 20 + x, y = geometry.get_full_geostationary_bounding_box_in_proj_coords(truncated_geos_area, nb_points) + assert len(x) == nb_points + assert len(y) == nb_points + + assert x[0] != x[-1] + assert y[0] != y[-1] + + expected_x = np.array([-5.43062255e+06, -5.16482897e+06, -4.39346593e+06, -3.19203985e+06, + -1.67815466e+06, 3.32529726e-10, 1.67815466e+06, 3.19203985e+06, + 4.39346593e+06, 5.16482897e+06, 5.43062255e+06, 5.16482897e+06, + 4.39346593e+06, 3.19203985e+06, 1.67815466e+06, 3.32529726e-10, + -1.67815466e+06, -3.19203985e+06, -4.39346593e+06, -5.16482897e+06]) + + expected_y = np.array([6.62789871e-10, 1.67242779e+06, 3.18114670e+06, 4.37847280e+06, + 5.14720348e+06, 5.41209002e+06, 5.14720348e+06, 4.37847280e+06, + 3.18114670e+06, 1.67242779e+06, -0.00000000e+00, -1.67242779e+06, + -3.18114670e+06, -4.37847280e+06, -5.14720348e+06, -5.41209002e+06, + -5.14720348e+06, -4.37847280e+06, -3.18114670e+06, -1.67242779e+06]) + + np.testing.assert_allclose(x, expected_x) + np.testing.assert_allclose(y, expected_y) + + def test_get_geostationary_bbox_works_with_truncated_area(self, truncated_geos_area): + """Ensure the geostationary bbox works when truncated.""" + lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20) + + expected_lon = np.array( + [-64.24072434653284, -68.69662326361153, -65.92516214783112, -60.726360278290336, + -47.39851775032484, 9.500000000000018, 66.39851775032487, 79.72636027829033, + 84.92516214783113, 87.69662326361151, 83.24072434653286]) + expected_lat = np.array( + [14.554922655532085, 17.768795771961937, 35.34328897185421, 52.597860701318254, 69.00533141646078, + 79.1481121862375, 69.00533141646076, 52.597860701318254, 35.34328897185421, 17.768795771961933, + 14.554922655532085]) + np.testing.assert_allclose(lon, expected_lon) + np.testing.assert_allclose(lat, expected_lat) + + def test_get_geostationary_bbox_works_with_truncated_area_proj_coords(self, truncated_geos_area): + """Ensure the geostationary bbox works when truncated.""" + x, y = geometry.get_geostationary_bounding_box_in_proj_coords(truncated_geos_area, 20) + + expected_x = np.array( + [-5209128.302753595, -5164828.965702432, -4393465.934674804, -3192039.8468840676, -1678154.6586309497, + 3.325297262895822e-10, 1678154.6586309501, 3192039.846884068, 4393465.934674805, 5164828.965702432, + 5209128.302753594]) + expected_y = np.array( + [1393687.2705, 1672427.7900638399, 3181146.6955466354, 4378472.798117005, 5147203.47659387, + 5412090.016106332, 5147203.476593869, 4378472.798117005, 3181146.695546635, 1672427.7900638392, + 1393687.2705]) + + np.testing.assert_allclose(x, expected_x) + np.testing.assert_allclose(y, expected_y) + + def test_get_geostationary_bbox_does_not_contain_inf(self, truncated_geos_area): + """Ensure the geostationary bbox does not contain np.inf.""" + lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20) + assert not any(np.isinf(lon)) + assert not any(np.isinf(lat)) + + def test_get_geostationary_bbox_returns_empty_lonlats_in_space(self, truncated_geos_area_in_space): + """Ensure the geostationary bbox is empty when in space.""" + lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area_in_space, 20) + + assert len(lon) == 0 + assert len(lat) == 0 def test_get_geostationary_bbox(self): """Get the geostationary bbox.""" @@ -2606,23 +2698,20 @@ def test_get_geostationary_bbox(self): geos_area.crs = CRS(proj_dict) geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.] - lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20) - # This musk be equal to lon. - elon = np.array([-79.23372832, -78.19662326, -75.42516215, -70.22636028, - -56.89851775, 0., 56.89851775, 70.22636028, - 75.42516215, 78.19662326, 79.23372832, 78.19662326, - 75.42516215, 70.22636028, 56.89851775, 0., - -56.89851775, -70.22636028, -75.42516215, -78.19662326, - -79.23372832]) - elat = np.array([0., 17.76879577, 35.34328897, 52.5978607, - 69.00533142, 79.14811219, 69.00533142, 52.5978607, - 35.34328897, 17.76879577, -0., -17.76879577, - -35.34328897, -52.5978607, -69.00533142, -79.14811219, - -69.00533142, -52.5978607, -35.34328897, -17.76879577, - 0.]) - - np.testing.assert_allclose(lon, elon, atol=1e-07) - np.testing.assert_allclose(lat, elat, atol=1e-07) + lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20) + expected_lon = np.array([-78.19662326, -75.42516215, -70.22636028, + -56.89851775, 0., 56.89851775, 70.22636028, + 75.42516215, 78.19662326, 79.23372832, 78.19662326, + 75.42516215, 70.22636028, 56.89851775, 0., + -56.89851775, -70.22636028, -75.42516215, -78.19662326, -79.23372832, ]) + expected_lat = np.array([17.76879577, 35.34328897, 52.5978607, + 69.00533142, 79.14811219, 69.00533142, 52.5978607, + 35.34328897, 17.76879577, -0., -17.76879577, + -35.34328897, -52.5978607, -69.00533142, -79.14811219, + -69.00533142, -52.5978607, -35.34328897, -17.76879577, 0.]) + + np.testing.assert_allclose(lon, expected_lon, atol=1e-07) + np.testing.assert_allclose(lat, expected_lat, atol=1e-07) geos_area = MagicMock() lon_0 = 10 @@ -2634,8 +2723,8 @@ def test_get_geostationary_bbox(self): geos_area.crs = CRS(proj_dict) geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.] - lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20) - np.testing.assert_allclose(lon, elon + lon_0) + lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20) + np.testing.assert_allclose(lon, expected_lon + lon_0) def test_get_geostationary_angle_extent(self): """Get max geostationary angles.""" @@ -2675,6 +2764,10 @@ def test_get_geostationary_angle_extent(self): np.testing.assert_allclose(expected, geometry.get_geostationary_angle_extent(geos_area)) + +class TestCrop(unittest.TestCase): + """Test the area helpers.""" + def test_sub_area(self): """Sub area slicing.""" area = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', @@ -3156,20 +3249,19 @@ def test_geostationary_projection(self): # Check boundary vertices n_vertices = 10 boundary = areadef.boundary(frequency=n_vertices, force_clockwise=False) - boundary.vertices.shape # Check boundary vertices is in correct order - expected_vertices = np.array([[-7.92337283e+01, 6.94302533e-15], - [-7.54251621e+01, 3.53432890e+01], + expected_vertices = np.array([[-7.54251621e+01, 3.53432890e+01], [-5.68985178e+01, 6.90053314e+01], [5.68985178e+01, 6.90053314e+01], [7.54251621e+01, 3.53432890e+01], - [7.92337283e+01, -6.94302533e-15], + [7.92337283e+01, -0.00000000e+00], [7.54251621e+01, -3.53432890e+01], [5.68985178e+01, -6.90053314e+01], [-5.68985178e+01, -6.90053314e+01], - [-7.54251621e+01, -3.53432890e+01]]) - assert np.allclose(expected_vertices, boundary.vertices) + [-7.54251621e+01, -3.53432890e+01], + [-7.92337283e+01, 6.94302533e-15]]) + np.testing.assert_allclose(expected_vertices, boundary.vertices) def test_global_platee_caree_projection(self): """Test boundary for global platee caree projection.""" diff --git a/pyresample/test/test_gradient.py b/pyresample/test/test_gradient.py index c28a190c5..59d9b0aad 100644 --- a/pyresample/test/test_gradient.py +++ b/pyresample/test/test_gradient.py @@ -599,7 +599,7 @@ def test_check_overlap(): assert check_overlap(poly1, poly2) is False -@mock.patch('pyresample.gradient.get_geostationary_bounding_box') +@mock.patch('pyresample.gradient.get_geostationary_bounding_box_in_lonlats') def test_get_border_lonlats(get_geostationary_bounding_box): """Test that correct methods are called in get_border_lonlats().""" from pyresample.gradient import get_border_lonlats diff --git a/setup.py b/setup.py index 4dd9ba43e..7a0bae933 100644 --- a/setup.py +++ b/setup.py @@ -29,13 +29,14 @@ requirements = ['setuptools>=3.2', 'pyproj>=3.0', 'configobj', 'pykdtree>=1.3.1', 'pyyaml', 'numpy>=1.10.0', + "shapely", ] if sys.version_info < (3, 10): requirements.append('importlib_metadata') -test_requires = ['rasterio', 'dask', 'xarray', 'cartopy', 'pillow', 'matplotlib', 'scipy', 'zarr', - 'pytest-lazy-fixtures'] +test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr', + 'pytest-lazy-fixtures', 'shapely'] extras_require = {'numexpr': ['numexpr'], 'quicklook': ['matplotlib', 'cartopy>=0.20.0', 'pillow'], 'rasterio': ['rasterio'], @@ -46,8 +47,6 @@ 'tests': test_requires} setup_requires = ['numpy>=1.10.0', 'cython'] -test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr', - 'pytest-lazy-fixture'] if sys.platform.startswith("win"): extra_compile_args = []