DAS-2278 - Handle spatial subsetting for SMAP L3 collections with pol…

…ar grid_mapping variables (#26) * Add master_geotransform attribute to polar grid_mapping variables * Update master_geotransform values * DAS-2278 - Create dimension arrays from geotransform * DAS-2278 - Add --platform linux/amd64 to Docker commands for compatibility * Update ICESat2 config item to align with recent trajectory subsetter config updates * Updates for draft PR comments * Modify for 3D variable cases * Modify functions related to getting crs and geotransform * Updates to varinfo config for errors discovered in testing * Updated unit tests for new and updated functions * Update CHANGELOG and docker service version * Update comment to clarify purpose of varinfo.get_missing_variable_attributes * Clarify CHANGELOG entry to better reflect purpose of the change
nasa · Feb 12, 2025 · 0437f41 · 0437f41
1 parent 212fda9
commit 0437f41
Show file tree

Hide file tree

Showing 11 changed files with 457 additions and 63 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,3 +1,12 @@
+## v1.1.4
+### 2025-02-12
+
+This version of HOSS adds support for SMAP L3 polar variables that are unable to have their
+dimension scale arrays created from their corresponding lat/lon variables. A 'master geotransform'
+attribute has been added to the grid mapping reference variable for the affected collections
+and function updates were made to create the dimension arrays from the master geotransform
+when it is present.
+
 ## v1.1.3
 ### 2025-01-29
 

diff --git a/bin/build-image b/bin/build-image
@@ -27,4 +27,4 @@ bin/clean-images
 #   version number from `docker/service_version.txt`.
 # - "latest", so the test Dockerfile can use the service image as a base image.
 #
-docker build -t ${image}:${tag} -t ${image}:latest -f docker/service.Dockerfile .
+docker build --platform linux/amd64 -t ${image}:${tag} -t ${image}:latest -f docker/service.Dockerfile .
diff --git a/bin/build-test b/bin/build-test
@@ -21,4 +21,4 @@ if [ ! -z "$old" ] && [ "$2" != "--no-delete" ]; then
 fi
 
 # Build the image
-docker build -t ${image}:${tag} -f docker/tests.Dockerfile .
+docker build --platform linux/amd64 -t ${image}:${tag} -f docker/tests.Dockerfile .
diff --git a/bin/run-test b/bin/run-test
@@ -23,7 +23,7 @@ mkdir -p coverage
 
 # Run the tests in a Docker container with mounted volumes for XML report
 # output and test coverage reporting
-docker run --rm \
+docker run --platform linux/amd64 --rm \
 	-v $(pwd)/test-reports:/home/tests/reports \
 	-v $(pwd)/coverage:/home/tests/coverage \
 	ghcr.io/nasa/harmony-opendap-subsetter-test "$@"
diff --git a/docker/service_version.txt b/docker/service_version.txt
@@ -1 +1 @@
-1.1.3
+1.1.4
diff --git a/hoss/coordinate_utilities.py b/hoss/coordinate_utilities.py
@@ -440,3 +440,44 @@ def interpolate_dim_values_from_sample_pts(
     dim_max = dim_values[1] + (dim_resolution * (dim_size - 1 - dim_indices[1]))
 
     return np.linspace(dim_min, dim_max, dim_size)
+
+
+def create_dimension_arrays_from_geotransform(
+    prefetch_dataset: Dataset,
+    latitude_coordinate: VariableFromDmr,
+    projected_dimension_names: list[str],
+    geotranform,
+) -> dict[str, np.ndarray]:
+    """Generate artificial 1D dimensions scales from geotransform"""
+    lat_arr = get_2d_coordinate_array(
+        prefetch_dataset,
+        latitude_coordinate.full_name_path,
+    )
+
+    # compute the x,y locations along a column and row
+    column_dimensions = [
+        col_row_to_xy(geotranform, col, 0) for col in range(lat_arr.shape[-1])
+    ]
+    row_dimensions = [
+        col_row_to_xy(geotranform, 0, row) for row in range(lat_arr.shape[-2])
+    ]
+
+    # pull out dimension values
+    x_values = np.array([x for x, y in column_dimensions], dtype=np.float64)
+    y_values = np.array([y for x, y in row_dimensions], dtype=np.float64)
+    projected_y, projected_x = projected_dimension_names[-2:]
+
+    return {projected_y: y_values, projected_x: x_values}
+
+
+def col_row_to_xy(geotransform, col: int, row: int) -> tuple[np.float64, np.float64]:
+    """Convert grid cell location to x,y coordinate."""
+    # Geotransform is defined from upper left corner as (0,0), so adjust
+    # input value to the center of grid at (.5, .5)
+    adj_col = col + 0.5
+    adj_row = row + 0.5
+
+    x = geotransform[0] + adj_col * geotransform[1] + adj_row * geotransform[2]
+    y = geotransform[3] + adj_col * geotransform[4] + adj_row * geotransform[5]
+
+    return x, y
diff --git a/hoss/hoss_config.json b/hoss/hoss_config.json
@@ -1,6 +1,6 @@
 {
   "Identification": "hoss_config",
-  "Version": 20,
+  "Version": 21,
   "CollectionShortNamePath": [
     "/HDF5_GLOBAL/short_name",
     "/NC_GLOBAL/short_name",
@@ -119,13 +119,27 @@
     {
       "Applicability": {
         "Mission": "SMAP",
-        "ShortNamePath": "SPL3FT(P|P_E)",
-        "VariablePattern": "(?i).*polar.*"
+        "ShortNamePath": "SPL3FTP",
+        "VariablePattern": "/Freeze_Thaw_Retrieval_Data_Polar/.*"
       },
       "Attributes": [
         {
           "Name": "grid_mapping",
-          "Value": "/EASE2_polar_projection"
+          "Value": "/EASE2_polar_projection_36km"
+        }
+      ],
+      "_Description": "SMAP L3 collections omit polar grid mapping information"
+    },
+    {
+      "Applicability": {
+        "Mission": "SMAP",
+        "ShortNamePath": "SPL3FTP_E",
+        "VariablePattern": "/Freeze_Thaw_Retrieval_Data_Polar/.*"
+      },
+      "Attributes": [
+        {
+          "Name": "grid_mapping",
+          "Value": "/EASE2_polar_projection_9km"
         }
       ],
       "_Description": "SMAP L3 collections omit polar grid mapping information"
@@ -134,7 +148,7 @@
       "Applicability": {
         "Mission": "SMAP",
         "ShortNamePath": "SPL3SMP_E",
-        "VariablePattern": "Soil_Moisture_Retrieval_Data_(A|P)M/.*"
+        "VariablePattern": "/Soil_Moisture_Retrieval_Data_(A|P)M/.*"
       },
       "Attributes": [
         {
@@ -148,12 +162,12 @@
       "Applicability": {
         "Mission": "SMAP",
         "ShortNamePath": "SPL3SMP_E",
-        "VariablePattern": "Soil_Moisture_Retrieval_Data_Polar_(A|P)M/.*"
+        "VariablePattern": "/Soil_Moisture_Retrieval_Data_Polar_(A|P)M/.*"
       },
       "Attributes": [
         {
           "Name": "grid_mapping",
-          "Value": "/EASE2_polar_projection"
+          "Value": "/EASE2_polar_projection_9km"
         }
       ],
       "_Description": "SMAP L3 collections omit polar grid mapping information"
@@ -166,15 +180,15 @@
       "Attributes": [
         {
           "Name": "grid_mapping",
-          "Value": "/EASE2_polar_projection"
+          "Value": "/EASE2_polar_projection_3km"
         }
       ],
       "_Description": "SMAP L3 collections omit polar grid mapping information"
     },
     {
       "Applicability": {
         "Mission": "SMAP",
-        "ShortNamePath": "SPL3SM(P|A|AP)|SPL2SMAP_S"
+        "ShortNamePath": "SPL3SM(P|A|AP)$|SPL2SMAP_S"
       },
       "Attributes": [
         {
@@ -217,8 +231,76 @@
     {
       "Applicability": {
         "Mission": "SMAP",
-        "ShortNamePath": "SPL3FT(A|P|P_E)|SPL3SM(P|P_E|A|AP)|SPL2SMAP_S",
-        "VariablePattern": "/EASE2_polar_projection"
+        "ShortNamePath": "SPL3FTA",
+        "VariablePattern": "/EASE2_polar_projection_3km"
+      },
+      "Attributes": [
+        {
+          "Name": "grid_mapping_name",
+          "Value": "lambert_azimuthal_equal_area"
+        },
+        {
+          "Name": "longitude_of_projection_origin",
+          "Value": 0.0
+        },
+        {
+          "Name": "latitude_of_projection_origin",
+          "Value": 90.0
+        },
+        {
+          "Name": "false_easting",
+          "Value": 0.0
+        },
+        {
+          "Name": "false_northing",
+          "Value": 0.0
+        },
+        {
+          "Name": "master_geotransform",
+          "Value": [-9000000, 3000, 0, 9000000, 0, -3000]
+        }
+      ],
+      "_Description": "Provide missing polar grid mapping attributes for SMAP L3 collections."
+    },
+    {
+      "Applicability": {
+        "Mission": "SMAP",
+        "ShortNamePath": "SPL3FTP_E|SPL3SMP_E",
+        "VariablePattern": "/EASE2_polar_projection_9km"
+      },
+      "Attributes": [
+        {
+          "Name": "grid_mapping_name",
+          "Value": "lambert_azimuthal_equal_area"
+        },
+        {
+          "Name": "longitude_of_projection_origin",
+          "Value": 0.0
+        },
+        {
+          "Name": "latitude_of_projection_origin",
+          "Value": 90.0
+        },
+        {
+          "Name": "false_easting",
+          "Value": 0.0
+        },
+        {
+          "Name": "false_northing",
+          "Value": 0.0
+        },
+        {
+          "Name": "master_geotransform",
+          "Value": [-9000000, 9000, 0, 9000000, 0, -9000]
+        }
+      ],
+      "_Description": "Provide missing polar grid mapping attributes for SMAP L3 collections."
+    },
+    {
+      "Applicability": {
+        "Mission": "SMAP",
+        "ShortNamePath": "SPL3FTP",
+        "VariablePattern": "/EASE2_polar_projection_36km"
       },
       "Attributes": [
         {
@@ -240,6 +322,10 @@
         {
           "Name": "false_northing",
           "Value": 0.0
+        },
+        {
+          "Name": "master_geotransform",
+          "Value": [-9000000, 36000, 0, 9000000, 0, -36000]
         }
       ],
       "_Description": "Provide missing polar grid mapping attributes for SMAP L3 collections."
@@ -318,7 +404,7 @@
       "Applicability": {
         "Mission": "SMAP",
         "ShortNamePath": "SPL3FT(A|P|P_E)",
-        "VariablePattern": "^/Freeze_Thaw_Retrieval_Data(?:_(Global|Polar))/(transition_direction$|transition_state_flag$)"
+        "VariablePattern": "^/Freeze_Thaw_Retrieval_Data(?:_(Global|Polar))?/(transition_direction$|transition_state_flag$)"
       },
       "Attributes": [
         {
@@ -332,7 +418,7 @@
       "Applicability": {
         "Mission": "SMAP",
         "ShortNamePath": "SPL3FT(A|P|P_E)",
-        "VariablePattern": "^/Freeze_Thaw_Retrieval_Data(?:_(Global|Polar))/((?!transition_state_flag$)(?!transition_direction$).)*$|/Radar_Data/.*"
+        "VariablePattern": "^/Freeze_Thaw_Retrieval_Data(?:_(Global|Polar))?/((?!transition_state_flag$)(?!transition_direction$).)*$|/Radar_Data/.*"
       },
       "Attributes": [
         {
@@ -526,7 +612,7 @@
     {
       "Applicability": {
         "Mission": "ICESat2",
-        "ShortNamePath": "ATL0[3-9]|ATL1[023]",
+        "ShortNamePath": "ATL03",
         "VariablePattern": "/gt[123][lr]/geolocation/.*"
       },
       "Attributes": [

diff --git a/hoss/projection_utilities.py b/hoss/projection_utilities.py
@@ -41,9 +41,16 @@
 
 
 def get_variable_crs(variable: str, varinfo: VarInfoFromDmr) -> CRS:
+    """Retrieves the grid mapping variable metadata attributes for a given
+    variable and creates a `pyproj.CRS` object from the grid mapping attributes.
+
+    """
+    return CRS.from_cf(get_grid_mapping_attributes(variable, varinfo))
+
+
+def get_grid_mapping_attributes(variable: str, varinfo: VarInfoFromDmr) -> Dict:
     """Check the metadata attributes for the variable to find the associated
-    grid mapping variable. Create a `pyproj.CRS` object from the grid
-    mapping variable metadata attributes.
+    grid mapping variable.
 
     All metadata attributes that contain references from one variable to
     another are stored in the `Variable.references` dictionary attribute
@@ -66,11 +73,11 @@ def get_variable_crs(variable: str, varinfo: VarInfoFromDmr) -> CRS:
             if grid_mapping_variable is not None:
                 cf_attributes = grid_mapping_variable.attributes
             else:
-                # check for any overrides
+                # check for configuration provided attributes
                 cf_attributes = varinfo.get_missing_variable_attributes(grid_mapping)
 
             if cf_attributes:
-                crs = CRS.from_cf(cf_attributes)
+                return cf_attributes
             else:
                 raise MissingGridMappingVariable(grid_mapping, variable)
 
@@ -80,7 +87,18 @@ def get_variable_crs(variable: str, varinfo: VarInfoFromDmr) -> CRS:
     else:
         raise MissingGridMappingMetadata(variable)
 
-    return crs
+
+def get_master_geotransform(
+    variable: str, varinfo: VarInfoFromDmr
+) -> Optional[List[int]]:
+    """Retrieves the `master_geotransform` attribute from the grid mapping
+    attributes of the given variable. If the `master_geotransform` attribute
+    doesn't exist, a `None` value will be returned.
+
+    """
+    return get_grid_mapping_attributes(variable, varinfo).get(
+        "master_geotransform", None
+    )
 
 
 def get_projected_x_y_variables(

diff --git a/hoss/spatial.py b/hoss/spatial.py
@@ -37,6 +37,7 @@
 )
 from hoss.coordinate_utilities import (
     create_dimension_arrays_from_coordinates,
+    create_dimension_arrays_from_geotransform,
     get_coordinate_variables,
     get_dimension_array_names,
     get_variables_with_anonymous_dims,
@@ -49,6 +50,7 @@
     get_dimension_index_range,
 )
 from hoss.projection_utilities import (
+    get_master_geotransform,
     get_projected_x_y_extents,
     get_projected_x_y_variables,
     get_variable_crs,
@@ -248,14 +250,22 @@ def get_x_y_index_ranges_from_coordinates(
     crs = get_variable_crs(non_spatial_variable, varinfo)
 
     projected_dimension_names = get_dimension_array_names(varinfo, non_spatial_variable)
-
-    dimension_arrays = create_dimension_arrays_from_coordinates(
-        prefetch_coordinate_datasets,
-        latitude_coordinate,
-        longitude_coordinate,
-        crs,
-        projected_dimension_names,
-    )
+    master_geotransform = get_master_geotransform(non_spatial_variable, varinfo)
+    if master_geotransform:
+        dimension_arrays = create_dimension_arrays_from_geotransform(
+            prefetch_coordinate_datasets,
+            latitude_coordinate,
+            projected_dimension_names,
+            master_geotransform,
+        )
+    else:
+        dimension_arrays = create_dimension_arrays_from_coordinates(
+            prefetch_coordinate_datasets,
+            latitude_coordinate,
+            longitude_coordinate,
+            crs,
+            projected_dimension_names,
+        )
 
     projected_y, projected_x = dimension_arrays.keys()