DAS-2437 Configure geographic spatial extent for SMAP polar grids

CHANGELOG.md

@@ -1,3 +1,9 @@
+## [v1.2.2] - 2026-01-22
+
+### Changed
+
+- SMAP Polar data subsetting constrained to North of Equator.
+
 ## [v1.2.1] - 2026-01-12
 
 ### Changed
@@ -231,6 +237,7 @@ Repository structure changes include:
 
 For more information on internal releases prior to NASA open-source approval,
 see legacy-CHANGELOG.md.
+[v1.2.2]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.2.2
 [v1.2.1]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.2.1
 [v1.2.0]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.2.0
 [v1.1.17]: https://github.com/nasa/harmony-opendap-subsetter/releases/tag/1.1.17

docker/service_version.txt

@@ -1 +1 @@
-1.2.1
+1.2.2

hoss/hoss_config.json

@@ -478,6 +478,10 @@
         {
           "Name": "master_geotransform",
           "Value": [-9000000, 3000, 0, 9000000, 0, -3000]
+        },
+        {
+          "Name": "geographic_spatial_extent",
+          "Value": [-180.0, 0, 180.0, 90.0]
         }
       ],
       "_Description": "Provide missing polar grid mapping attributes for SMAP L3 collections."
@@ -512,6 +516,10 @@
         {
           "Name": "master_geotransform",
           "Value": [-9000000, 9000, 0, 9000000, 0, -9000]
+        },
+        {
+          "Name": "geographic_spatial_extent",
+          "Value": [-180.0, 0, 180.0, 90.0]
         }
       ],
       "_Description": "Provide missing polar grid mapping attributes for SMAP L3 collections."
@@ -546,6 +554,10 @@
         {
           "Name": "master_geotransform",
           "Value": [-9000000, 36000, 0, 9000000, 0, -36000]
+        },
+        {
+          "Name": "geographic_spatial_extent",
+          "Value": [-180.0, 0, 180.0, 90.0]
         }
       ],
       "_Description": "Provide missing polar grid mapping attributes for SMAP L3 collections."

hoss/projection_utilities.py

@@ -102,6 +102,27 @@ def get_master_geotransform(
     )
 
 
+def get_geographic_spatial_extent(
+    variable: str, varinfo: VarInfoFromDmr
+) -> BBox | None:
+    """Retrieves the `geographic_spatial_extent` attribute from the grid mapping
+    attributes of the given variable. If the `geographic_spatial_extent` attribute
+    doesn't exist, a `None` value will be returned.
+
+    """
+    spatial_extent = get_grid_mapping_attributes(variable, varinfo).get(
+        "geographic_spatial_extent", None
+    )
+    if spatial_extent is not None:
+        return BBox(
+            spatial_extent[0],  # west
+            spatial_extent[1],  # south
+            spatial_extent[2],  # east
+            spatial_extent[3],  # north
+        )
+    return None
+
+
 def get_projected_x_y_variables(
     varinfo: VarInfoFromDmr, variable: str
 ) -> Tuple[Optional[str]]:
@@ -182,6 +203,7 @@ def get_projected_x_y_extents(
     crs: CRS,
     shape_file: str = None,
     bounding_box: BBox = None,
+    geographic_spatial_extent: BBox = None,
 ) -> Dict[str, float]:
     """Retrieve the minimum and maximum values for a projected grid as derived
     from either a bounding box or GeoJSON shape file, both of which are
@@ -202,6 +224,7 @@ def get_projected_x_y_extents(
                    'y_max': 5500}
 
     """
+
     grid_lats, grid_lons = get_grid_lat_lons(  # pylint: disable=unpacking-non-sequence
         x_values, y_values, crs
     )
@@ -219,13 +242,27 @@ def get_projected_x_y_extents(
         geographic_resolution, shape_file=shape_file, bounding_box=bounding_box
     )
 
-    # To avoid out-of-limits projection, we need to clip the bounding perimeter to
-    # the source file's geographic extents
     granule_extent = BBox(
         np.min(grid_lons), np.min(grid_lats), np.max(grid_lons), np.max(grid_lats)
     )
 
