diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c1720d..66f7ffa 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/bin/extract-release-notes.sh b/bin/extract-release-notes.sh index 3da532c..bb15f23 100755 --- a/bin/extract-release-notes.sh +++ b/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" diff --git a/docker/service_version.txt b/docker/service_version.txt index 6d7de6e..9084fa2 100644 --- a/docker/service_version.txt +++ b/docker/service_version.txt @@ -1 +1 @@ -1.0.2 +1.1.0 diff --git a/docs/HyBIG-Example-Usage.ipynb b/docs/HyBIG-Example-Usage.ipynb index a72fa04..ecf34fd 100644 --- a/docs/HyBIG-Example-Usage.ipynb +++ b/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", diff --git a/harmony_browse_image_generator/sizes.py b/harmony_browse_image_generator/sizes.py index d7c9f25..f6e9473 100644 --- a/harmony_browse_image_generator/sizes.py +++ b/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: diff --git a/pip_requirements.txt b/pip_requirements.txt index f0352af..0353a86 100644 --- a/pip_requirements.txt +++ b/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 diff --git a/tests/fixtures/RGB.byte.small.tif b/tests/fixtures/RGB.byte.small.tif new file mode 100644 index 0000000..21f3b51 Binary files /dev/null and b/tests/fixtures/RGB.byte.small.tif differ diff --git a/tests/test_adapter.py b/tests/test_adapter.py index 838657c..454e23d 100644 --- a/tests/test_adapter.py +++ b/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, ) diff --git a/tests/unit/test_browse.py b/tests/unit/test_browse.py index f9bf619..52ad50f 100644 --- a/tests/unit/test_browse.py +++ b/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( [ [ diff --git a/tests/unit/test_sizes.py b/tests/unit/test_sizes.py index d309ced..4899850 100644 --- a/tests/unit/test_sizes.py +++ b/tests/unit/test_sizes.py @@ -1,5 +1,6 @@ """Tests covering the size module.""" +from pathlib import Path from unittest import TestCase from unittest.mock import MagicMock, patch @@ -13,7 +14,7 @@ from harmony_browse_image_generator.exceptions import HyBIGValueError from harmony_browse_image_generator.sizes import ( METERS_PER_DEGREE, - best_guess_scale_extent, + ScaleExtent, best_guess_target_dimensions, choose_scale_extent, choose_target_dimensions, @@ -27,7 +28,6 @@ get_cells_per_tile, get_rasterio_parameters, get_target_grid_parameters, - icd_defined_extent_from_crs, needs_tiling, resolution_in_target_crs_units, ) @@ -120,39 +120,44 @@ def test_grid_parameters_from_harmony_message_has_complete_information(self): self.assertDictEqual(expected_parameters, actual_parameters) def test_grid_parameters_from_harmony_no_message_information(self): - """input granule is in preferred_crs on a 25km grid - - But has grid extents different from the ICD defaults. - - """ + """input granule is in preferred_crs on a 25km grid""" + crs = CRS.from_epsg(sp_seaice_grid['epsg']) + height = sp_seaice_grid['height'] + width = sp_seaice_grid['width'] + img_transform = Affine.translation( + sp_seaice_grid['left'], sp_seaice_grid['top'] + ) * Affine.scale(sp_seaice_grid['xres'], -1 * sp_seaice_grid['yres']) - preferred_scale_extent = { - 'xmin': -4194304.0, - 'ymin': -4194304.0, - 'xmax': 4194304.0, - 'ymax': 4194304.0, + left, bottom, right, top = rasterio.transform.array_bounds( + height, width, img_transform + ) + image_scale_extent = { + 'xmin': left, + 'ymin': bottom, + 'xmax': right, + 'ymax': top, } + expected_height = round( - (preferred_scale_extent['ymax'] - preferred_scale_extent['ymin']) + (image_scale_extent['ymax'] - image_scale_extent['ymin']) / sp_seaice_grid['yres'] - ) # 336 + ) # 332 expected_width = round( - (preferred_scale_extent['xmax'] - preferred_scale_extent['xmin']) + (image_scale_extent['xmax'] - image_scale_extent['xmin']) / sp_seaice_grid['xres'] - ) # 336 + ) # 316 expected_y_resolution = ( - preferred_scale_extent['ymax'] - preferred_scale_extent['ymin'] + image_scale_extent['ymax'] - image_scale_extent['ymin'] ) / expected_height expected_x_resolution = ( - preferred_scale_extent['xmax'] - preferred_scale_extent['xmin'] + image_scale_extent['xmax'] - image_scale_extent['xmin'] ) / expected_width expected_transform = Affine.translation( - preferred_scale_extent['xmin'], preferred_scale_extent['ymax'] + image_scale_extent['xmin'], image_scale_extent['ymax'] ) * Affine.scale(expected_x_resolution, -1 * expected_y_resolution) - crs = CRS.from_epsg(sp_seaice_grid['epsg']) expected_parameters = { 'height': expected_height, 'width': expected_width, @@ -160,14 +165,11 @@ def test_grid_parameters_from_harmony_no_message_information(self): 'transform': expected_transform, } - height = sp_seaice_grid['height'] - width = sp_seaice_grid['width'] with rasterio_test_file( height=height, width=width, crs=crs, - transform=Affine.translation(sp_seaice_grid['left'], sp_seaice_grid['top']) - * Affine.scale(sp_seaice_grid['xres'], -1 * sp_seaice_grid['yres']), + transform=img_transform, ) as tmp_file: message = Message({'format': {}}) with open_rasterio(tmp_file) as rio_data_array: @@ -397,6 +399,8 @@ def test_create_tile_output_parameters( class TestChooseScaleExtent(TestCase): """Test for correct scale extents.""" + fixtures = Path(__file__).resolve().parent.parent / 'fixtures' + def test_scale_extent_in_harmony_message(self): """Basic case of user supplied scaleExtent.""" message = Message( @@ -416,30 +420,43 @@ def test_scale_extent_in_harmony_message(self): 'ymax': 500.0, } crs = None - actual_scale_extent = choose_scale_extent(message, crs) + actual_scale_extent = choose_scale_extent(message, crs, None) self.assertDictEqual(expected_scale_extent, actual_scale_extent) - @patch('harmony_browse_image_generator.sizes.best_guess_scale_extent') - @patch('harmony_browse_image_generator.sizes.icd_defined_extent_from_crs') - def test_preferred_crs(self, mock_icd_extents, mock_best_guess): - """Preferred CRS will default to ICD expected.""" - target_crs = CRS.from_string(PREFERRED_CRS['north']) - message = Message({}) - choose_scale_extent(message, target_crs) - - mock_icd_extents.assert_called_once_with(target_crs) - mock_best_guess.assert_not_called() - - @patch('harmony_browse_image_generator.sizes.best_guess_scale_extent') - @patch('harmony_browse_image_generator.sizes.icd_defined_extent_from_crs') - def test_guess_extents_from_crs(self, mock_icd_extents, mock_best_guess): - """un-preferred CRS will try to guess scale extents.""" - target_crs = CRS.from_string('EPSG:6931') - message = Message({}) - choose_scale_extent(message, target_crs) - - mock_icd_extents.assert_not_called() - mock_best_guess.assert_called_once_with(target_crs) + def test_scale_extent_from_input_image_and_no_crs_transformation(self): + """Ensure no change of output extent when src_crs == target_crs""" + + with open_rasterio( + self.fixtures / 'RGB.byte.small.tif', mode='r', mask_and_scale=True + ) as in_array: + source_crs = in_array.rio.crs + left, bottom, right, top = in_array.rio.bounds() + expected_scale_extent = ScaleExtent( + {'xmin': left, 'ymin': bottom, 'xmax': right, 'ymax': top} + ) + + actual_scale_extent = choose_scale_extent({}, source_crs, in_array) + self.assertEqual(actual_scale_extent, expected_scale_extent) + + def test_scale_extent_from_input_image_with_crs_transformation(self): + """Ensure no change of output extent when src_crs == target_crs""" + target_crs = CRS.from_string(PREFERRED_CRS['global']) + with open_rasterio( + self.fixtures / 'RGB.byte.small.tif', mode='r', mask_and_scale=True + ) as in_array: + + left, bottom, right, top = ( + -78.95864996539397, + 23.568866283727235, + -76.59780097339339, + 25.550618627487918, + ) + expected_scale_extent = ScaleExtent( + {'xmin': left, 'ymin': bottom, 'xmax': right, 'ymax': top} + ) + + actual_scale_extent = choose_scale_extent({}, target_crs, in_array) + self.assertEqual(actual_scale_extent, expected_scale_extent) class TestChooseTargetDimensions(TestCase): @@ -492,97 +509,6 @@ def test_message_has_just_one_dimension(self, mock_best_guess_target_dimensions) ) -class TestICDDefinedExtentFromCrs(TestCase): - """Ensure standard CRS yield standard scaleExtents. - - The expected scale extents are pulled from the ICD document. - - """ - - def test_global_preferred_crs(self): - crs = CRS.from_string(PREFERRED_CRS['global']) - expected_scale_extent = { - 'xmin': -180.0, - 'ymin': -90.0, - 'xmax': 180.0, - 'ymax': 90.0, - } - actual_scale_extent = icd_defined_extent_from_crs(crs) - self.assertDictEqual(expected_scale_extent, actual_scale_extent) - - def test_south_preferred_crs(self): - crs = CRS.from_string(PREFERRED_CRS['south']) - expected_scale_extent = { - 'xmin': -4194304.0, - 'ymin': -4194304.0, - 'xmax': 4194304.0, - 'ymax': 4194304.0, - } - actual_scale_extent = icd_defined_extent_from_crs(crs) - self.assertDictEqual(expected_scale_extent, actual_scale_extent) - - def test_north_preferred_crs(self): - crs = CRS.from_string(PREFERRED_CRS['north']) - expected_scale_extent = { - 'xmin': -4194304.0, - 'ymin': -4194304.0, - 'xmax': 4194304.0, - 'ymax': 4194304.0, - } - actual_scale_extent = icd_defined_extent_from_crs(crs) - self.assertDictEqual(expected_scale_extent, actual_scale_extent) - - def test_non_preferred(self): - crs = CRS.from_string('EPSG:4311') - with self.assertRaises(HyBIGValueError) as excepted: - icd_defined_extent_from_crs(crs) - - self.assertEqual(excepted.exception.message, 'Invalid input CRS: EPSG:4311') - pass - - -class TestBestGuessScaleExtent(TestCase): - def test_unprojected_crs_has_area_of_use(self): - crs = CRS.from_epsg(4326) - expected_extent = {'xmin': -180.0, 'ymin': -90.0, 'xmax': 180.0, 'ymax': 90.0} - actual_extent = best_guess_scale_extent(crs) - self.assertDictEqual(expected_extent, actual_extent) - - def test_projected_crs_has_area_of_use(self): - """North Pole LAEA Canada""" - crs = CRS.from_epsg(3573) - expected_extent = { - 'xmin': -4859208.992805643, - 'ymin': -4886873.230107171, - 'xmax': 4859208.992805643, - 'ymax': 4886873.230107171, - } - - actual_extent = best_guess_scale_extent(crs) - self.assertDictEqual(expected_extent, actual_extent) - - def test_unprojected_crs_has_no_area_of_use(self): - """Use EPSG 4269 created with a dictionary.""" - crs = CRS.from_dict( - {'proj': 'longlat', 'datum': 'NAD83', 'no_defs': None, 'type': 'crs'} - ) - - expected_extent = {'xmin': -180.0, 'ymin': -90.0, 'xmax': 180.0, 'ymax': 90.0} - - actual_extent = best_guess_scale_extent(crs) - - self.assertDictEqual(expected_extent, actual_extent) - - def test_unprojected_nonstandard_crs_has_no_area_of_use(self): - """Use EPSG 4269.""" - crs = CRS.from_epsg(4269) - - # west=167.65, south=14.92, east=-40.73, north=86.45 - expected_extent = {'xmin': 167.65, 'ymin': 14.92, 'xmax': -40.73, 'ymax': 86.45} - actual_extent = best_guess_scale_extent(crs) - self.assertDictEqual(expected_extent, actual_extent) - - class TestBestGuessTargetDimensions(TestCase): def test_projected_crs(self): """A coarse resolution image uses the input granules' height and width