DAS-2232 - added get_x_y_index_ranges_from_coordinates as a separate …

…method

hoss/dimension_utilities.py

@@ -21,16 +21,17 @@
 from netCDF4 import Dataset
 from numpy import ndarray
 from numpy.ma.core import MaskedArray
+from pyproj import CRS
 from varinfo import VariableFromDmr, VarInfoFromDmr
 
 from hoss.exceptions import (
     InvalidNamedDimension,
     InvalidRequestedRange,
+    IrregularCoordinateDatasets,
     MissingCoordinateDataset,
     MissingValidCoordinateDataset,
 )
 from hoss.projection_utilities import (
-    get_variable_crs,
     get_x_y_extents_from_geographic_points,
 )
 from hoss.utilities import (
@@ -83,12 +84,13 @@ def get_prefetch_variables(
 
     """
     required_dimensions = varinfo.get_required_dimensions(required_variables)
-    if not required_dimensions:
+    if required_dimensions:
+        bounds = varinfo.get_references_for_attribute(required_dimensions, 'bounds')
+        required_dimensions.update(bounds)
+    else:
         coordinate_variables = get_coordinate_variables(varinfo, required_variables)
-        required_dimensions = set(coordinate_variables)
-
-    bounds = varinfo.get_references_for_attribute(required_dimensions, 'bounds')
-    required_dimensions.update(bounds)
+        if coordinate_variables:
+            required_dimensions = set(coordinate_variables)
 
     logger.info(
         'Variables being retrieved in prefetch request: '
@@ -152,9 +154,7 @@ def get_coordinate_variables(
 
 
 def update_dimension_variables(
-    prefetch_dataset: Dataset,
-    coordinates: Set[str],
-    varinfo: VarInfoFromDmr,
+    prefetch_dataset: Dataset, coordinates: Set[str], varinfo: VarInfoFromDmr, crs: CRS
 ) -> Dict[str, ndarray]:
     """Generate artificial 1D dimensions variable for each
     2D dimension or coordinate variable
@@ -167,12 +167,9 @@ def update_dimension_variables(
     (5) Generate the x-y dimscale array and return to the calling method
 
     """
-    for coordinate_name in coordinates:
-        coordinate_variable = varinfo.get_variable(coordinate_name)
-        if prefetch_dataset[coordinate_variable.full_name_path][:].ndim > 1:
-            col_size = prefetch_dataset[coordinate_variable.full_name_path][:].shape[0]
-            row_size = prefetch_dataset[coordinate_variable.full_name_path][:].shape[1]
-        crs = get_variable_crs(coordinate_name, varinfo)
+    row_size, col_size = get_row_col_sizes_from_coordinate_datasets(
+        prefetch_dataset, coordinates, varinfo
+    )
 
     geo_grid_corners = get_geo_grid_corners(prefetch_dataset, coordinates, varinfo)
 
@@ -198,6 +195,27 @@ def update_dimension_variables(
     return {'projected_y': y_dim, 'projected_x': x_dim}
 
 
+def get_row_col_sizes_from_coordinate_datasets(
+    prefetch_dataset: Dataset,
+    coordinates: Set[str],
+    varinfo: VarInfoFromDmr,
+) -> Tuple[int, int]:
+    """
+    This method returns the row and column sizes of the coordinate datasets
+
+    """
+    lat_arr, lon_arr = get_lat_lon_arrays(prefetch_dataset, coordinates, varinfo)
+    if lat_arr.ndim > 1:
+        col_size = lat_arr.shape[0]
+        row_size = lat_arr.shape[1]
+    if (lon_arr.shape[0] != lat_arr.shape[0]) or (lon_arr.shape[1] != lat_arr.shape[1]):
+        raise IrregularCoordinateDatasets(lon_arr.shape, lat_arr.shape)
+    if lat_arr.ndim and lon_arr.ndim == 1:
+        col_size = lat_arr.size
+        row_size = lon_arr.size
+    return row_size, col_size
+
+
 def is_latitude_ascending(
     prefetch_dataset: Dataset,
     coordinates: Set[str],

hoss/exceptions.py

@@ -138,6 +138,21 @@ def __init__(self, referring_variable):
         )
 
 
+class IrregularCoordinateDatasets(CustomError):
+    """This exception is raised when HOSS tries to get latitude and longitude
+    datasets and they are missing or empty. These datasets are referred to
+    in the science variables with coordinate attributes.
+
+    """
+
+    def __init__(self, longitude_shape, latitude_shape):
+        super().__init__(
+            'IrregularCoordinateDatasets',
+            f'Longitude dataset shape: "{longitude_shape}"'
+            f'does not match the latitude dataset shape: "{latitude_shape}"',
+        )
+
+
 class UnsupportedShapeFileFormat(CustomError):
     """This exception is raised when the shape file included in the input
     Harmony message is not GeoJSON.

hoss/spatial.py

@@ -80,25 +80,14 @@ def get_spatial_index_ranges(
     around the exterior of the user-defined GeoJSON shape, to ensure the
     correct extents are derived.
 
-    For projected grids which do not follow CF standards, the projected
-    dimension scales are computed based on the values in the coordinate
-    datasets if they are available. The geocorners are obtained from the
-    coordinate datasets and converted to projected meters based on the crs
-    of the product. The dimension scales are then computed based on the
-    grid size and grid resolution
-
     """
     bounding_box = get_harmony_message_bbox(harmony_message)
     index_ranges = {}
-    coordinate_variables = []
 
     geographic_dimensions = varinfo.get_geographic_spatial_dimensions(
         required_variables
     )
     projected_dimensions = varinfo.get_projected_spatial_dimensions(required_variables)
-    if (not geographic_dimensions) and (not projected_dimensions):
-        coordinate_variables = get_coordinate_variables(varinfo, required_variables)
-
     non_spatial_variables = required_variables.difference(
         varinfo.get_spatial_dimensions(required_variables)
     )
@@ -115,7 +104,8 @@ def get_spatial_index_ranges(
                 index_ranges[dimension] = get_geographic_index_range(
                     dimension, varinfo, dimensions_file, bounding_box
                 )
-        if projected_dimensions or coordinate_variables:
+
+        if projected_dimensions:
             for non_spatial_variable in non_spatial_variables:
                 index_ranges.update(
                     get_projected_x_y_index_ranges(
@@ -125,20 +115,35 @@ def get_spatial_index_ranges(
                         index_ranges,
                         bounding_box=bounding_box,
                         shape_file_path=shape_file_path,
-                        coordinates=coordinate_variables,
                     )
                 )
-        return index_ranges
+
+        if (not geographic_dimensions) and (not projected_dimensions):
+            coordinate_variables = get_coordinate_variables(varinfo, required_variables)
+            if coordinate_variables:
+                for non_spatial_variable in non_spatial_variables:
+                    index_ranges.update(
+                        get_x_y_index_ranges_from_coordinates(
+                            non_spatial_variable,
+                            varinfo,
+                            dimensions_file,
+                            coordinate_variables,
+                            index_ranges,
+                            bounding_box=bounding_box,
+                            shape_file_path=shape_file_path,
+                        )
+                    )
+
+    return index_ranges
 
 
 def get_projected_x_y_index_ranges(
     non_spatial_variable: str,
     varinfo: VarInfoFromDmr,
-    dimension_datasets: Dataset,
+    dimensions_file: Dataset,
     index_ranges: IndexRanges,
     bounding_box: BBox = None,
     shape_file_path: str = None,
-    coordinates: Set[str] = set(),
 ) -> IndexRanges:
     """This function returns a dictionary containing the minimum and maximum
     index ranges for a pair of projection x and y coordinates, e.g.:
@@ -157,47 +162,116 @@ def get_projected_x_y_index_ranges(
     projected coordinate points.
 
     """
-    if not coordinates:
-        projected_x, projected_y = get_projected_x_y_variables(
-            varinfo, non_spatial_variable
+    projected_x, projected_y = get_projected_x_y_variables(
+        varinfo, non_spatial_variable
+    )
+
+    if (
+        projected_x is not None
+        and projected_y is not None
+        and not set((projected_x, projected_y)).issubset(set(index_ranges.keys()))
+    ):
+        crs = get_variable_crs(non_spatial_variable, varinfo)
+
+        x_y_extents = get_projected_x_y_extents(
+            dimensions_file[projected_x][:],
+            dimensions_file[projected_y][:],
+            crs,
+            shape_file=shape_file_path,
+            bounding_box=bounding_box,
         )
 
-    else:
-        projected_x = 'projected_x'
-        projected_y = 'projected_y'
-        override_dimension_datasets = update_dimension_variables(
-            dimension_datasets,
-            coordinates,
-            varinfo,
+        x_bounds = get_dimension_bounds(projected_x, varinfo, dimensions_file)
+        y_bounds = get_dimension_bounds(projected_y, varinfo, dimensions_file)
+
+        x_index_ranges = get_dimension_index_range(
+            dimensions_file[projected_x][:],
+            x_y_extents['x_min'],
+            x_y_extents['x_max'],
+            bounds_values=x_bounds,
         )
-        dimension_datasets = override_dimension_datasets
+
+        y_index_ranges = get_dimension_index_range(
+            dimensions_file[projected_y][:],
+            x_y_extents['y_min'],
+            x_y_extents['y_max'],
+            bounds_values=y_bounds,
+        )
+
+        x_y_index_ranges = {projected_x: x_index_ranges, projected_y: y_index_ranges}
+    else:
+        x_y_index_ranges = {}
+
+    return x_y_index_ranges
+
+
+def get_x_y_index_ranges_from_coordinates(
+    non_spatial_variable: str,
+    varinfo: VarInfoFromDmr,
+    dimension_datasets: Dataset,
+    coordinates: Set[str],
+    index_ranges: IndexRanges,
+    bounding_box: BBox = None,
+    shape_file_path: str = None,
+) -> IndexRanges:
+    """This function returns a dictionary containing the minimum and maximum
+    index ranges for a pair of projection x and y coordinates, e.g.:
+
+    index_ranges = {'/x': (20, 42), '/y': (31, 53)}
+
+    This method is called when the CF standards are not followed in the source
+    granule and only coordinate datasets are provided.
+    The coordinate datasets along with the crs is used to calculate the x-y
+    projected dimension scales.
+    The dimensions of the input, non-spatial variable are checked
+    for associated projection x and y coordinates. If these are present,
+    and they have not already been added to the `index_ranges` cache, the
+    extents of the input spatial subset are determined in these projected
+    coordinates. This requires the derivation of a minimum resolution of
+    the target grid in geographic coordinates. Points must be placed along
+    the exterior of the spatial subset shape. All points are then projected
+    from a geographic Coordinate Reference System (CRS) to the target grid
+    CRS. The minimum and maximum values are then derived from these
+    projected coordinate points.
+
+    """
+    crs = get_variable_crs(non_spatial_variable, varinfo)
+
+    projected_x = 'projected_x'
+    projected_y = 'projected_y'
+    override_dimension_datasets = update_dimension_variables(
+        dimension_datasets, coordinates, varinfo, crs
+    )
 
     if (
         projected_x is not None
         and projected_y is not None
         and not set((projected_x, projected_y)).issubset(set(index_ranges.keys()))
     ):
-        crs = get_variable_crs(non_spatial_variable, varinfo)
 
         x_y_extents = get_projected_x_y_extents(
-            dimension_datasets[projected_x][:],
-            dimension_datasets[projected_y][:],
+            override_dimension_datasets[projected_x][:],
+            override_dimension_datasets[projected_y][:],
             crs,
             shape_file=shape_file_path,
             bounding_box=bounding_box,
         )
 
-        x_bounds = get_dimension_bounds(projected_x, varinfo, dimension_datasets)
-        y_bounds = get_dimension_bounds(projected_y, varinfo, dimension_datasets)
+        x_bounds = get_dimension_bounds(
+            projected_x, varinfo, override_dimension_datasets
+        )
+        y_bounds = get_dimension_bounds(
+            projected_y, varinfo, override_dimension_datasets
+        )
 
         x_index_ranges = get_dimension_index_range(
-            dimension_datasets[projected_x][:],
+            override_dimension_datasets[projected_x][:],
             x_y_extents['x_min'],
             x_y_extents['x_max'],
             bounds_values=x_bounds,
         )
         y_index_ranges = get_dimension_index_range(
-            dimension_datasets[projected_y][:],
+            override_dimension_datasets[projected_y][:],
             x_y_extents['y_min'],
             x_y_extents['y_max'],
             bounds_values=y_bounds,

hoss/subset.py

@@ -29,7 +29,6 @@
 from hoss.temporal import get_temporal_index_ranges
 from hoss.utilities import (
     download_url,
-    format_dictionary_string,
     format_variable_set_string,
     get_opendap_nc4,
 )