Skip to content

Commit

Permalink
DAS-2232 - moved new methods to a separate file
Browse files Browse the repository at this point in the history
  • Loading branch information
@sudha-murthy
sudha-murthy committed Oct 4, 2024
1 parent 1e7bc35 commit 36e15c7
Show file tree
Hide file tree
Showing 3 changed files with 188 additions and 48 deletions.
5 changes: 2 additions & 3 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down
82 changes: 37 additions & 45 deletions hoss/coordinate_utilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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')
Expand Down
149 changes: 149 additions & 0 deletions tests/unit/test_spatial.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
)


Expand Down Expand Up @@ -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(
Expand Down

0 comments on commit 36e15c7

Please sign in to comment.