diff --git a/CHANGELOG.md b/CHANGELOG.md index 03075cd..f233352 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,12 @@ -## [unreleased] - 2025-10-20 +## [v1.1.13] - 2025-10-21 ### Changed -- Updated tests to be less dependent on architecture when comparing floats. -- Fixed test that modified a source file fixture. -- Changed infrastructure so that local and Docker runs of the tests produce output in same locations. +- Updates evaluation of bbox or polygon constraint to exclude areas outside the + projected target grid. An error exception occurs if spatial constraint is entirely +- Updates tests to be less dependent on architecture when comparing floats. +- Fixes test that modified a source file fixture. +- Changes infrastructure so that local and Docker runs of the tests produce output in same locations. - GitHub once again captures the artifacts from the tests and coverage. @@ -172,6 +174,7 @@ Repository structure changes include: For more information on internal releases prior to NASA open-source approval, see legacy-CHANGELOG.md. +[v1.1.13]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.1.13 [v1.1.12]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.1.12 [v1.1.11]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.1.11 [v1.1.10]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.1.10 diff --git a/docker/service_version.txt b/docker/service_version.txt index ccad953..9ea63db 100644 --- a/docker/service_version.txt +++ b/docker/service_version.txt @@ -1 +1 @@ -1.1.12 +1.1.13 diff --git a/hoss/exceptions.py b/hoss/exceptions.py index 0e01611..403f788 100644 --- a/hoss/exceptions.py +++ b/hoss/exceptions.py @@ -61,6 +61,19 @@ def __init__(self): ) +class InvalidGranuleDimensions(CustomError): + """This exception is raised when the granule dimensions are not valid for + the specific projection of the granule. + + """ + + def __init__(self): + super().__init__( + 'InvalidGranuleDimensions', + 'The dimensions used for the granule appear to be invalid for the crs', + ) + + class MissingGridMappingMetadata(CustomError): """This exception is raised when HOSS tries to obtain the `grid_mapping` metadata attribute for a projected variable and it is not present in diff --git a/hoss/projection_utilities.py b/hoss/projection_utilities.py index 3c55f32..db28e8d 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -29,7 +29,9 @@ from hoss.bbox_utilities import BBox, flatten_list from hoss.exceptions import ( + InvalidGranuleDimensions, InvalidInputGeoJSON, + InvalidRequestedRange, MissingGridMappingMetadata, MissingGridMappingVariable, MissingSpatialSubsetInformation, @@ -203,12 +205,16 @@ def get_projected_x_y_extents( grid_lats, grid_lons = get_grid_lat_lons( # pylint: disable=unpacking-non-sequence x_values, y_values, crs ) + if not np.all(np.isfinite(grid_lats)) or not np.all(np.isfinite(grid_lons)): + raise InvalidGranuleDimensions + # When projected, the perimeter of a bounding box or polygon in geographic # terms will become curved. To determine the X, Y extent of the requested # bounding area, we need a perimeter with a suitable density of points, such that # we catch that curve when projected. The source file resolution is used to define # the necessary number of points (density). geographic_resolution = get_geographic_resolution(grid_lons, grid_lats) + densified_perimeter = get_densified_perimeter( geographic_resolution, shape_file=shape_file, bounding_box=bounding_box ) @@ -221,13 +227,22 @@ def get_projected_x_y_extents( clipped_perimeter = get_filtered_points(densified_perimeter, granule_extent) - return get_x_y_extents_from_geographic_points(clipped_perimeter, crs) + granule_extent_projected_meters = { + "x_min": np.min(x_values), + "x_max": np.max(x_values), + "y_min": np.min(y_values), + "y_max": np.max(y_values), + } + return get_x_y_extents_from_geographic_perimeter( + clipped_perimeter, crs, granule_extent_projected_meters + ) def get_filtered_points( points_in_requested_extent: List[Coordinates], granule_extent: BBox ) -> List[Coordinates]: - """Returns lat/lon values cropped to the extent of the granule""" + """Returns lat/lon values clipped to the extent of the granule""" + requested_lons, requested_lats = zip(*points_in_requested_extent) # all lon values are clipped within the granule lon extent @@ -458,29 +473,99 @@ def get_resolved_line( return list(zip(new_x, new_y)) -def get_x_y_extents_from_geographic_points( - points: List[Coordinates], crs: CRS +def get_x_y_extents_from_geographic_perimeter( + points: List[Coordinates], + crs: CRS, + granule_extent_projected_meters: dict[str, float], ) -> Dict[str, float]: """Take an input list of (longitude, latitude) coordinates that define the exterior of the input GeoJSON shape or bounding box, and project those points to the target grid. Then return the minimum and maximum values - of those projected coordinates. + of those projected coordinates. Check first for perimeter exceeding grid on + all axes (whole grid extents returned). Then remove any points that are + outside the grid before finding the min and max extent. """ + # get the x,y projected values from the geographic points point_longitudes, point_latitudes = zip(*points) from_geo_transformer = Transformer.from_crs(4326, crs) points_x, points_y = ( # pylint: disable=unpacking-non-sequence from_geo_transformer.transform(point_latitudes, point_longitudes) ) - # isfinite checks for NaN and infinty values returned for certain projections - points_x = np.asarray(points_x) - points_y = np.asarray(points_y) - finite_x = points_x[np.isfinite(points_x)] - finite_y = points_y[np.isfinite(points_y)] + finite_x, finite_y = remove_non_finite_projected_values(points_x, points_y) + + # Check if perimeter exceeds the grid extents on all axes. If true, return + # whole grid extents and skips the code that follows (which fails in + # this case). + if perimeter_surrounds_grid(finite_x, finite_y, granule_extent_projected_meters): + return granule_extent_projected_meters + + # Remove any points that are outside the grid + finite_x, finite_y = remove_points_outside_grid_extents( + finite_x, finite_y, granule_extent_projected_meters + ) + return { 'x_min': np.min(finite_x), 'x_max': np.max(finite_x), 'y_min': np.min(finite_y), 'y_max': np.max(finite_y), } + + +def remove_non_finite_projected_values( + points_x: list[float], points_y: list[float] +) -> tuple[np.ndarray, np.ndarray]: + """Removes any NaN and infinity values and returns the results as numpy arrays""" + # isfinite checks for NaN and infinty values returned for certain projections + points_x = np.asarray(points_x) + points_y = np.asarray(points_y) + + finite_mask = np.isfinite(points_x) & np.isfinite(points_y) + finite_x = points_x[finite_mask] + finite_y = points_y[finite_mask] + return finite_x, finite_y + + +def perimeter_surrounds_grid( + finite_x: np.ndarray, finite_y: np.ndarray, granule_extent: dict[str, float] +) -> bool: + """Returns True if perimeter exceeds the grid extents on all axes. + Returns False if does not. + + """ + + if ( + np.min(finite_x) < granule_extent["x_min"] + and np.max(finite_x) > granule_extent["x_max"] + and np.min(finite_y) < granule_extent["y_min"] + and np.max(finite_y) > granule_extent["y_max"] + ): + return True + + return False + + +def remove_points_outside_grid_extents( + finite_x: np.ndarray, finite_y: np.ndarray, granule_extent: dict[str, float] +) -> tuple[np.ndarray, np.ndarray]: + """Remove any points that are outside the grid and are invalid and raise an + exception if the resulting grid is empty. + + """ + + mask = ( + (finite_x >= granule_extent["x_min"]) + & (finite_x <= granule_extent["x_max"]) + & (finite_y >= granule_extent["y_min"]) + & (finite_y <= granule_extent["y_max"]) + ) + + finite_x = finite_x[mask] + finite_y = finite_y[mask] + + if finite_x.size == 0 or finite_y.size == 0: + raise InvalidRequestedRange + + return finite_x, finite_y diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 2695b9d..3758049 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -20,7 +20,9 @@ from hoss.bbox_utilities import BBox from hoss.exceptions import ( + InvalidGranuleDimensions, InvalidInputGeoJSON, + InvalidRequestedRange, MissingGridMappingMetadata, MissingGridMappingVariable, MissingSpatialSubsetInformation, @@ -40,9 +42,12 @@ get_resolved_geometry, get_resolved_line, get_variable_crs, - get_x_y_extents_from_geographic_points, + get_x_y_extents_from_geographic_perimeter, is_projection_x_dimension, is_projection_y_dimension, + perimeter_surrounds_grid, + remove_non_finite_projected_values, + remove_points_outside_grid_extents, ) from tests import utilities from tests.utilities import assert_float_dict_almost_equal @@ -400,10 +405,10 @@ def test_get_projected_x_y_extents_whole_earth(self): } ) expected_output = { - 'x_min': -12702459.818865139, - 'x_max': 12702459.818865139, - 'y_min': -12702440.623773243, - 'y_max': 12702440.710450241, + 'x_min': -8982000, + 'x_max': 8982000, + 'y_min': -8982000, + 'y_max': 8982000, } with self.subTest('Whole Earth LAEA - Bounding box input'): assert_float_dict_almost_equal( @@ -421,6 +426,52 @@ def test_get_projected_x_y_extents_whole_earth(self): expected_output, ) + def test_get_projected_x_y_extents_edge_case(self): + """Ensure that the expected values for the x and y dimension extents + are recovered for a polar projected grid and when a bounding box + for any edge cases are requested. + + """ + bbox = BBox(-89, -79, -40, -59) + + x_values = np.linspace(-9000000, 9000000, 2000) + y_values = np.linspace(9000000, -9000000, 2000) + + crs = CRS.from_cf( + { + 'false_easting': 0.0, + 'false_northing': 0.0, + 'longitude_of_central_meridian': 0.0, + 'latitude_of_projection_origin': 90.0, + 'grid_mapping_name': 'lambert_azimuthal_equal_area', + } + ) + + expected_output = { + 'x_min': -8993061.78423412, + 'x_max': -8350580.505440015, + 'y_min': -8997181.591145469, + 'y_max': -8354987.361637551, + } + with self.subTest( + 'LAEA - Bounding box which is close to the edge of granule extent' + ): + assert_float_dict_almost_equal( + get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox), + expected_output, + ) + + bbox = BBox(-90, -87, -75, -85) + with self.subTest('LAEA - Bounding box which is outside the granule extent'): + with self.assertRaises(InvalidRequestedRange): + get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox) + + with self.subTest('When the granule has invalid dimensions'): + with self.assertRaises(InvalidGranuleDimensions): + x_values1 = np.linspace(-9200000, 9200000, 500) + y_values1 = np.linspace(9200000, -9200000, 500) + get_projected_x_y_extents(x_values1, y_values1, crs, bounding_box=bbox) + def test_get_filtered_points(self): """Ensure that the coordinates returned are clipped to the granule extent or the bbox extent whichever is the smaller of the two. @@ -1148,12 +1199,19 @@ def test_get_resolved_line(self): get_resolved_line(point_one, point_two, resolution), expected_output ) - def test_get_x_y_extents_from_geographic_points(self): + def test_get_x_y_extents_from_geographic_perimeter(self): """Ensure that a list of coordinates is transformed to a specified projection, and that the expected extents in the projected x and y dimensions are returned. """ + + granule_extent = { + "x_min": -9000000, + "x_max": 9000000, + "y_min": -9000000, + "y_max": 9000000, + } points = [(-180, 75), (-90, 75), (0, 75), (90, 75)] crs = CRS.from_epsg(6931) expected_x_y_extents = { @@ -1164,15 +1222,22 @@ def test_get_x_y_extents_from_geographic_points(self): } assert_float_dict_almost_equal( - get_x_y_extents_from_geographic_points(points, crs), expected_x_y_extents + get_x_y_extents_from_geographic_perimeter(points, crs, granule_extent), + expected_x_y_extents, ) - def test_get_x_y_extents_from_geographic_points_full_earth_laea(self): + def test_get_x_y_extents_from_geographic_perimeter_full_earth_laea(self): """Ensure that a list of coordinates is transformed to the specified laea projection, and valid values in the projected x and y dimensions are returned even for edge cases like whole earth. """ + granule_extent_projected_meters = { + "x_min": -9000000, + "x_max": 9000000, + "y_min": -9000000, + "y_max": 9000000, + } crs = CRS.from_cf( { 'false_easting': 0.0, @@ -1184,10 +1249,144 @@ def test_get_x_y_extents_from_geographic_points_full_earth_laea(self): ) points1 = [(-180, -90), (-180, 90), (180, 90), (180, -90)] - x_y_extents = get_x_y_extents_from_geographic_points(points1, crs) + x_y_extents = get_x_y_extents_from_geographic_perimeter( + points1, crs, granule_extent_projected_meters + ) self.assertTrue(all(not math.isinf(value) for value in x_y_extents.values())) + def test_remove_non_finite_projected_values(self): + """Ensure that only valid values in the x,y list are returned and any + NaN or inf values are removed. + """ + points_x = [ + float('-inf'), + -900.3, + -800.2, + -700.1, + float('nan'), + 500.0, + 600.9, + 700.0, + 800.0, + ] + points_y = [ + 89.0, + float('nan'), + 69.1, + 40.5, + -10.6, + 50.9, + 70.2, + 80.4, + float('inf'), + ] + expected_x_values = np.array([-800.2, -700.1, 500.0, 600.9, 700.0]) + expected_y_values = np.array([69.1, 40.5, 50.9, 70.2, 80.4]) + valid_x, valid_y = remove_non_finite_projected_values(points_x, points_y) + np.array_equal(valid_x, expected_x_values) + np.array_equal(valid_y, expected_y_values) + + def test_check_perimeter_exceeds_grid_extents(self): + """Ensure that True value is returned when the input perimeter array exceeds grid + extents and a False value is returned when the perimeter is within the grid extents + + """ + granule_extent = { + "x_min": -1000.0, + "x_max": 1000.0, + "y_min": -1000.0, + "y_max": 1000.0, + } + with self.subTest("Perimeter exceeds grid extents in all axes"): + self.assertTrue( + perimeter_surrounds_grid( + np.array( + [-2000.1, -1000.1, -900, -500, 200.3, 700.1, 800.1, 1200.5] + ), + np.array( + [1100.1, 400.9, 200, -100.8, -300.3, -500.1, -600.1, -1000.5] + ), + granule_extent, + ) + ) + with self.subTest("Perimeter does not surround grid"): + self.assertFalse( + perimeter_surrounds_grid( + np.array( + [-2000.1, -1000.1, -900, -500, 200.3, 700.1, 800.1, 900.0] + ), + np.array( + [400.1, 300.9, 200, -100.8, -300.3, -500.1, -600.1, -700.5] + ), + granule_extent, + ) + ) + + def test_remove_points_outside_grid_extents(self): + """Ensure that any point outside the grid extents and removed and the grid returned + only has values within the grid extent. + + """ + points_x = np.array( + [ + -1100.1, + -1000.0, + -800.3, + -700.4, + -700.5, + -700.6, + -700.7, + -800.8, + -900.9, + -1000.0, + -1100.1, + ] + ) + points_y = np.array( + [ + 600.1, + 700.2, + 800.3, + 900.4, + 1000.0, + 1100.6, + 1000.0, + 900.8, + 800.9, + 700.0, + 600.1, + ] + ) + expected_x = np.array( + [-1000.0, -800.3, -700.4, -700.5, -700.7, -800.8, -900.9, -1000.0] + ) + expected_y = np.array( + [700.2, 800.3, 900.4, 1000.0, 1000.0, 900.8, 800.9, 700.0] + ) + + granule_extent = { + "x_min": -1000.0, + "x_max": 1000.0, + "y_min": -1000.0, + "y_max": 1000.0, + } + + with self.subTest("Perimeter contains some valid points within grid extent"): + x, y = remove_points_outside_grid_extents( + points_x, points_y, granule_extent + ) + np.array_equal(expected_x, x) + np.array_equal(expected_y, y) + + with self.subTest("Perimeter has no valid points within the grid extent"): + with self.assertRaises(InvalidRequestedRange): + remove_points_outside_grid_extents( + np.array([-2000.1, -1000.1, -900, -500, 200.3, 700.1, 1200.5]), + np.array([600, 800.9, 1000.2, 1100.8, 1100.1, 1000.1, 900.5]), + granule_extent, + ) + @patch('hoss.projection_utilities.get_grid_mapping_attributes') def test_get_master_geotransform(self, mock_get_grid_mapping_attributes): """Ensure that the `master_geotransform` attribute is returned. If it doesn't