DAS-2115: Change how the default scale_extent is computed (#9)

CHANGELOG.md

@@ -1,17 +1,24 @@
-## v1.0.2
-### 2024-04-05
+## v1.1.0 - 2024-04-30
+
+### Changed
+Changes the computation for an output image's default scale extent. Previously
+we considered ICD preferred ScaleExtents as covering the entire globe or pole.
+This change now takes just the input image bounds and converts them to the target crs
+and uses that transformed boundry as the default region to make a scale extent from.
+
+Upgraded harmony-service-lib to v1.0.26
+
+## v1.0.2 - 2024-04-05
 
 This version of HyBIG correctly handles missing/bad input data marked by _FillValue or NoData.
 Anytime a bad value occurs in the input raster, the output png image will set to transparent.
 
-## v1.0.1
-### 2024-04-05
+## v1.0.1 - 2024-04-05
 
 This version of HyBIG updates the repository to use `black` code formatting
 throughout. There should be no functional change to the service.
 
-## v1.0.0
-### 2024-01-22
+## v1.0.0 - 2024-01-22
 This version of the Harmony Browse Image Generator (HyBIG) contains all
 functionality previously released internally to EOSDIS as
 sds/harmony-browse-image-generator:0.0.11.

bin/extract-release-notes.sh

@@ -19,4 +19,4 @@ VERSION_PATTERN="^## v"
 result=$(awk "/$VERSION_PATTERN/{c++; if(c==2) exit;} c==1" "$CHANGELOG_FILE")
 
 # Print the result
-echo "$result" |  grep -v "^#"
+echo "$result" |  grep -v "$VERSION_PATTERN"

docker/service_version.txt

@@ -1 +1 @@
-1.0.2
+1.1.0

docs/HyBIG-Example-Usage.ipynb

@@ -233,6 +233,7 @@
     "    collection=aster_collection,\n",
     "    granule_id=aster_granule,\n",
     "    scale_extent=scale_extent,\n",
+    "    crs='EPSG:4326',\n",
     "    format='image/jpeg',\n",
     ")\n",
     "\n",
@@ -312,6 +313,7 @@
     "    collection=measures_collection,\n",
     "    granule_id=measures_granule,\n",
     "    scale_size=scale_sizes,\n",
+    "    crs='EPSG:4326',\n",
     "    format='image/png',\n",
     ")\n",
     "\n",
@@ -435,6 +437,7 @@
     "    granule_id=measures_granule,\n",
     "    scale_extent=iceland_extent,\n",
     "    scale_size=iceland_scale_size,\n",
+    "    crs='EPSG:4326',\n",
     "    format='image/png',\n",
     ")\n",
     "\n",

harmony_browse_image_generator/sizes.py

@@ -125,19 +125,22 @@ def get_target_grid_parameters(message: Message, data_array: DataArray) -> GridP
 
     """
     target_crs = choose_target_crs(message.format.srs, data_array)
-    target_scale_extent = choose_scale_extent(message, target_crs)
+    target_scale_extent = choose_scale_extent(message, target_crs, data_array)
     target_dimensions = choose_target_dimensions(
         message, data_array, target_scale_extent, target_crs
     )
     return get_rasterio_parameters(target_crs, target_scale_extent, target_dimensions)
 
 
-def choose_scale_extent(message: Message, target_crs: CRS) -> ScaleExtent:
-    """Return the scaleExtent of the target image.
+def choose_scale_extent(
+    message: Message, target_crs: CRS, data_array: DataArray
+) -> ScaleExtent:
+    """Return the scaleExtent for the target image.
 
-    Check the message for a defined scale extent and returns that or returns
-    the best alternative based on the target_crs either returning it from a
-    lookup based on the ICD, or computed with pyproj.area_of_use
+    Returns a scale extent found in the input Message.
+
+    Otherwise, computes a bounding box in the target CRS based on the input
+    granule extent.
 
     """
     if has_scale_extents(message):
@@ -150,11 +153,11 @@ def choose_scale_extent(message: Message, target_crs: CRS) -> ScaleExtent:
                 'ymax': message.format.scaleExtent.y.max,
             }
         )
-    elif is_preferred_crs(target_crs):
-        scale_extent = icd_defined_extent_from_crs(target_crs)
     else:
-        # compute a best guess area based on the CRS's region of interest
-        scale_extent = best_guess_scale_extent(target_crs)
+        left, bottom, right, top = data_array.rio.transform_bounds(target_crs)
+        scale_extent = ScaleExtent(
+            {'xmin': left, 'ymin': bottom, 'xmax': right, 'ymax': top}
+        )
     return scale_extent
 
 
@@ -322,68 +325,6 @@ def needs_tiling(grid_parameters: GridParams) -> bool:
     return grid_parameters['height'] * grid_parameters['width'] > MAX_UNTILED_GRIDCELLS
 
 
