diff --git a/test/grid/geometry/test_geometry.py b/test/grid/geometry/test_geometry.py index c50a765f2..e602e1e0d 100644 --- a/test/grid/geometry/test_geometry.py +++ b/test/grid/geometry/test_geometry.py @@ -1,11 +1,13 @@ import numpy as np -import numpy.testing as nt import pytest +import xarray as xr import uxarray as ux from uxarray.constants import ERROR_TOLERANCE, INT_FILL_VALUE -from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz, _xyz_to_lonlat_rad -from uxarray.grid.geometry import _pole_point_inside_polygon_cartesian +from uxarray.grid.coordinates import _lonlat_rad_to_xyz, _normalize_xyz +from uxarray.grid.geometry import haversine_distance, _pole_point_inside_polygon_cartesian +from uxarray.grid.point_in_face import _face_contains_point +from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array def test_antimeridian_crossing(): @@ -197,13 +199,6 @@ def test_haversine_distance_creation(): np.testing.assert_allclose(distance, expected_distance, atol=ERROR_TOLERANCE) -from uxarray.constants import ERROR_TOLERANCE -from uxarray.grid.coordinates import _lonlat_rad_to_xyz -from uxarray.grid.utils import _get_cartesian_face_edge_nodes_array -from uxarray.grid.point_in_face import _face_contains_point -from uxarray.grid.geometry import haversine_distance - - def test_engine_geodataframe(gridpath): uxgrid = ux.open_grid(gridpath("ugrid", "geoflow-small", "grid.nc")) @@ -215,3 +210,56 @@ def test_periodic_elements_geodataframe(gridpath): uxgrid = ux.open_grid(gridpath("ugrid", "geoflow-small", "grid.nc")) for periodic_elements in ["ignore", "exclude", "split"]: gdf = uxgrid.to_geodataframe(periodic_elements=periodic_elements) + + +def test_geodataframe_crashes_when_nnode_equals_nface(): + """Ensure UxDataArray.to_geodataframe() crashes for non-face_centered data, even if n_node==n_face. + regression test for issue #1616. + """ + # Below is a visualization of the example here, with ni node i, fj face j: + # (n1)----------------(n2) + # | %% f1 =%/| + # | %% =% / | + # | %% =% / | + # | %% =% f2 / | + # | (n3)----(n4) f3| + # | f0 %% =% f4 \ | + # | %% =% \ | + # | %% f5 =% \| + # (n0)----------------(n5) + XA, XB, XC, XF = 0, 10, 20, 30 + YA, YB, YF = 0, 10, 25 + node_lonlat = [ + [XA, YA], # n0 + [XA, YF], # n1 + [XF, YF], # n2 + [XB, YB], # n3 + [XC, YB], # n4 + [XF, YA], # n5 + ] + node_lon = [n[0] for n in node_lonlat] + node_lat = [n[1] for n in node_lonlat] + face_node_connectivity = [ + [0, 1, 3], # f0 + [1, 3, 2], # f1 + [3, 4, 2], # f2 + [4, 5, 2], # f3 + [3, 5, 4], # f4 + [0, 5, 3], # f5 + ] + # convert to xarray + ds_vars = { + 'node_lon': xr.DataArray(node_lon, dims=['n_node']), + 'node_lat': xr.DataArray(node_lat, dims=['n_node']), + 'face_node_connectivity': xr.DataArray(face_node_connectivity, dims=['n_face', 'n_max_face_nodes']), + } + ds = xr.Dataset(ds_vars) + # convert to uxarray + uxgrid = ux.Grid(ds, 'UGRID') # putting 'UGRID' avoids raising warning but does not affect results. + assert uxgrid.n_node == uxgrid.n_face + # make array of values, convert to uxarray + vals = xr.DataArray([100,200,300,400,500,600], dims=['n_node']) + uxarr = ux.UxDataArray(vals, uxgrid=uxgrid) + # ensure to_geodataframe crashes for non-face_centered data, even if n_node==n_face + with pytest.raises(ValueError): + uxarr.to_geodataframe() diff --git a/test/grid/integrate/test_basic.py b/test/grid/integrate/test_basic.py index e45df9dfd..c90db58e7 100644 --- a/test/grid/integrate/test_basic.py +++ b/test/grid/integrate/test_basic.py @@ -1,6 +1,7 @@ import numpy as np import numpy.testing as nt import pytest +import xarray as xr import uxarray as ux from uxarray.constants import ERROR_TOLERANCE @@ -26,3 +27,55 @@ def test_multi_dim(gridpath): integral = uxda.integrate() assert integral.ndim == len(dims) - 1 nt.assert_almost_equal(integral, np.ones((5, 5)) * 4 * np.pi) + +def test_integrate_crashes_when_nnode_equals_nface(): + """Ensure UxDataArray.integrate() crashes for non-face_centered data, even if n_node==n_face. + regression test for issue #1616. + """ + # Below is a visualization of the example here, with (ni)=node i; fj=face j: + # (n1)----------------(n2) + # | %% f1 =%/| + # | %% =% / | + # | %% =% / | + # | %% =% f2 / | + # | (n3)----(n4) f3| + # | f0 %% =% f4 \ | + # | %% =% \ | + # | %% f5 =% \| + # (n0)----------------(n5) + XA, XB, XC, XF = 0, 10, 20, 30 + YA, YB, YF = 0, 10, 25 + node_lonlat = [ + [XA, YA], # n0 + [XA, YF], # n1 + [XF, YF], # n2 + [XB, YB], # n3 + [XC, YB], # n4 + [XF, YA], # n5 + ] + node_lon = [n[0] for n in node_lonlat] + node_lat = [n[1] for n in node_lonlat] + face_node_connectivity = [ + [0, 1, 3], # f0 + [1, 3, 2], # f1 + [3, 4, 2], # f2 + [4, 5, 2], # f3 + [3, 5, 4], # f4 + [0, 5, 3], # f5 + ] + # convert to xarray + ds_vars = { + 'node_lon': xr.DataArray(node_lon, dims=['n_node']), + 'node_lat': xr.DataArray(node_lat, dims=['n_node']), + 'face_node_connectivity': xr.DataArray(face_node_connectivity, dims=['n_face', 'n_max_face_nodes']), + } + ds = xr.Dataset(ds_vars) + # convert to uxarray + uxgrid = ux.Grid(ds, 'UGRID') # putting 'UGRID' avoids raising warning but does not affect results. + assert uxgrid.n_node == uxgrid.n_face + # make array of values, convert to uxarray + vals = xr.DataArray([100,200,300,400,500,600], dims=['n_node']) + uxarr = ux.UxDataArray(vals, uxgrid=uxgrid) + # ensure integrate() crashes for non-face_centered data, even if n_node==n_face + with pytest.raises(ValueError): + uxarr.integrate() diff --git a/uxarray/core/dataarray.py b/uxarray/core/dataarray.py index 8f41a1627..69b12c0de 100644 --- a/uxarray/core/dataarray.py +++ b/uxarray/core/dataarray.py @@ -254,7 +254,7 @@ def to_geodataframe( f"for face-centered data." ) - if self.values.size == self.uxgrid.n_face: + if self._face_centered(): gdf, non_nan_polygon_indices = self.uxgrid.to_geodataframe( periodic_elements=periodic_elements, projection=projection, @@ -291,22 +291,11 @@ def to_geodataframe( gdf[var_name] = _data - elif self.values.size == self.uxgrid.n_node: - raise DataCenteringError( - f"Data Variable with size {self.values.size} does not match the number of faces " - f"({self.uxgrid.n_face}. Current size matches the number of nodes. Consider running " - f"``UxDataArray.topological_mean(destination='face') to aggregate the data onto the faces." - ) - elif self.values.size == self.uxgrid.n_edge: - raise DataCenteringError( - f"Data Variable with size {self.values.size} does not match the number of faces " - f"({self.uxgrid.n_face}. Current size matches the number of edges." - ) else: - # data is not mapped to - raise DimensionError( - f"Data Variable with size {self.values.size} does not match the number of faces " - f"({self.uxgrid.n_face}." + raise DataCenteringError( + f"to_geodataframe() expects face_centered data; got {self.data_location} data " + f"(with sizes={dict(**self.sizes)}). Consider running " + "``UxDataArray.topological_mean(destination='face')`` to aggregate the data onto faces." ) return gdf @@ -613,28 +602,25 @@ def integrate( >>> uxds = ux.open_dataset("grid.ug", "centroid_pressure_data_ug") >>> integral = uxds["psi"].integrate() """ - if self.values.shape[-1] == self.uxgrid.n_face: + # TODO: support integration regardless of n_face dimension position, + # and remove the self.dims[-1] == "n_face" check. + # (uxarray/xarray features should be agnostic to dimension positions.) + if self._face_centered() and self.dims[-1] == "n_face": face_areas = self.uxgrid.face_areas.values # perform dot product between face areas and last dimension of data integral = np.einsum("i,...i", face_areas, self.values) - elif self.values.shape[-1] == self.uxgrid.n_node: - raise DataCenteringError( - "Integrating data mapped to each node not yet supported." - ) - - elif self.values.shape[-1] == self.uxgrid.n_edge: + elif not self._face_centered(): raise DataCenteringError( - "Integrating data mapped to each edge not yet supported." + "Integration of non-face_centered data is not yet supported. " + f"(Got {self.data_location} data with sizes={dict(**self.sizes)})" ) else: raise DimensionError( - f"The final dimension of the data variable does not match the number of nodes, edges, " - f"or faces. Expected one of " - f"{self.uxgrid.n_node}, {self.uxgrid.n_edge}, or {self.uxgrid.n_face}, " - f"but received {self.values.shape[-1]}" + "Integration of data with n_face not as the final dimension is not yet supported." + f"Got face_centered data, but the final dimension was {self.dims[-1]}, not 'n_face'." ) # construct a uxda with integrated quantity