DAS-2232 - Functionality added to support SMAP L3 products

CHANGELOG.md

@@ -1,3 +1,9 @@
+## v1.1.0
+### 2024-09-10
+
+This version of HOSS provides support for products without CF compliance like SMAP L3
+Methods added to get dimension scales from coordinate attributes and grid mapping with overrides
+in the json file
 
 ## v1.0.5
 ### 2024-08-19

docker/service_version.txt

@@ -1 +1 @@
-1.0.5
+1.1.0

docs/requirements.txt

@@ -16,5 +16,5 @@
 #
 harmony-py~=0.4.10
 netCDF4~=1.6.4
-notebook~=7.0.4
+notebook~=7.2.2
 xarray~=2023.9.0

hoss/dimension_utilities.py

@@ -19,11 +19,15 @@
 from harmony.message_utility import rgetattr
 from harmony.util import Config
 from netCDF4 import Dataset
+from numpy import ndarray
 from numpy.ma.core import MaskedArray
 from varinfo import VariableFromDmr, VarInfoFromDmr
 
-from hoss.bbox_utilities import flatten_list
 from hoss.exceptions import InvalidNamedDimension, InvalidRequestedRange
+from hoss.projection_utilities import (
+    get_variable_crs,
+    get_x_y_extents_from_geographic_points,
+)
 from hoss.utilities import (
     format_variable_set_string,
     get_opendap_nc4,
@@ -55,6 +59,42 @@ def is_index_subset(message: Message) -> bool:
     )
 
 
+def get_override_projected_dimensions(
+    varinfo: VarInfoFromDmr,
+    override_variable_name: str,
+) -> str:
+    """returns the x-y projection variable names that would
+    match the geo coordinate names
+
+    """
+    projection_variable_name = None
+    override_variable = varinfo.get_variable(override_variable_name)
+    if override_variable.is_latitude():
+        projection_variable_name = 'projected_y'
+    elif override_variable.is_longitude():
+        projection_variable_name = 'projected_x'
+    return projection_variable_name
+
+
+def get_override_dimensions(
+    varinfo: VarInfoFromDmr,
+    required_variables: Set[str],
+) -> set[str]:
+    """Determine the dimensions that need to be "pre-fetched" from OPeNDAP
+    If dimensions are not available, get variables that can be
+    used to generate dimensions
+    """
+    # check if coordinate variables are provided
+    # this should be be a configurable in hoss_config.json
+    try:
+        override_dimensions = varinfo.get_references_for_attribute(
+            required_variables, 'coordinates'
+        )
+        return override_dimensions
+    except AttributeError:
+        return set()
+
+
 def prefetch_dimension_variables(
     opendap_url: str,
     varinfo: VarInfoFromDmr,
@@ -63,7 +103,7 @@ def prefetch_dimension_variables(
     logger: Logger,
     access_token: str,
     config: Config,
-) -> str:
+) -> tuple[str, set[str]] | str:
     """Determine the dimensions that need to be "pre-fetched" from OPeNDAP in
     order to derive index ranges upon them. Initially, this was just
     spatial and temporal dimensions, but to support generic dimension
@@ -73,38 +113,133 @@ def prefetch_dimension_variables(
 
     """
     required_dimensions = varinfo.get_required_dimensions(required_variables)
+    if len(required_dimensions) == 0:
+        required_dimensions = get_override_dimensions(varinfo, required_variables)
+        logger.info('coordinates: ' f'{required_dimensions}')
 
-    # Iterate through all requested dimensions and extract a list of bounds
-    # references for each that has any. This will produce a list of lists,
-    # which should be flattened into a single list and then combined into a set
-    # to remove duplicates.
-    bounds = set(
-        flatten_list(
-            [
-                list(varinfo.get_variable(dimension).references.get('bounds'))
-                for dimension in required_dimensions
-                if varinfo.get_variable(dimension).references.get('bounds') is not None
-            ]
-        )
-    )
-
+    bounds = varinfo.get_references_for_attribute(required_dimensions, 'bounds')
+    logger.info('bounds: ' f'{bounds}')
     required_dimensions.update(bounds)
 
     logger.info(
         'Variables being retrieved in prefetch request: '
         f'{format_variable_set_string(required_dimensions)}'
     )
