From 5c97edb4d3a4b71fb605372cb765b246fbe40066 Mon Sep 17 00:00:00 2001 From: joeyschultz <83600443+joeyschultz@users.noreply.github.com> Date: Wed, 30 Oct 2024 13:37:45 -0400 Subject: [PATCH] DAS-2216 - Add quick fixes to support TEMPO level 2 data --- CHANGELOG.md | 24 +++++-- bin/project_local_granule.py | 6 +- docker/service_version.txt | 2 +- swath_projector/earthdata_varinfo_config.json | 43 ++++++++++- swath_projector/interpolation.py | 2 + swath_projector/swath_geometry.py | 16 ++--- swath_projector/utilities.py | 61 +++++++++++++--- tests/unit/test_interpolation.py | 5 ++ tests/unit/test_swath_geometry.py | 12 ++-- tests/unit/test_utilities.py | 72 ++++++++++++++++++- 10 files changed, 202 insertions(+), 41 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3504dd7..c85e6c2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,14 @@ # Changelog +## [v1.2.0] - 2024-10-10 + +### Changed + +- [[DAS-2216](https://bugs.earthdata.nasa.gov/browse/DAS-2216)] + The Swath Projector has been updated with quick fixes to add support for TEMPO level 2 data. These changes include optional transposing of arrays based on dimension sizes, addition of rows_per_scan parameter for ewa interpolation, and updates to the configuration file for TEMPO_O3TOT_L2 to correctly locate coordinate variables and exclude science variables with dimensions that do no match those of the coordinate variables. + ## [v1.1.1] - 2024-09-16 + ### Changed - [[TRT-558](https://bugs.earthdata.nasa.gov/browse/TRT-558)] @@ -12,6 +20,7 @@ also been renamed to `earthdata_varinfo_config.json`. ## [v1.1.0] - 2024-08-29 + ### Changed - [[DAS-1934](https://bugs.earthdata.nasa.gov/browse/DAS-1934)] @@ -37,14 +46,15 @@ include updated documentation and files outlined by the Repository structure changes include: -* Migrating `pymods` directory to `swath_projector`. -* Migrating `swotrepr.py` to `swath_projector/adapter.py`. -* Addition of `swath_projector/main.py`. +- Migrating `pymods` directory to `swath_projector`. +- Migrating `swotrepr.py` to `swath_projector/adapter.py`. +- Addition of `swath_projector/main.py`. For more information on internal releases prior to NASA open-source approval, see legacy-CHANGELOG.md. -[v1.1.1]:(https://github.com/nasa/harmony-swath-projector/releases/tag/1.1.0) -[v1.1.0]:(https://github.com/nasa/harmony-swath-projector/releases/tag/1.0.1) -[v1.0.1]:(https://github.com/nasa/harmony-swath-projector/releases/tag/1.0.1) -[v1.0.0]:(https://github.com/nasa/harmony-swath-projector/releases/tag/1.0.0) +[v1.2.0]: (https://github.com/nasa/harmony-swath-projector/releases/tag/1.2.0) +[v1.1.1]: (https://github.com/nasa/harmony-swath-projector/releases/tag/1.1.1) +[v1.1.0]: (https://github.com/nasa/harmony-swath-projector/releases/tag/1.1.0) +[v1.0.1]: (https://github.com/nasa/harmony-swath-projector/releases/tag/1.0.1) +[v1.0.0]: (https://github.com/nasa/harmony-swath-projector/releases/tag/1.0.0) diff --git a/bin/project_local_granule.py b/bin/project_local_granule.py index 8924d08..4bb6445 100644 --- a/bin/project_local_granule.py +++ b/bin/project_local_granule.py @@ -125,8 +125,8 @@ def project_granule( { 'url': local_file_path, 'temporal': { - 'start': '2021-01-03T23:45:00.000Z', - 'end': '2020-01-04T00:00:00.000Z', + 'start': '2020-01-03T23:45:00.000Z', + 'end': '2025-01-04T00:00:00.000Z', }, 'bbox': [-180, -90, 180, 90], } @@ -141,5 +141,5 @@ def project_granule( reprojector = SwathProjectorAdapter(message, config=config(False)) - with patch('swotrepr.shutil.rmtree', side_effect=rmtree_side_effect): + with patch('swath_projector.adapter.shutil.rmtree', side_effect=rmtree_side_effect): reprojector.invoke() diff --git a/docker/service_version.txt b/docker/service_version.txt index 524cb55..26aaba0 100644 --- a/docker/service_version.txt +++ b/docker/service_version.txt @@ -1 +1 @@ -1.1.1 +1.2.0 diff --git a/swath_projector/earthdata_varinfo_config.json b/swath_projector/earthdata_varinfo_config.json index 66e5e23..7af6c25 100644 --- a/swath_projector/earthdata_varinfo_config.json +++ b/swath_projector/earthdata_varinfo_config.json @@ -1,12 +1,35 @@ { "Identification": "Swath Projector VarInfo configuration", - "Version": 3, + "Version": 4, "CollectionShortNamePath": [ - "ShortName" + "ShortName", + "collection_shortname" ], "Mission": { - "VNP10": "VIIRS" + "VNP10": "VIIRS", + "TEMPO_O3TOT_L2": "TEMPO" }, + "ExcludedScienceVariables": [ + { + "Applicability": { + "Mission": "TEMPO", + "ShortNamePath": "TEMPO_O3TOT_L2" + }, + "VariablePattern": [ + "/support_data/a_priori_layer_o3", + "/support_data/cal_adjustment", + "/support_data/dNdR", + "/support_data/layer_efficiency", + "/support_data/lut_wavelength", + "/support_data/N_value", + "/support_data/N_value_residual", + "/support_data/ozone_sensitivity_ratio", + "/support_data/step_1_N_value_residual", + "/support_data/step_2_N_value_residual", + "/support_data/temp_sensitivity_ratio" + ] + } + ], "MetadataOverrides": [ { "Applicability": { @@ -21,6 +44,20 @@ } ], "_Description": "VNP10 SnowData variables have incorrect relative paths for coordinates." + }, + { + "Applicability": { + "Mission": "TEMPO", + "ShortNamePath": "TEMPO_O3TOT_L2", + "VariablePattern": "^/product/.*|^/support_data/.*|^/geolocation/(?!time$).*" + }, + "Attributes": [ + { + "Name": "coordinates", + "Value": "/geolocation/latitude, /geolocation/longitude" + } + ], + "_Description": "TEMPO_O3TOT_L2 variables only contain basenames for coordinates, which are found in sibling hierarchical groups. This rule fully qualifies the paths to these coordinates. Some variables in these groups are excluded via 'ExcludedScienceVariables'" } ] } diff --git a/swath_projector/interpolation.py b/swath_projector/interpolation.py index 4ab6d7b..9fc92e3 100644 --- a/swath_projector/interpolation.py +++ b/swath_projector/interpolation.py @@ -25,6 +25,7 @@ from swath_projector.utilities import ( create_coordinates_key, get_coordinate_variable, + get_rows_per_scan, get_scale_and_offset, get_variable_file_path, get_variable_numeric_fill_value, @@ -268,6 +269,7 @@ def get_ewa_results( ewa_information['target_area'], variable['values'], maximum_weight_mode=maximum_weight_mode, + rows_per_scan=get_rows_per_scan(variable['values'].shape[0]), ) if variable['fill_value'] is not None: diff --git a/swath_projector/swath_geometry.py b/swath_projector/swath_geometry.py index e6da482..55890cc 100644 --- a/swath_projector/swath_geometry.py +++ b/swath_projector/swath_geometry.py @@ -13,7 +13,7 @@ def get_projected_resolution( - projection: Proj, longitudes: Variable, latitudes: Variable + projection: Proj, longitudes: np.ma.MaskedArray, latitudes: np.ma.MaskedArray ) -> Tuple[float]: """Find the resolution of the target grid in the projected coordinates, x and y. First the perimeter points are found. These are then projected @@ -40,7 +40,7 @@ def get_projected_resolution( def get_extents_from_perimeter( - projection: Proj, longitudes: Variable, latitudes: Variable + projection: Proj, longitudes: np.ma.MaskedArray, latitudes: np.ma.MaskedArray ) -> Tuple[float]: """Find the swath extents in the target CRS. First the perimeter points of unfilled valid pixels are found. These are then projected to the target @@ -59,19 +59,17 @@ def get_extents_from_perimeter( def get_projected_coordinates( coordinates_mask: np.ma.core.MaskedArray, projection: Proj, - longitudes: Variable, - latitudes: Variable, + longitudes: np.ma.MaskedArray, + latitudes: np.ma.MaskedArray, ) -> Tuple[np.ndarray]: """Get the required coordinate points projected in the target Coordinate Reference System (CRS). """ if len(longitudes.shape) == 1: - coordinates = get_all_coordinates(longitudes[:], latitudes[:], coordinates_mask) + coordinates = get_all_coordinates(longitudes, latitudes, coordinates_mask) else: - coordinates = get_perimeter_coordinates( - longitudes[:], latitudes[:], coordinates_mask - ) + coordinates = get_perimeter_coordinates(longitudes, latitudes, coordinates_mask) return reproject_coordinates(coordinates, projection) @@ -135,7 +133,7 @@ def get_absolute_resolution(polygon_area: float, n_pixels: int) -> float: def get_valid_coordinates_mask( - longitudes: Variable, latitudes: Variable + longitudes: np.ma.MaskedArray, latitudes: np.ma.MaskedArray ) -> np.ma.core.MaskedArray: """Get a `numpy` N-d array containing boolean values (0 or 1) indicating whether the elements of both longitude and latitude are valid at that diff --git a/swath_projector/utilities.py b/swath_projector/utilities.py index 35f8ed2..468a14c 100644 --- a/swath_projector/utilities.py +++ b/swath_projector/utilities.py @@ -31,7 +31,8 @@ def get_variable_values( As the variable data are returned as a `numpy.ma.MaskedArray`, the will return no data in the filled pixels. To ensure that the data are correctly handled, the fill value is applied to masked pixels using the - `filled` method. + `filled` method. The variable values are transposed if the `along-track` + dimension size is less than the `across-track` dimension size. """ # TODO: Remove in favour of apply2D or process_subdimension. @@ -42,27 +43,38 @@ def get_variable_values( if len(variable[:].shape) == 1: return make_array_two_dimensional(variable[:]) elif 'time' in input_file.variables and 'time' in variable.dimensions: - # Assumption: Array = (1, y, x) - return variable[0][:].filled(fill_value=fill_value) + # Assumption: Array = (time, along-track, across-track) + return transpose_if_xdim_less_than_ydim(variable[0][:]).filled( + fill_value=fill_value + ) else: - # Assumption: Array = (y, x) - return variable[:].filled(fill_value=fill_value) + # Assumption: Array = (along-track, across-track) + return transpose_if_xdim_less_than_ydim(variable[:]).filled( + fill_value=fill_value + ) def get_coordinate_variable( dataset: Dataset, coordinates_tuple: Tuple[str], coordinate_substring -) -> Optional[Variable]: +) -> Optional[np.ma.MaskedArray]: """Search the coordinate dataset names for a match to the substring, which will be either "lat" or "lon". Return the corresponding variable - from the dataset. Only the base variable name is used, as the group - path may contain either of the strings as part of other words. + data from the dataset. Only the base variable name is used, as the group + path may contain either of the strings as part of other words. The + coordinate variables are transposed if the `along-track`dimension size is + less than the `across-track` dimension size. """ for coordinate in coordinates_tuple: if coordinate_substring in coordinate.split('/')[-1] and variable_in_dataset( coordinate, dataset ): - return dataset[coordinate] + # QuickFix (DAS-2216) for short and wide swaths + if dataset[coordinate].ndim == 1: + return dataset[coordinate][:] + + return transpose_if_xdim_less_than_ydim(dataset[coordinate][:]) + raise MissingCoordinatesError(coordinates_tuple) @@ -216,3 +228,34 @@ def make_array_two_dimensional(one_dimensional_array: np.ndarray) -> np.ndarray: """ return np.expand_dims(one_dimensional_array, 1) + + +def transpose_if_xdim_less_than_ydim( + variable_values: np.ma.MaskedArray, +) -> np.ma.MaskedArray: + """Return transposed variable when variable is wider than tall. + + QuickFix (DAS-2216): We presume that a swath has more rows than columns and + if that's not the case we transpose it so that it does. + """ + if len(variable_values.shape) != 2: + raise ValueError( + f'Input variable must be 2 dimensional, but got {len(variable_values.shape)} dimensions.' + ) + if variable_values.shape[0] < variable_values.shape[1]: + return np.ma.transpose(variable_values).copy() + + return variable_values + + +def get_rows_per_scan(total_rows: int) -> int: + """ + Finds the smallest divisor of the total number of rows. If no divisor is + found, return the total number of rows. + """ + if total_rows < 2: + return 1 + for row_number in range(2, int(total_rows**0.5) + 1): + if total_rows % row_number == 0: + return row_number + return total_rows diff --git a/tests/unit/test_interpolation.py b/tests/unit/test_interpolation.py index 1203af2..5c7a094 100644 --- a/tests/unit/test_interpolation.py +++ b/tests/unit/test_interpolation.py @@ -384,6 +384,7 @@ def test_resample_ewa( self.mock_target_area, mock_values, maximum_weight_mode=False, + rows_per_scan=2, # Added in QuickFix DAS-2216 to be fixed in DAS-2220 ) mock_write_output.assert_called_once_with( self.mock_target_area, @@ -423,6 +424,7 @@ def test_resample_ewa( self.mock_target_area, mock_values, maximum_weight_mode=False, + rows_per_scan=2, # Added in QuickFix DAS-2216 to be fixed in DAS-2220 ) mock_write_output.assert_called_once_with( self.mock_target_area, @@ -491,6 +493,7 @@ def test_resample_ewa_nn( self.mock_target_area, mock_values, maximum_weight_mode=True, + rows_per_scan=2, # Added in QuickFix DAS-2216 to be fixed in DAS-2220 ) mock_write_output.assert_called_once_with( self.mock_target_area, @@ -530,6 +533,7 @@ def test_resample_ewa_nn( self.mock_target_area, mock_values, maximum_weight_mode=True, + rows_per_scan=2, # Added in QuickFix DAS-2216 to be fixed in DAS-2220 ) mock_write_output.assert_called_once_with( self.mock_target_area, @@ -581,6 +585,7 @@ def test_resample_ewa_nn( harmony_target_area, mock_values, maximum_weight_mode=True, + rows_per_scan=2, # Added in QuickFix DAS-2216 to be fixed in DAS-2220 ) # The Harmony target area should be given to the output function diff --git a/tests/unit/test_swath_geometry.py b/tests/unit/test_swath_geometry.py index 68bae25..896b380 100644 --- a/tests/unit/test_swath_geometry.py +++ b/tests/unit/test_swath_geometry.py @@ -65,8 +65,8 @@ def setUpClass(cls): def setUp(self): self.test_dataset = Dataset(self.test_path) - self.longitudes = self.test_dataset['lon'] - self.latitudes = self.test_dataset['lat'] + self.longitudes = self.test_dataset['lon'][:] + self.latitudes = self.test_dataset['lat'][:] def tearDown(self): self.test_dataset.close() @@ -101,8 +101,8 @@ def test_get_projected_resolution_1d(self): """Ensure the calculated one-dimensional resolution is correct.""" resolution = get_projected_resolution( self.geographic_projection, - self.test_dataset['lon_1d'], - self.test_dataset['lat_1d'], + self.test_dataset['lon_1d'][:], + self.test_dataset['lat_1d'][:], ) self.assertAlmostEqual(resolution, 5.0) @@ -167,9 +167,7 @@ def test_get_perimeter_coordinates(self): np.logical_not(valid_pixels), np.ones(self.longitudes.shape) ) - coordinates = get_perimeter_coordinates( - self.longitudes[:], self.latitudes[:], mask - ) + coordinates = get_perimeter_coordinates(self.longitudes, self.latitudes, mask) self.assertCountEqual(coordinates, expected_points) diff --git a/tests/unit/test_utilities.py b/tests/unit/test_utilities.py index abfcb1b..fe3eae3 100644 --- a/tests/unit/test_utilities.py +++ b/tests/unit/test_utilities.py @@ -10,12 +10,14 @@ construct_absolute_path, create_coordinates_key, get_coordinate_variable, + get_rows_per_scan, get_scale_and_offset, get_variable_file_path, get_variable_numeric_fill_value, get_variable_values, make_array_two_dimensional, qualify_reference, + transpose_if_xdim_less_than_ydim, variable_in_dataset, ) @@ -79,7 +81,9 @@ def test_get_variable_values(self): wind_speed_values = get_variable_values(dataset, wind_speed, None) self.assertIsInstance(wind_speed_values, np.ndarray) self.assertEqual(len(wind_speed_values.shape), 2) - self.assertEqual(wind_speed_values.shape, wind_speed.shape) + self.assertEqual( + wind_speed_values.shape, np.ma.transpose(wind_speed).shape + ) with self.subTest('Masked values are set to fill value.'): fill_value = 210 @@ -135,7 +139,7 @@ def test_get_coordinate_variables(self): dataset, coordinates_tuple, coordinate ) - self.assertIsInstance(coordinates, Variable) + self.assertIsInstance(coordinates, np.ma.MaskedArray) with self.subTest( 'Non existent coordinate variable "latitude" returns MissingCoordinatesError' @@ -358,3 +362,67 @@ def test_make_array_two_dimensional(self): self.assertEqual(len(output_array.shape), 2) np.testing.assert_array_equal(output_array, expected_output) + + +class TestTransposeIfXdimLessThanYdim(TestCase): + + def test_wider_than_tall(self): + """Test case where x dim < y dim and should transpose.""" + input_array = np.ma.array([[1, 2, 3], [4, 5, 6]]) + expected_output = np.ma.array([[1, 4], [2, 5], [3, 6]]) + result = transpose_if_xdim_less_than_ydim(input_array) + np.testing.assert_array_equal(result, expected_output) + self.assertEqual(result.shape, (3, 2)) + + def test_taller_than_wide(self): + """Test case where x < y and should not transpose.""" + input_array = np.ma.array([[1, 2], [3, 4], [5, 6]]) + result = transpose_if_xdim_less_than_ydim(input_array) + np.testing.assert_array_equal(result, input_array) + self.assertEqual(result.shape, (3, 2)) + + def test_square_array(self): + """Test case where y dim == x dim and should not transpose.""" + input_array = np.ma.array([[1, 2], [3, 4]]) + result = transpose_if_xdim_less_than_ydim(input_array) + np.testing.assert_array_equal(result, input_array) + self.assertEqual(result.shape, (2, 2)) + + def test_1d_array(self): + """Test case with a 1D array""" + input_array = np.ma.array([1, 2, 3]) + with self.assertRaisesRegex(ValueError, 'variable must be 2 dimensional'): + transpose_if_xdim_less_than_ydim(input_array) + + def test_3d_array(self): + """Test case with a 3D array""" + input_array = np.ma.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]]) + with self.assertRaisesRegex(ValueError, 'variable must be 2 dimensional'): + transpose_if_xdim_less_than_ydim(input_array) + + def test_masked_array(self): + """Test case with a masked array""" + input_array = np.ma.array( + [[1, 2, 3], [4, 5, 6]], mask=[[True, False, False], [False, True, False]] + ) + expected_output = np.ma.array( + [[1, 4], [2, 5], [3, 6]], + mask=[[True, False], [False, True], [False, False]], + ) + result = transpose_if_xdim_less_than_ydim(input_array) + np.testing.assert_array_equal(result, expected_output) + np.testing.assert_array_equal(result.mask, expected_output.mask) + + +class TestGetRowsPerScan(TestCase): + def test_number_less_than_2(self): + self.assertEqual(get_rows_per_scan(1), 1) + + def test_even_composite_number(self): + self.assertEqual(get_rows_per_scan(4), 2) + + def test_odd_composite_number(self): + self.assertEqual(get_rows_per_scan(9), 3) + + def test_prime_number(self): + self.assertEqual(get_rows_per_scan(3), 3)