DAS-2424 Update calculation of projected extents from requested spatial range to remove points outside the granule extent

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

docker/service_version.txt

@@ -1 +1 @@
-1.1.12
+1.1.13

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_perimeter(
+        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
@@ -458,13 +464,15 @@ 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, 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. 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)
@@ -478,6 +486,45 @@ 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)
+        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),
+        }
+    # 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
+
     return {
         'x_min': np.min(finite_x),
         'x_max': np.max(finite_x),

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,
@@ -40,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,
 )
@@ -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.
@@ -1145,12 +1183,14 @@ 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.
 
         """
+        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,15 +1201,18 @@ 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_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.
 
         """
+        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_perimeter(
+            points1, crs, x_values, y_values
+        )
 
         self.assertTrue(all(not math.isinf(value) for value in x_y_extents.values()))