-    clipped_perimeter = get_filtered_points(densified_perimeter, granule_extent)
+    # If there is a configuration for geographic spatial extent apply that.
+    if geographic_spatial_extent is not None:
+        granule_extent = BBox(
+            np.max([granule_extent.west, geographic_spatial_extent.west]),
+            np.max([granule_extent.south, geographic_spatial_extent.south]),
+            np.min([granule_extent.east, geographic_spatial_extent.east]),
+            np.min([granule_extent.north, geographic_spatial_extent.north]),
+        )
+        requested_lons, requested_lats = np.array(densified_perimeter).T
+        clipped_lons, clipped_lats = remove_points_outside_grid_extents(
+            requested_lons, requested_lats, granule_extent
+        )
+        clipped_perimeter = list(zip(clipped_lons, clipped_lats))
+    else:
+        # To avoid out-of-limits projection, we need to clip the bounding perimeter to
+        # the source file's geographic extents
+        clipped_perimeter = get_filtered_points(densified_perimeter, granule_extent)
 
     granule_extent_projected_meters = {
         "x_min": np.min(x_values),
@@ -548,7 +585,9 @@ def perimeter_surrounds_grid(
 
 
 def remove_points_outside_grid_extents(
-    finite_x: np.ndarray, finite_y: np.ndarray, granule_extent: dict[str, float]
+    finite_x: np.ndarray,
+    finite_y: np.ndarray,
+    granule_extent: Union[BBox, 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.
@@ -559,12 +598,21 @@ def remove_points_outside_grid_extents(
     # The points are checked to make sure they are within
     # all 4 extents
 
-    mask = (
-        (finite_x >= granule_extent['x_min'] - tolerance)
-        & (finite_x <= granule_extent['x_max'] + tolerance)
-        & (finite_y >= granule_extent['y_min'] - tolerance)
-        & (finite_y <= granule_extent['y_max'] + tolerance)
-    )
+    if isinstance(granule_extent, BBox):
+        # finite_x is the requested lons and finite_y is the requested lats
+        mask = (
+            (finite_x >= granule_extent.west - tolerance)
+            & (finite_x <= granule_extent.east + tolerance)
+            & (finite_y >= granule_extent.south - tolerance)
+            & (finite_y <= granule_extent.north + tolerance)
+        )
+    else:
+        mask = (
+            (finite_x >= granule_extent['x_min'] - tolerance)
+            & (finite_x <= granule_extent['x_max'] + tolerance)
+            & (finite_y >= granule_extent['y_min'] - tolerance)
+            & (finite_y <= granule_extent['y_max'] + tolerance)
+        )
 
     finite_x = finite_x[mask]
     finite_y = finite_y[mask]

hoss/spatial.py

@@ -52,6 +52,7 @@
 )
 from hoss.exceptions import InvalidRequestedRange
 from hoss.projection_utilities import (
+    get_geographic_spatial_extent,
     get_master_geotransform,
     get_projected_x_y_extents,
     get_projected_x_y_variables,
@@ -288,13 +289,16 @@ def get_x_y_index_ranges_from_coordinates(
     projected_y, projected_x = dimension_arrays.keys()
 
     if not set((projected_x, projected_y)).issubset(set(index_ranges.keys())):
-
+        geographic_spatial_extent = get_geographic_spatial_extent(
+            non_spatial_variable, varinfo
+        )
         x_y_extents = get_projected_x_y_extents(
             dimension_arrays[projected_x][:],
             dimension_arrays[projected_y][:],
             crs,
             shape_file=shape_file_path,
             bounding_box=bounding_box,
+            geographic_spatial_extent=geographic_spatial_extent,
         )
 
         x_index_ranges = get_dimension_index_range(

tests/unit/test_projection_utilities.py

@@ -32,6 +32,7 @@
     get_densified_perimeter,
     get_filtered_points,
     get_geographic_resolution,
+    get_geographic_spatial_extent,
     get_grid_lat_lons,
     get_grid_mapping_attributes,
     get_master_geotransform,
@@ -385,7 +386,7 @@ def test_get_projected_x_y_extents(self):
     def test_get_projected_x_y_extents_whole_earth(self):
         """Ensure that the expected values for the x and y dimension extents
         are recovered for a polar projected grid and when a whole earth
-        bounding box or shape is requested.
+        bounding box or shape is requested and geographic extent is not constrained.
 
         """
         whole_earth_bbox = BBox(-180.0, -90.0, 180.0, 90.0)
@@ -472,6 +473,75 @@ def test_get_projected_x_y_extents_edge_case(self):
                 y_values1 = np.linspace(9200000, -9200000, 500)
                 get_projected_x_y_extents(x_values1, y_values1, crs, bounding_box=bbox)
 
+    def test_get_projected_x_y_extents_with_configured_geographic_extent(self):
+        """Ensure that the expected values for the x and y dimension extents
+        are returned when the geographic spatial extent is configured
+
+        The dimension values used below mimic the SPL3FTP collection's polar
+        grid. The values returned should not contain extents outside the
+        configured spatial extent.
+
+        """
+
+        x_values = np.linspace(-8982000, 8982000, 500)
+        y_values = np.linspace(8982000, -8982000, 500)
+        crs = CRS.from_cf(
+            {
+                'false_easting': 0.0,
+                'false_northing': 0.0,
+                'latitude_of_projection_origin': 90.0,
+                'longitude_of_projection_origin': 0.0,
+                'grid_mapping_name': 'lambert_azimuthal_equal_area',
+            }
+        )
+
+        bounding_box = BBox(12, -38, 36, 68)
+        geographic_spatial_extent = BBox(-180.0, 0, 180, 90.0)
+
+        expected_output = {
+            'x_min': 507518.9840003274,
+            'x_max': 5288134.075113788,
+            'y_min': -8800111.32936258,
+            'y_max': -1974835.9575607865,
+        }
+
+        with self.subTest('geo spatial extent configured'):
+            geographic_spatial_extent = BBox(-180.0, 0, 180, 90.0)
+            expected_output = {
+                'x_min': 507518.9840003274,
+                'x_max': 5288134.075113788,
+                'y_min': -8800111.32936258,
+                'y_max': -1974835.9575607865,
+            }
+            assert_float_dict_almost_equal(
+                get_projected_x_y_extents(
+                    x_values,
+                    y_values,
+                    crs,
+                    bounding_box=bounding_box,
+                    geographic_spatial_extent=geographic_spatial_extent,
+                ),
+                expected_output,
+            )
+
+        with self.subTest('geo spatial extent not configured'):
+            expected_output = {
+                'x_min': 507518.9840003274,
+                'x_max': 6525593.1802023165,
+                'y_min': -8981708.47358046,
+                'y_max': -1974835.9575607865,
+            }
+            assert_float_dict_almost_equal(
+                get_projected_x_y_extents(
+                    x_values,
+                    y_values,
+                    crs,
+                    bounding_box=bounding_box,
+                    geographic_spatial_extent=None,
+                ),
+                expected_output,
+            )
+
     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.
@@ -1473,3 +1543,29 @@ def test_get_master_geotransform(self, mock_get_grid_mapping_attributes):
             }
             result = get_master_geotransform("test_variable", varinfo)
             self.assertIsNone(result)
+
+    def test_get_geographic_spatial_extent(self):
+        """Ensure that the `geographic_spatial_extent` attribute is returned if
+        it exists. If it doesn't exist the return value should be `None`.
+
+        """
+
+        varinfo = VarInfoFromDmr(
+            'tests/data/SC_SPL3FTP_004.dmr',
+            'SPL3FTP',
+            'hoss/hoss_config.json',
+        )
+
+        with self.subTest('grid mapping attribute contains geographic spatial extent'):
+            result = get_geographic_spatial_extent(
+                "/Freeze_Thaw_Retrieval_Data_Polar/altitude_dem", varinfo
+            )
+            self.assertEqual(result, BBox(west=-180.0, south=0, east=180.0, north=90.0))
+
+        with self.subTest(
+            'grid mapping attribute does not contain geographic spatial extent'
+        ):
+            result = get_geographic_spatial_extent(
+                "/Freeze_Thaw_Retrieval_Data_Global/altitude_dem", varinfo
+            )
+            self.assertIsNone(result)