-def icd_defined_extent_from_crs(crs: CRS) -> ScaleExtent:
-    """return the predefined scaleExtent for a GIBS image.
-
-    looks up which projetion is being used and returns the scaleExtent.
-
-    """
-    if crs.to_string() == PREFERRED_CRS['global']:
-        scale_extent = ScaleExtent(
-            {'xmin': -180.0, 'ymin': -90.0, 'xmax': 180.0, 'ymax': 90.0}
-        )
-    elif crs.to_string() in [PREFERRED_CRS['north'], PREFERRED_CRS['south']]:
-        # both north and south preferred CRSs have same extents.
-        scale_extent = ScaleExtent(
-            {
-                'xmin': -4194304.0,
-                'ymin': -4194304.0,
-                'xmax': 4194304.0,
-                'ymax': 4194304.0,
-            }
-        )
-    else:
-        raise HyBIGValueError(f'Invalid input CRS: {crs.to_string()}')
-
-    return scale_extent
-
-
-def best_guess_scale_extent(in_crs: CRS) -> ScaleExtent:
-    """Guess the best scale extent.
-
-    This routine will try to guess what a user intended if they did not include
-    a scaleExtent and also used a non-preferred CRS. We convert the CRS into a
-    pyproj crs check for an area of use.  If this exists we return the bounds,
-    projecting them if the crs is a projected crs.
-
-    if no area_of_use exists, we return the ICD defined scale extent that
-    relates to the closest prefered CRS.
-
-    """
-    crs = pyCRS(in_crs.to_wkt(version='WKT2'))
-    if crs.area_of_use is None:
-        best_crs = choose_best_crs_from_metadata(crs)
-        scale_extent = icd_defined_extent_from_crs(best_crs)
-    elif crs.is_projected:
-        transformer = Transformer.from_crs(crs.geodetic_crs, crs, always_xy=True)
-        projected_bounds = transformer.transform_bounds(*crs.area_of_use.bounds)
-        scale_extent = {
-            'xmin': projected_bounds[0],
-            'ymin': projected_bounds[1],
-            'xmax': projected_bounds[2],
-            'ymax': projected_bounds[3],
-        }
-    else:
-        scale_extent = {
-            'xmin': crs.area_of_use.bounds[0],
-            'ymin': crs.area_of_use.bounds[1],
-            'xmax': crs.area_of_use.bounds[2],
-            'ymax': crs.area_of_use.bounds[3],
-        }
-
-    return scale_extent
-
-
 def best_guess_target_dimensions(
     data_array: DataArray, scale_extent: ScaleExtent, target_crs: CRS
 ) -> Dimensions:

pip_requirements.txt

@@ -1,4 +1,4 @@
-harmony-service-lib~=1.0.25
+harmony-service-lib~=1.0.26
 pystac~=0.5.6
 matplotlib==3.7.1
 rasterio==1.3.6

tests/test_adapter.py

@@ -10,7 +10,7 @@
 from harmony.message import Message
 from harmony.util import config
 from pystac import Catalog
-from rasterio.transform import from_bounds
+from rasterio.transform import array_bounds, from_bounds
 from rasterio.warp import Resampling
 from rioxarray import open_rasterio
 
@@ -195,14 +195,20 @@ def move_tif(*args, **kwargs):
         # Scale Extent from icd
         # dimensions from input data
         rio_data_array = open_rasterio(self.red_tif_fixture)
-        icd_scale_extent = {'xmin': -180.0, 'ymin': -90.0, 'xmax': 180.0, 'ymax': 90.0}
-        expected_width = round((180 - -180) / rio_data_array.rio.transform().a)
-        expected_height = round((90 - -90) / -rio_data_array.rio.transform().e)
+        left, bottom, right, top = array_bounds(
+            rio_data_array.rio.width,
+            rio_data_array.rio.height,
+            rio_data_array.rio.transform(),
+        )
+        image_scale_extent = {'xmin': left, 'ymin': bottom, 'xmax': right, 'ymax': top}
+
+        expected_width = round((right - left) / rio_data_array.rio.transform().a)
+        expected_height = round((top - bottom) / -rio_data_array.rio.transform().e)
         expected_transform = from_bounds(
-            icd_scale_extent['xmin'],
-            icd_scale_extent['ymin'],
-            icd_scale_extent['xmax'],
-            icd_scale_extent['ymax'],
+            image_scale_extent['xmin'],
+            image_scale_extent['ymin'],
+            image_scale_extent['xmax'],
+            image_scale_extent['ymax'],
             expected_width,
             expected_height,
         )
@@ -424,15 +430,20 @@ def move_tif(*args, **kwargs):
         # Scale Extent from icd
         # dimensions from input data
         rio_data_array = open_rasterio(self.red_tif_fixture)
-        icd_scale_extent = {'xmin': -180.0, 'ymin': -90.0, 'xmax': 180.0, 'ymax': 90.0}
+        left, bottom, right, top = array_bounds(
+            rio_data_array.rio.width,
+            rio_data_array.rio.height,
+            rio_data_array.rio.transform(),
+        )
+        image_scale_extent = {'xmin': left, 'ymin': bottom, 'xmax': right, 'ymax': top}
 
-        expected_width = round((180 - -180) / rio_data_array.rio.transform().a)
-        expected_height = round((90 - -90) / -rio_data_array.rio.transform().e)
+        expected_width = round((right - left) / rio_data_array.rio.transform().a)
+        expected_height = round((top - bottom) / -rio_data_array.rio.transform().e)
         expected_transform = from_bounds(
-            icd_scale_extent['xmin'],
-            icd_scale_extent['ymin'],
-            icd_scale_extent['xmax'],
-            icd_scale_extent['ymax'],
+            image_scale_extent['xmin'],
+            image_scale_extent['ymin'],
+            image_scale_extent['xmax'],
+            image_scale_extent['ymax'],
             expected_width,
             expected_height,
         )

tests/unit/test_browse.py

@@ -16,6 +16,7 @@
 from rasterio import Affine
 from rasterio.crs import CRS
 from rasterio.io import DatasetReader, DatasetWriter
+from rasterio.transform import array_bounds
 from rasterio.warp import Resampling
 from xarray import DataArray
 
@@ -167,6 +168,8 @@ def test_create_browse_imagery_with_mocks(
         da_mock.rio.count = 1
         in_dataset_mock.colormap = Mock(side_effect=ValueError)
 
+        da_mock.rio.transform_bounds.return_value = array_bounds(4, 4, file_transform)
+
         expected_raster = np.array(
             [
                 [