From c9f3cc4c1bd1b1c63c790ea5bceb9dc1be379339 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Fri, 17 Oct 2025 14:07:06 -0400 Subject: [PATCH 01/11] DAS-2424 Update get_x_y_extents_from_geographic_points to remove points outside granule extent --- CHANGELOG.md | 9 ++++ docker/service_version.txt | 2 +- hoss/projection_utilities.py | 48 +++++++++++++++++++-- tests/unit/test_projection_utilities.py | 57 ++++++++++++++++++++++--- 4 files changed, 105 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8501142..d55e490 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,11 @@ +## [v1.1.13] - 2025-10-16 + +### Changed + +- Updates determination of projected min-max in the bounding box calculations to correct + edge cases where the area outside the granules gets included for lambert-azimuthal + grid projections. An exception is thrown for invalid spatial extent ranges. + ## [v1.1.12] - 2025-10-15 ### Added @@ -162,6 +170,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/projection_utilities.py b/hoss/projection_utilities.py index 3c55f32..da2860b 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -30,6 +30,7 @@ from hoss.bbox_utilities import BBox, flatten_list from hoss.exceptions import ( InvalidInputGeoJSON, + InvalidRequestedRange, MissingGridMappingMetadata, MissingGridMappingVariable, MissingSpatialSubsetInformation, @@ -203,12 +204,14 @@ def get_projected_x_y_extents( grid_lats, grid_lons = get_grid_lat_lons( # pylint: disable=unpacking-non-sequence x_values, y_values, crs ) + # 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 +224,16 @@ 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) + return get_x_y_extents_from_geographic_points( + clipped_perimeter, crs, x_values, y_values + ) 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 @@ -459,12 +465,13 @@ def get_resolved_line( def get_x_y_extents_from_geographic_points( - points: List[Coordinates], crs: CRS + points: List[Coordinates], crs: CRS, x_values: np.ndarray, y_values: np.ndarray ) -> 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. Remove any points that are outside the + granule extent. """ point_longitudes, point_latitudes = zip(*points) @@ -478,6 +485,39 @@ def get_x_y_extents_from_geographic_points( finite_x = points_x[np.isfinite(points_x)] finite_y = points_y[np.isfinite(points_y)] + + if ( + np.min(finite_x) < np.min(x_values) + and np.max(finite_x) > np.max(x_values) + and np.min(finite_y) < np.min(y_values) + and np.max(finite_y) > np.max(y_values) + ): + + return { + 'x_min': np.min(x_values), + 'x_max': np.max(x_values), + 'y_min': np.min(y_values), + 'y_max': np.max(y_values), + } + + delete_x_min = np.where(finite_x < np.min(x_values)) + finite_x = np.delete(finite_x, delete_x_min) + finite_y = np.delete(finite_y, delete_x_min) + + delete_x_max = np.where(finite_x > np.max(x_values)) + finite_x = np.delete(finite_x, delete_x_max) + finite_y = np.delete(finite_y, delete_x_max) + + delete_y_min = np.where(finite_y < np.min(y_values)) + finite_y = np.delete(finite_y, delete_y_min) + finite_x = np.delete(finite_x, delete_y_min) + + delete_y_max = np.where(finite_y > np.max(y_values)) + finite_y = np.delete(finite_y, delete_y_max) + finite_x = np.delete(finite_x, delete_y_max) + if finite_x.size == 0 or finite_y.size == 0: + raise InvalidRequestedRange + return { 'x_min': np.min(finite_x), 'x_max': np.max(finite_x), diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 54652c1..6106cd5 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -21,6 +21,7 @@ from hoss.bbox_utilities import BBox from hoss.exceptions import ( InvalidInputGeoJSON, + InvalidRequestedRange, MissingGridMappingMetadata, MissingGridMappingVariable, MissingSpatialSubsetInformation, @@ -397,10 +398,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'): self.assertDictEqual( @@ -418,6 +419,43 @@ 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 an edge case is 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 input - edge case 1'): + self.assertDictEqual( + get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox), + expected_output, + ) + + bbox = BBox(-89, -79, -40, -59) + with self.subTest('LAEA - Bounding box input - edge case 2'): + self.assertRaises(InvalidRequestedRange) + 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. @@ -1151,6 +1189,8 @@ def test_get_x_y_extents_from_geographic_points(self): dimensions are returned. """ + x_values = np.linspace(-9000000, 9000000, 2000) + y_values = np.linspace(9000000, -9000000, 2000) points = [(-180, 75), (-90, 75), (0, 75), (90, 75)] crs = CRS.from_epsg(6931) expected_x_y_extents = { @@ -1161,7 +1201,8 @@ def test_get_x_y_extents_from_geographic_points(self): } self.assertDictEqual( - get_x_y_extents_from_geographic_points(points, crs), expected_x_y_extents + get_x_y_extents_from_geographic_points(points, crs, x_values, y_values), + expected_x_y_extents, ) def test_get_x_y_extents_from_geographic_points_full_earth_laea(self): @@ -1170,6 +1211,8 @@ def test_get_x_y_extents_from_geographic_points_full_earth_laea(self): dimensions are returned even for edge cases like whole earth. """ + x_values = np.linspace(-9000000, 9000000, 2000) + y_values = np.linspace(9000000, -9000000, 2000) crs = CRS.from_cf( { 'false_easting': 0.0, @@ -1181,7 +1224,9 @@ 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_points( + points1, crs, x_values, y_values + ) self.assertTrue(all(not math.isinf(value) for value in x_y_extents.values())) From e341b9b733c28b8eba8635b08f943fda1651aa1c Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Fri, 17 Oct 2025 18:02:37 -0400 Subject: [PATCH 02/11] DAS-2424 Add comments and change a function name --- hoss/projection_utilities.py | 17 ++++++++++++----- tests/unit/test_projection_utilities.py | 10 +++++----- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/hoss/projection_utilities.py b/hoss/projection_utilities.py index da2860b..41c8f79 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -224,7 +224,7 @@ def get_projected_x_y_extents( clipped_perimeter = get_filtered_points(densified_perimeter, granule_extent) - return get_x_y_extents_from_geographic_points( + return get_x_y_extents_from_geographic_perimeter( clipped_perimeter, crs, x_values, y_values ) @@ -464,14 +464,15 @@ def get_resolved_line( return list(zip(new_x, new_y)) -def get_x_y_extents_from_geographic_points( +def get_x_y_extents_from_geographic_perimeter( points: List[Coordinates], crs: CRS, x_values: np.ndarray, y_values: np.ndarray ) -> 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. Remove any points that are outside the - granule extent. + of those projected coordinates. Check first for permiter 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. """ point_longitudes, point_latitudes = zip(*points) @@ -486,6 +487,11 @@ def get_x_y_extents_from_geographic_points( finite_x = points_x[np.isfinite(points_x)] finite_y = points_y[np.isfinite(points_y)] + # Check if perimeter exceeds the grid extents on all axes. If true, return + # whole grid extents. This handles the case where the perimeter wholly + # encloses the grid (e.g., whole-earth bbox, any lesser grid, polar + # use-cases in particular), and skips the code that follows (which fails in + # this case). if ( np.min(finite_x) < np.min(x_values) and np.max(finite_x) > np.max(x_values) @@ -499,7 +505,8 @@ def get_x_y_extents_from_geographic_points( 'y_min': np.min(y_values), 'y_max': np.max(y_values), } - + # Remove any points that are outside the grid and are invalid before + # determining min/max extents of perimeter inside the grid. delete_x_min = np.where(finite_x < np.min(x_values)) finite_x = np.delete(finite_x, delete_x_min) finite_y = np.delete(finite_y, delete_x_min) diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 6106cd5..1aefb34 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -41,7 +41,7 @@ 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, ) @@ -1183,7 +1183,7 @@ 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. @@ -1201,11 +1201,11 @@ def test_get_x_y_extents_from_geographic_points(self): } self.assertDictEqual( - get_x_y_extents_from_geographic_points(points, crs, x_values, y_values), + get_x_y_extents_from_geographic_perimeter(points, crs, x_values, y_values), 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. @@ -1224,7 +1224,7 @@ 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( + x_y_extents = get_x_y_extents_from_geographic_perimeter( points1, crs, x_values, y_values ) From 3b642b7c1711a9b48f4b4dc7a8c4d07659488666 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Mon, 20 Oct 2025 12:34:46 -0400 Subject: [PATCH 03/11] DAS-2424 Correct incomplete unit unit and fix comments --- CHANGELOG.md | 8 ++++---- hoss/projection_utilities.py | 2 +- tests/unit/test_projection_utilities.py | 5 +++-- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d55e490..1bcb91d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,10 @@ -## [v1.1.13] - 2025-10-16 +## [v1.1.13] - 2025-10-20 ### Changed -- Updates determination of projected min-max in the bounding box calculations to correct - edge cases where the area outside the granules gets included for lambert-azimuthal - grid projections. An exception is thrown for invalid spatial extent ranges. +- Updates evaluation of bbox or polygon constraint to exclude areas outside the + projected target grid. An exception occurs if spatial constraint is entirely + outside grid. ## [v1.1.12] - 2025-10-15 diff --git a/hoss/projection_utilities.py b/hoss/projection_utilities.py index 41c8f79..5925e4a 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -470,7 +470,7 @@ def get_x_y_extents_from_geographic_perimeter( """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. Check first for permiter exceeding grid on + 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. diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 1aefb34..f76e4c5 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -452,9 +452,10 @@ def test_get_projected_x_y_extents_edge_case(self): expected_output, ) - bbox = BBox(-89, -79, -40, -59) + bbox = BBox(-90, -87, -75, -85) with self.subTest('LAEA - Bounding box input - edge case 2'): - self.assertRaises(InvalidRequestedRange) + with self.assertRaises(InvalidRequestedRange): + get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox) def test_get_filtered_points(self): """Ensure that the coordinates returned are clipped to the granule extent or From 7cab64eb3a461838345e5378e54ce5e7798c4fed Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Mon, 20 Oct 2025 12:38:46 -0400 Subject: [PATCH 04/11] DAS-2424 Correct CHANGELOG comments --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1bcb91d..c78f51b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,7 +3,7 @@ ### Changed - Updates evaluation of bbox or polygon constraint to exclude areas outside the - projected target grid. An exception occurs if spatial constraint is entirely + projected target grid. An error exception occurs if spatial constraint is entirely outside grid. ## [v1.1.12] - 2025-10-15 From 222e79b9d84b32a1e281774e1ae786075a696427 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Mon, 20 Oct 2025 12:46:27 -0400 Subject: [PATCH 05/11] DAS-2424 Correct test case title --- tests/unit/test_projection_utilities.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index f76e4c5..27c9391 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -422,7 +422,7 @@ def test_get_projected_x_y_extents_whole_earth(self): 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 an edge case is requested. + for any edge cases are requested. """ bbox = BBox(-89, -79, -40, -59) @@ -446,14 +446,16 @@ def test_get_projected_x_y_extents_edge_case(self): 'y_min': -8997181.591145469, 'y_max': -8354987.361637551, } - with self.subTest('LAEA - Bounding box input - edge case 1'): + with self.subTest( + 'LAEA - Bounding box which is close to the edge of granule extent' + ): self.assertDictEqual( 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 input - edge case 2'): + 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) From 93198a687c0986f23175c20d45a7dc033d53d020 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Tue, 21 Oct 2025 03:12:52 -0400 Subject: [PATCH 06/11] DAS-2424 - modularize code into smaller functions and add unit tests --- hoss/projection_utilities.py | 122 ++++++++++++++++-------- tests/unit/test_projection_utilities.py | 103 ++++++++++++++++++-- 2 files changed, 180 insertions(+), 45 deletions(-) diff --git a/hoss/projection_utilities.py b/hoss/projection_utilities.py index 5925e4a..1bad149 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -224,8 +224,14 @@ def get_projected_x_y_extents( clipped_perimeter = get_filtered_points(densified_perimeter, granule_extent) + 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, x_values, y_values + clipped_perimeter, crs, granule_extent_projected_meters ) @@ -465,7 +471,9 @@ def get_resolved_line( def get_x_y_extents_from_geographic_perimeter( - points: List[Coordinates], crs: CRS, x_values: np.ndarray, y_values: np.ndarray + 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 @@ -475,55 +483,32 @@ def get_x_y_extents_from_geographic_perimeter( 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_any_invalid_projected_values(points_x, points_y) # Check if perimeter exceeds the grid extents on all axes. If true, return - # whole grid extents. This handles the case where the perimeter wholly - # encloses the grid (e.g., whole-earth bbox, any lesser grid, polar - # use-cases in particular), and skips the code that follows (which fails in + # whole grid extents and skips the code that follows (which fails in # this case). - if ( - np.min(finite_x) < np.min(x_values) - and np.max(finite_x) > np.max(x_values) - and np.min(finite_y) < np.min(y_values) - and np.max(finite_y) > np.max(y_values) + if check_perimeter_exceeds_grid_extents( + finite_x, finite_y, granule_extent_projected_meters ): - return { - 'x_min': np.min(x_values), - 'x_max': np.max(x_values), - 'y_min': np.min(y_values), - 'y_max': np.max(y_values), + 'x_min': granule_extent_projected_meters["x_min"], + 'x_max': granule_extent_projected_meters["x_max"], + 'y_min': granule_extent_projected_meters["y_min"], + 'y_max': granule_extent_projected_meters["y_max"], } - # Remove any points that are outside the grid and are invalid before - # determining min/max extents of perimeter inside the grid. - delete_x_min = np.where(finite_x < np.min(x_values)) - finite_x = np.delete(finite_x, delete_x_min) - finite_y = np.delete(finite_y, delete_x_min) - - delete_x_max = np.where(finite_x > np.max(x_values)) - finite_x = np.delete(finite_x, delete_x_max) - finite_y = np.delete(finite_y, delete_x_max) - - delete_y_min = np.where(finite_y < np.min(y_values)) - finite_y = np.delete(finite_y, delete_y_min) - finite_x = np.delete(finite_x, delete_y_min) - - delete_y_max = np.where(finite_y > np.max(y_values)) - finite_y = np.delete(finite_y, delete_y_max) - finite_x = np.delete(finite_x, delete_y_max) - if finite_x.size == 0 or finite_y.size == 0: - raise InvalidRequestedRange + + # 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), @@ -531,3 +516,62 @@ def get_x_y_extents_from_geographic_perimeter( 'y_min': np.min(finite_y), 'y_max': np.max(finite_y), } + + +def remove_any_invalid_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 check_perimeter_exceeds_grid_extents( + finite_x: np.ndarray, finite_y: np.ndarray, granule_extent: dict[str, float] +) -> bool: + """Check if perimeter exceeds the grid extents on all axes. If true, return + whole grid extents. This handles the case where the perimeter wholly + encloses the grid (e.g., whole-earth bbox, any lesser grid, polar + use-cases in particular) + + """ + + 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 27c9391..6c88619 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -27,6 +27,7 @@ MissingSpatialSubsetInformation, ) from hoss.projection_utilities import ( + check_perimeter_exceeds_grid_extents, get_bbox_polygon, get_densified_perimeter, get_filtered_points, @@ -44,6 +45,8 @@ get_x_y_extents_from_geographic_perimeter, is_projection_x_dimension, is_projection_y_dimension, + remove_any_invalid_projected_values, + remove_points_outside_grid_extents, ) @@ -1192,8 +1195,13 @@ def test_get_x_y_extents_from_geographic_perimeter(self): dimensions are returned. """ - x_values = np.linspace(-9000000, 9000000, 2000) - y_values = np.linspace(9000000, -9000000, 2000) + + 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 = { @@ -1204,7 +1212,7 @@ def test_get_x_y_extents_from_geographic_perimeter(self): } self.assertDictEqual( - get_x_y_extents_from_geographic_perimeter(points, crs, x_values, y_values), + get_x_y_extents_from_geographic_perimeter(points, crs, granule_extent), expected_x_y_extents, ) @@ -1214,8 +1222,12 @@ def test_get_x_y_extents_from_geographic_perimeter_full_earth_laea(self): dimensions are returned even for edge cases like whole earth. """ - x_values = np.linspace(-9000000, 9000000, 2000) - y_values = np.linspace(9000000, -9000000, 2000) + granule_extent_projected_meters = { + "x_min": -9000000, + "x_max": 9000000, + "y_min": -9000000, + "y_max": 9000000, + } crs = CRS.from_cf( { 'false_easting': 0.0, @@ -1228,11 +1240,90 @@ def test_get_x_y_extents_from_geographic_perimeter_full_earth_laea(self): points1 = [(-180, -90), (-180, 90), (180, 90), (180, -90)] x_y_extents = get_x_y_extents_from_geographic_perimeter( - points1, crs, x_values, y_values + points1, crs, granule_extent_projected_meters ) self.assertTrue(all(not math.isinf(value) for value in x_y_extents.values())) + def test_remove_any_invalid_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_any_invalid_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, "x_max": 1000, "y_min": -1000, "y_max": 1000} + self.assertTrue( + check_perimeter_exceeds_grid_extents( + 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, + ) + ) + + self.assertFalse( + check_perimeter_exceeds_grid_extents( + 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, -1000, -800, -700, -700, -700, -700, -800, -900, -1000, -1100] + ) + points_y = np.array([600, 700, 800, 900, 1000, 1100, 1000, 900, 800, 700, 600]) + expected_x = np.array([-1000, -800, -700, -700, -700, -800, -900, -1000]) + expected_y = np.array([700, 800, 900, 1000, 1000, 900, 800, 700]) + + granule_extent = {"x_min": -1000, "x_max": 1000, "y_min": -1000, "y_max": 1000} + + 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.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 From 77aedb750a6eaf7f14be05158dd602a5b662073b Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Tue, 21 Oct 2025 14:04:23 -0400 Subject: [PATCH 07/11] DAS-2424 Update comments and function names --- hoss/projection_utilities.py | 23 ++--- tests/unit/test_projection_utilities.py | 117 +++++++++++++++++------- 2 files changed, 92 insertions(+), 48 deletions(-) diff --git a/hoss/projection_utilities.py b/hoss/projection_utilities.py index 1bad149..f9d9d2c 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -490,20 +490,13 @@ def get_x_y_extents_from_geographic_perimeter( from_geo_transformer.transform(point_latitudes, point_longitudes) ) - finite_x, finite_y = remove_any_invalid_projected_values(points_x, 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 check_perimeter_exceeds_grid_extents( - finite_x, finite_y, granule_extent_projected_meters - ): - return { - 'x_min': granule_extent_projected_meters["x_min"], - 'x_max': granule_extent_projected_meters["x_max"], - 'y_min': granule_extent_projected_meters["y_min"], - 'y_max': granule_extent_projected_meters["y_max"], - } + 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( @@ -518,7 +511,7 @@ def get_x_y_extents_from_geographic_perimeter( } -def remove_any_invalid_projected_values( +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""" @@ -532,13 +525,11 @@ def remove_any_invalid_projected_values( return finite_x, finite_y -def check_perimeter_exceeds_grid_extents( +def perimeter_surrounds_grid( finite_x: np.ndarray, finite_y: np.ndarray, granule_extent: dict[str, float] ) -> bool: - """Check if perimeter exceeds the grid extents on all axes. If true, return - whole grid extents. This handles the case where the perimeter wholly - encloses the grid (e.g., whole-earth bbox, any lesser grid, polar - use-cases in particular) + """Returns True if perimeter exceeds the grid extents on all axes. + Returns False if does not. """ diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 6c88619..612d69a 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -27,7 +27,6 @@ MissingSpatialSubsetInformation, ) from hoss.projection_utilities import ( - check_perimeter_exceeds_grid_extents, get_bbox_polygon, get_densified_perimeter, get_filtered_points, @@ -45,7 +44,8 @@ get_x_y_extents_from_geographic_perimeter, is_projection_x_dimension, is_projection_y_dimension, - remove_any_invalid_projected_values, + perimeter_surrounds_grid, + remove_non_finite_projected_values, remove_points_outside_grid_extents, ) @@ -1245,7 +1245,7 @@ def test_get_x_y_extents_from_geographic_perimeter_full_earth_laea(self): self.assertTrue(all(not math.isinf(value) for value in x_y_extents.values())) - def test_remove_any_invalid_projected_values(self): + 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. """ @@ -1273,7 +1273,7 @@ def test_remove_any_invalid_projected_values(self): ] 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_any_invalid_projected_values(points_x, points_y) + 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) @@ -1282,22 +1282,36 @@ def test_check_perimeter_exceeds_grid_extents(self): extents and a False value is returned when the perimeter is within the grid extents """ - granule_extent = {"x_min": -1000, "x_max": 1000, "y_min": -1000, "y_max": 1000} - self.assertTrue( - check_perimeter_exceeds_grid_extents( - 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, + 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, + ) ) - ) - - self.assertFalse( - check_perimeter_exceeds_grid_extents( - 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, + 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 @@ -1305,24 +1319,63 @@ def test_remove_points_outside_grid_extents(self): """ points_x = np.array( - [-1100, -1000, -800, -700, -700, -700, -700, -800, -900, -1000, -1100] + [ + -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] ) - points_y = np.array([600, 700, 800, 900, 1000, 1100, 1000, 900, 800, 700, 600]) - expected_x = np.array([-1000, -800, -700, -700, -700, -800, -900, -1000]) - expected_y = np.array([700, 800, 900, 1000, 1000, 900, 800, 700]) - - granule_extent = {"x_min": -1000, "x_max": 1000, "y_min": -1000, "y_max": 1000} - 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) + granule_extent = { + "x_min": -1000.0, + "x_max": 1000.0, + "y_min": -1000.0, + "y_max": 1000.0, + } - 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, + 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): From f1e75ec272b6badf792981120312e6ad50223061 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Tue, 21 Oct 2025 14:42:21 -0400 Subject: [PATCH 08/11] DAS-2424 Correct Merge conflict --- tests/unit/test_projection_utilities.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 762ebf6..47280b5 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -1215,7 +1215,8 @@ def test_get_x_y_extents_from_geographic_perimeter(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_perimeter_full_earth_laea(self): From 5373aafccdf504e92cc473ed6c27a72b51c4dc8a Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Tue, 21 Oct 2025 15:14:32 -0400 Subject: [PATCH 09/11] DAS-2424 Update date in changelog and assert function in unit test --- CHANGELOG.md | 2 +- tests/unit/test_projection_utilities.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ef9bfb8..f233352 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,4 @@ -## [v1.1.13] - 2025-10-20 +## [v1.1.13] - 2025-10-21 ### Changed diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 47280b5..49bde5d 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -455,7 +455,7 @@ def test_get_projected_x_y_extents_edge_case(self): with self.subTest( 'LAEA - Bounding box which is close to the edge of granule extent' ): - self.assertDictEqual( + assert_float_dict_almost_equal( get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox), expected_output, ) From af686f7b6b544481cffb86acd129949a0c26c406 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Tue, 21 Oct 2025 17:18:43 -0400 Subject: [PATCH 10/11] DAS-2424 Add exception for invalid granule dimensions --- hoss/exceptions.py | 13 +++++++++++++ hoss/projection_utilities.py | 3 +++ tests/unit/test_projection_utilities.py | 12 ++++++++++++ 3 files changed, 28 insertions(+) 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 f9d9d2c..db28e8d 100644 --- a/hoss/projection_utilities.py +++ b/hoss/projection_utilities.py @@ -29,6 +29,7 @@ from hoss.bbox_utilities import BBox, flatten_list from hoss.exceptions import ( + InvalidGranuleDimensions, InvalidInputGeoJSON, InvalidRequestedRange, MissingGridMappingMetadata, @@ -204,6 +205,8 @@ 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 diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 49bde5d..6cbbae5 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -20,6 +20,7 @@ from hoss.bbox_utilities import BBox from hoss.exceptions import ( + InvalidGranuleDimensions, InvalidInputGeoJSON, InvalidRequestedRange, MissingGridMappingMetadata, @@ -455,6 +456,11 @@ def test_get_projected_x_y_extents_edge_case(self): with self.subTest( 'LAEA - Bounding box which is close to the edge of granule extent' ): + outputs = get_projected_x_y_extents( + x_values, y_values, crs, bounding_box=bbox + ) + print(f'outputs={outputs}') + print(f'expected_output={expected_output}') assert_float_dict_almost_equal( get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox), expected_output, @@ -465,6 +471,12 @@ def test_get_projected_x_y_extents_edge_case(self): 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. From da880a328ac46bc8530b6b64a03da35707e59781 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Tue, 21 Oct 2025 17:42:23 -0400 Subject: [PATCH 11/11] DAS-2424 Remove print statements --- tests/unit/test_projection_utilities.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tests/unit/test_projection_utilities.py b/tests/unit/test_projection_utilities.py index 6cbbae5..3758049 100644 --- a/tests/unit/test_projection_utilities.py +++ b/tests/unit/test_projection_utilities.py @@ -456,11 +456,6 @@ def test_get_projected_x_y_extents_edge_case(self): with self.subTest( 'LAEA - Bounding box which is close to the edge of granule extent' ): - outputs = get_projected_x_y_extents( - x_values, y_values, crs, bounding_box=bbox - ) - print(f'outputs={outputs}') - print(f'expected_output={expected_output}') assert_float_dict_almost_equal( get_projected_x_y_extents(x_values, y_values, crs, bounding_box=bbox), expected_output,