-
     required_dimensions_nc4 = get_opendap_nc4(
         opendap_url, required_dimensions, output_dir, logger, access_token, config
     )
-
     # Create bounds variables if necessary.
     add_bounds_variables(required_dimensions_nc4, required_dimensions, varinfo, logger)
-
     return required_dimensions_nc4
 
 
+def is_variable_one_dimensional(
+    prefetch_dataset: Dataset, dimension_variable: VariableFromDmr
+) -> bool:
+    """Check if a dimension variable is 1D"""
+    dimensions_array = prefetch_dataset[dimension_variable.full_name_path][:]
+    return dimensions_array.ndim == 1
+
+
+def update_dimension_variables(
+    prefetch_dataset: Dataset,
+    required_dimensions: Set[str],
+    varinfo: VarInfoFromDmr,
+) -> Dict[str, ndarray]:
+    """Augment a NetCDF4 file with artificial 1D dimensions variable for each
+    2D dimension variable"
+
+    For each dimension variable:
+    (1) Check if the dimension variable is 1D.
+    (2) If it is not 1D and is 2D create a dimensions array from
+        within the `write_1D_dimensions`
+        function.
+    (3) Then write the 1D dimensions variable to the NetCDF4 URL.
+
+    """
+    for dimension_name in required_dimensions:
+        dimension_variable = varinfo.get_variable(dimension_name)
+        if not is_variable_one_dimensional(prefetch_dataset, dimension_variable):
+            col_size = prefetch_dataset[dimension_variable.full_name_path][:].shape[0]
+            row_size = prefetch_dataset[dimension_variable.full_name_path][:].shape[1]
+            crs = get_variable_crs(dimension_name, varinfo)
+
+    geo_grid_corners = get_geo_grid_corners(
+        prefetch_dataset, required_dimensions, varinfo
+    )
+
+    x_y_extents = get_x_y_extents_from_geographic_points(geo_grid_corners, crs)
+
+    # get grid size and resolution
+    x_min = x_y_extents['x_min']
+    x_max = x_y_extents['x_max']
+    y_min = x_y_extents['y_min']
+    y_max = x_y_extents['y_max']
+    x_resolution = (x_max - x_min) / row_size
+    y_resolution = (y_max - y_min) / col_size
+
+    # create the xy dim scales
+    x_dim = np.arange(x_min, x_max, x_resolution)
+
+    # ascending versus descending..should be based on the coordinate grid
+    y_dim = np.arange(y_max, y_min, -y_resolution)
+    return {'projected_y': y_dim, 'projected_x': x_dim}
+
+
+def get_geo_grid_corners(
+    prefetch_dataset: Dataset,
+    required_dimensions: Set[str],
+    varinfo: VarInfoFromDmr,
+) -> list[Tuple[float]]:
+    """
+    This method is used to return the lat lon corners from a 2D
+    coordinate dataset. This does a check for values below -180
+    which could be fill values. Does not check if there are fill
+    values in the corner points to go down to the next row and col
+
+    """
+    for dimension_name in required_dimensions:
+        dimension_variable = varinfo.get_variable(dimension_name)
+        if dimension_variable.is_latitude():
+            lat_arr = prefetch_dataset[dimension_variable.full_name_path][:]
+        elif dimension_variable.is_longitude():
+            lon_arr = prefetch_dataset[dimension_variable.full_name_path][:]
+
+    # skip fill values when calculating min values
+    # topleft =  minlon, maxlat
+    # bottomright = maxlon, minlat
+    top_left_row_idx = 0
+    top_left_col_idx = 0
+    lon_row = lon_arr[top_left_row_idx, :]
+    lon_row_valid_indices = np.where(lon_row >= -180.0)[0]
+    top_left_col_idx = lon_row_valid_indices[lon_row[lon_row_valid_indices].argmin()]
+    minlon = lon_row[top_left_col_idx]
+    top_right_col_idx = lon_row_valid_indices[lon_row[lon_row_valid_indices].argmax()]
+    maxlon = lon_row[top_right_col_idx]
+
+    lat_col = lat_arr[:, top_right_col_idx]
+    lat_col_valid_indices = np.where(lat_col >= -180.0)[0]
+    bottom_right_row_idx = lat_col_valid_indices[
+        lat_col[lat_col_valid_indices].argmin()
+    ]
+    minlat = lat_col[bottom_right_row_idx]
+    top_right_row_idx = lat_col_valid_indices[lat_col[lat_col_valid_indices].argmax()]
+    maxlat = lat_col[top_right_row_idx]
+
+    topleft_corner = [minlon, maxlat]
+    topright_corner = [maxlon, maxlat]
+    bottomright_corner = [maxlon, minlat]
+    bottomleft_corner = [minlon, minlat]
+
+    geo_grid_corners = [
+        topleft_corner,
+        topright_corner,
+        bottomright_corner,
+        bottomleft_corner,
+    ]
+    return geo_grid_corners
+
+
 def add_bounds_variables(
     dimensions_nc4: str,
     required_dimensions: Set[str],
@@ -409,7 +544,9 @@ def get_dimension_indices_from_bounds(
 
 
 def add_index_range(
-    variable_name: str, varinfo: VarInfoFromDmr, index_ranges: IndexRanges
+    variable_name: str,
+    varinfo: VarInfoFromDmr,
+    index_ranges: IndexRanges,
 ) -> str:
     """Append the index ranges of each dimension for the specified variable.
     If there are no dimensions with listed index ranges, then the full
@@ -422,22 +559,48 @@ def add_index_range(
 
     """
     variable = varinfo.get_variable(variable_name)
-
     range_strings = []
-
-    for dimension in variable.dimensions:
-        dimension_range = index_ranges.get(dimension)
-
-        if dimension_range is not None and dimension_range[0] <= dimension_range[1]:
-            range_strings.append(f'[{dimension_range[0]}:{dimension_range[1]}]')
+    if variable.dimensions == []:
+        override_dimensions = get_override_dimensions(varinfo, [variable_name])
+        if len(override_dimensions) > 0:
+            for override in reversed(list(override_dimensions)):
+                dimension = get_override_projected_dimensions(varinfo, override)
+                dimension_range = index_ranges.get(dimension)
+                if (
+                    dimension_range is not None
+                    and dimension_range[0] <= dimension_range[1]
+                ):
+                    range_strings.append(f'[{dimension_range[0]}:{dimension_range[1]}]')
+                else:
+                    range_strings.append('[]')
         else:
-            range_strings.append('[]')
+            # if the override is the variable
+            override = get_override_projected_dimensions(varinfo, variable_name)
+            if override is not None and override in index_ranges.keys():
+                for dimension in reversed(index_ranges.keys()):
+                    dimension_range = index_ranges.get(dimension)
+                    if (
+                        dimension_range is not None
+                        and dimension_range[0] <= dimension_range[1]
+                    ):
+                        range_strings.append(
+                            f'[{dimension_range[0]}:{dimension_range[1]}]'
+                        )
+            else:
+                range_strings.append('[]')
+    else:
+        for dimension in variable.dimensions:
+            dimension_range = index_ranges.get(dimension)
+
+            if dimension_range is not None and dimension_range[0] <= dimension_range[1]:
+                range_strings.append(f'[{dimension_range[0]}:{dimension_range[1]}]')
+            else:
+                range_strings.append('[]')
 
     if all(range_string == '[]' for range_string in range_strings):
         indices_string = ''
     else:
         indices_string = ''.join(range_strings)
-
     return f'{variable_name}{indices_string}'
 
 
@@ -549,7 +712,10 @@ def get_dimension_bounds(
     be returned.
 
     """
-    bounds = varinfo.get_variable(dimension_name).references.get('bounds')
+    try:
+        bounds = varinfo.get_variable(dimension_name).references.get('bounds')
+    except AttributeError:
+        bounds = None
 
     if bounds is not None:
         try: