From 36e15c7ba107e815461773711916a1f8bc5cd504 Mon Sep 17 00:00:00 2001 From: sudhamurthy Date: Fri, 4 Oct 2024 11:46:36 -0400 Subject: [PATCH] DAS-2232 - moved new methods to a separate file --- CHANGELOG.md | 5 +- hoss/coordinate_utilities.py | 82 +++++++++---------- tests/unit/test_spatial.py | 149 +++++++++++++++++++++++++++++++++++ 3 files changed, 188 insertions(+), 48 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eb8b4b4..a10326a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,9 +1,8 @@ ## v1.1.0 ### 2024-09-10 -This version of HOSS provides support for products that do not comply with CF standards like SMAP L3. -The SMAP L3 products do not have grid mapping variable to provide the projection information. -They also do not have x-y dimension scales for the projected grid. +This version of HOSS provides support for gridded products that do not contain CF-Convention +compliant grid mapping variables and 1-D dimension variables, such as SMAP L3. New methods are added to retrieve dimension scales from coordinate attributes and grid mappings, using overrides specified in the hoss_config.json configuration file. - `get_coordinate_variables` gets coordinate datasets when the dimension scales are not present in diff --git a/hoss/coordinate_utilities.py b/hoss/coordinate_utilities.py index f1cd0fd..1eb53c4 100644 --- a/hoss/coordinate_utilities.py +++ b/hoss/coordinate_utilities.py @@ -8,7 +8,6 @@ import numpy as np from netCDF4 import Dataset from numpy import ndarray -from numpy.ma.core import MaskedArray from pyproj import CRS from varinfo import VariableFromDmr, VarInfoFromDmr @@ -133,18 +132,30 @@ def update_dimension_variables( (5) Generate the x-y dimscale array and return to the calling method """ - # if there is more than one grid, determine the lat/lon coordinate - # associated with this variable - row_size, col_size = get_row_col_sizes_from_coordinate_datasets( + lat_arr, lon_arr = get_lat_lon_arrays( prefetch_dataset, latitude_coordinate, longitude_coordinate, ) + if not lat_arr.size: + raise MissingCoordinateDataset('latitude') + if not lon_arr.size: + raise MissingCoordinateDataset('longitude') + + lat_fill, lon_fill = get_fill_values_for_coordinates( + latitude_coordinate, longitude_coordinate + ) + + row_size, col_size = get_row_col_sizes_from_coordinate_datasets( + lat_arr, + lon_arr, + ) geo_grid_corners = get_geo_grid_corners( - prefetch_dataset, - latitude_coordinate, - longitude_coordinate, + lat_arr, + lon_arr, + lat_fill, + lon_fill, ) x_y_extents = get_x_y_extents_from_geographic_points(geo_grid_corners, crs) @@ -158,11 +169,7 @@ def update_dimension_variables( y_resolution = (y_max - y_min) / col_size # create the xy dim scales - lat_asc, lon_asc = is_lat_lon_ascending( - prefetch_dataset, - latitude_coordinate, - longitude_coordinate, - ) + lat_asc, lon_asc = is_lat_lon_ascending(lat_arr, lon_arr, lat_fill, lon_fill) if lon_asc: x_dim = np.arange(x_min, x_max, x_resolution) @@ -178,17 +185,14 @@ def update_dimension_variables( def get_row_col_sizes_from_coordinate_datasets( - prefetch_dataset: Dataset, - latitude_coordinate: VariableFromDmr, - longitude_coordinate: VariableFromDmr, + lat_arr: ndarray, + lon_arr: ndarray, ) -> 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, latitude_coordinate, longitude_coordinate - ) + if lat_arr.ndim > 1: col_size = lat_arr.shape[0] row_size = lat_arr.shape[1] @@ -201,26 +205,27 @@ def get_row_col_sizes_from_coordinate_datasets( def is_lat_lon_ascending( - prefetch_dataset: Dataset, - latitude_coordinate: VariableFromDmr, - longitude_coordinate: VariableFromDmr, + lat_arr: ndarray, + lon_arr: ndarray, + lat_fill: float, + lon_fill: float, ) -> tuple[bool, bool]: """ Checks if the latitude and longitude cooordinate datasets have values that are ascending """ - lat_arr, lon_arr = get_lat_lon_arrays( - prefetch_dataset, latitude_coordinate, longitude_coordinate - ) + lat_col = lat_arr[:, 0] lon_row = lon_arr[0, :] - first_index, last_index = np.ma.flatnotmasked_edges(lat_col) - latitude_ascending = lat_col.size == 1 or lat_col[first_index] < lat_col[last_index] + lat_col_valid_indices = get_valid_indices(lon_row, lat_fill, 'latitude') + latitude_ascending = ( + lat_col[lat_col_valid_indices[1]] > lat_col[lat_col_valid_indices[0]] + ) - first_index, last_index = np.ma.flatnotmasked_edges(lat_col) + lon_row_valid_indices = get_valid_indices(lon_row, lon_fill, 'longitude') longitude_ascending = ( - lon_row.size == 1 or lon_row[first_index] < lon_row[last_index] + lon_row[lon_row_valid_indices[1]] > lon_row[lon_row_valid_indices[0]] ) return latitude_ascending, longitude_ascending @@ -245,9 +250,10 @@ def get_lat_lon_arrays( def get_geo_grid_corners( - prefetch_dataset: Dataset, - latitude_coordinate: VariableFromDmr, - longitude_coordinate: VariableFromDmr, + lat_arr: ndarray, + lon_arr: ndarray, + lat_fill: float, + lon_fill: float, ) -> list[Tuple[float]]: """ This method is used to return the lat lon corners from a 2D @@ -258,24 +264,10 @@ def get_geo_grid_corners( still needs to be addressed. It will raise an exception in those cases. """ - lat_arr, lon_arr = get_lat_lon_arrays( - prefetch_dataset, - latitude_coordinate, - longitude_coordinate, - ) - - if not lat_arr.size: - raise MissingCoordinateDataset('latitude') - if not lon_arr.size: - raise MissingCoordinateDataset('longitude') top_left_row_idx = 0 top_left_col_idx = 0 - lat_fill, lon_fill = get_fill_values_for_coordinates( - latitude_coordinate, longitude_coordinate - ) - # get the first row from the longitude dataset lon_row = lon_arr[top_left_row_idx, :] lon_row_valid_indices = get_valid_indices(lon_row, lon_fill, 'longitude') diff --git a/tests/unit/test_spatial.py b/tests/unit/test_spatial.py index e5ee6d3..e0c72da 100644 --- a/tests/unit/test_spatial.py +++ b/tests/unit/test_spatial.py @@ -17,6 +17,7 @@ get_longitude_in_grid, get_projected_x_y_index_ranges, get_spatial_index_ranges, + get_x_y_index_ranges_from_coordinates, ) @@ -182,6 +183,154 @@ def test_get_spatial_index_ranges_geographic(self): {'/latitude': (5, 44), '/longitude': (160, 199)}, ) + @patch('hoss.spatial.get_dimension_index_range') + @patch('hoss.spatial.get_projected_x_y_extents') + def test_get_x_y_index_ranges_from_coordinates( + self, mock_get_x_y_extents, mock_get_dimension_index_range + ): + """Ensure that x and y index ranges are only requested only when there are + no projected dimensions and when there are coordinate datasets, + and the values have not already been calculated. + + The example used in this test is for the SMAP SPL3SMP collection, + (SMAP L3 Radiometer Global Daily 36 km EASE-Grid Soil Moisture) + which has a Equal-Area Scalable Earth Grid (EASE-Grid 2.0) CRS for + a projected grid which is lambert_cylindrical_equal_area projection + + """ + smap_varinfo = VarInfoFromDmr( + 'tests/data/SC_SPL3SMP_009.dmr', + 'SPL3SMP', + 'hoss/hoss_config.json', + ) + smap_file_path = 'tests/data/SC_SPL3SMP_009_prefetch.nc4' + expected_index_ranges = {'projected_x': (487, 595), 'projected_y': (9, 38)} + bbox = BBox(2, 54, 42, 72) + + latitude_coordinate = smap_varinfo.get_variable( + '/Soil_Moisture_Retrieval_Data_AM/latitude' + ) + longitude_coordinate = smap_varinfo.get_variable( + '/Soil_Moisture_Retrieval_Data_AM/longitude' + ) + + crs = CRS.from_cf( + { + 'false_easting': 0.0, + 'false_northing': 0.0, + 'longitude_of_central_meridian': 0.0, + 'standard_parallel': 30.0, + 'grid_mapping_name': 'lambert_cylindrical_equal_area', + } + ) + + x_y_extents = { + 'x_min': 192972.56050179302, + 'x_max': 4052423.7705376535, + 'y_min': 5930779.396449475, + 'y_max': 6979878.9118312765, + } + + mock_get_x_y_extents.return_value = x_y_extents + + # When ranges are derived, they are first calculated for x, then y: + mock_get_dimension_index_range.side_effect = [(487, 595), (9, 38)] + + with self.subTest( + 'Projected grid from coordinates gets expected dimension ranges' + ): + with Dataset(smap_file_path, 'r') as smap_prefetch: + self.assertDictEqual( + get_x_y_index_ranges_from_coordinates( + '/Soil_Moisture_Retrieval_Data_AM/surface_flag', + smap_varinfo, + smap_prefetch, + latitude_coordinate, + longitude_coordinate, + {}, + bounding_box=bbox, + shape_file_path=None, + ), + expected_index_ranges, + ) + + # Assertions don't like direct comparisons of numpy arrays, so + # have to extract the call arguments and compare those + mock_get_x_y_extents.assert_called_once_with( + ANY, ANY, crs, shape_file=None, bounding_box=bbox + ) + + actual_x_values = mock_get_x_y_extents.call_args_list[0][0][0] + actual_y_values = mock_get_x_y_extents.call_args_list[0][0][1] + + assert_array_equal(actual_x_values, smap_prefetch['projected_x'][:]) + assert_array_equal(actual_y_values, smap_prefetch['projected_y'][:]) + + self.assertEqual(mock_get_dimension_index_range.call_count, 2) + mock_get_dimension_index_range.assert_has_calls( + [ + call( + ANY, + x_y_extents['x_min'], + x_y_extents['x_max'], + bounds_values=None, + ), + call( + ANY, + x_y_extents['y_min'], + x_y_extents['y_max'], + bounds_values=None, + ), + ] + ) + assert_array_equal( + mock_get_dimension_index_range.call_args_list[0][0][0], + smap_prefetch['projected_x'][:], + ) + assert_array_equal( + mock_get_dimension_index_range.call_args_list[1][0][0], + smap_prefetch['projected_y'][:], + ) + + mock_get_x_y_extents.reset_mock() + mock_get_dimension_index_range.reset_mock() + + with self.subTest( + 'Non projected grid with no coordinates not try to get index ranges' + ): + with Dataset(smap_file_path, 'r') as smap_prefetch: + self.assertDictEqual( + get_x_y_index_ranges_from_coordinates( + 'projected_x', + smap_varinfo, + smap_prefetch, + latitude_coordinate, + longitude_coordinate, + {}, + bounding_box=bbox, + ), + {}, + ) + + mock_get_x_y_extents.assert_not_called() + mock_get_dimension_index_range.assert_not_called() + + with self.subTest('Function does not rederive known index ranges'): + with Dataset(smap_file_path, 'r') as smap_prefetch: + self.assertDictEqual( + get_projected_x_y_index_ranges( + 'Soil_Moisture_Retrieval_Data_AM/surface_flag', + smap_varinfo, + smap_prefetch, + expected_index_ranges, + bounding_box=bbox, + ), + {}, + ) + + mock_get_x_y_extents.assert_not_called() + mock_get_dimension_index_range.assert_not_called() + @patch('hoss.spatial.get_dimension_index_range') @patch('hoss.spatial.get_projected_x_y_extents') def test_get_projected_x_y_index_ranges(