diff --git a/docs/user_guide/examples/tutorial_schism.ipynb b/docs/user_guide/examples/tutorial_schism.ipynb index 406938c96..a81b4cb97 100644 --- a/docs/user_guide/examples/tutorial_schism.ipynb +++ b/docs/user_guide/examples/tutorial_schism.ipynb @@ -26,7 +26,7 @@ "1. The mesh coordinates are in a **projected coordinate system (meters)**, not longitude/latitude, so\n", " we build the grid as a *flat* (Cartesian) mesh.\n", "2. Velocities are stored at **mesh nodes** over a vertical column of layers, ordered **bottom → surface**.\n", - "3. SCHISM uses a **localized vertical grid (LSC2)**: the number of valid vertical levels varies from\n", + "3. SCHISM uses a **localised vertical grid (LSC2)**: the number of valid vertical levels varies from\n", " node to node, and levels below the seabed are stored as `NaN`. We use this to flag particles that\n", " end up below the local bathymetry and stop advecting them.\n", "\n", @@ -147,7 +147,7 @@ "\n", "* **Layer ordering.** SCHISM stores levels **bottom → surface**, while Parcels expects depth increasing\n", " downward from the surface, so we reverse the layer axis.\n", - "* **Vertical coordinate.** SCHISM's LSC2 vertical grid (*Localized Sigma Coordinates with Shaved\n", + "* **Vertical coordinate.** SCHISM's LSC2 vertical grid (*Localised Sigma Coordinates with Shaved\n", " cells*) varies with horizontal position: each node has its own layer depths, and even its own\n", " *number* of levels (the true depths are written to a separate `zCoordinates` output, not used here).\n", " **Parcels does not currently support a vertical grid that varies with lateral position.** `UxGrid`\n", diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 9355c145e..ca6f3aea7 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -70,12 +70,12 @@ def __init__(self, models: list[ModelData]): for model in models: if not isinstance(model, ModelData): raise ValueError(f"Expected `model` to be a ModelData object. Got {model}") - # assert_compatible_calendars(fields) self.models = models self.constant_model: StructuredModelData | None = None self._fields: dict[str, Field | VectorField] | None = None self.reconstruct_fields() + assert_compatible_calendars(self.fields.values()) self.context: dict[str, float] = {} assert_models_have_same_mesh(self.models) diff --git a/src/parcels/_core/particlefile.py b/src/parcels/_core/particlefile.py index fb02885c6..ff397f17f 100644 --- a/src/parcels/_core/particlefile.py +++ b/src/parcels/_core/particlefile.py @@ -166,7 +166,6 @@ def write(self, pset: ParticleSet | ParticleSetView, t, indices=None): ) if isinstance(t, (np.timedelta64, np.datetime64)): - print(t) if self._time_interval is None: raise ValueError( "The time interval is unknown, so datetime/timedelta times cannot be converted. " @@ -174,10 +173,13 @@ def write(self, pset: ParticleSet | ParticleSetView, t, indices=None): ) t = timedelta_to_float(t - self._time_interval.left) vars_to_write = _get_vars_to_write(pclass) - if indices is None: - indices_to_write = _to_write_particles(particle_data, t) - else: + if indices is not None: indices_to_write = indices + elif isinstance(pset, ParticleSetView): + # A view already selects the particles to write, so `t` is not used to filter + indices_to_write = np.flatnonzero(pset._index) + else: + indices_to_write = _to_write_particles(particle_data, t) self._writer.write_table( pa.table({v.name: pa.array(particle_data[v.name][indices_to_write]) for v in vars_to_write}), diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 80bada25d..b58783086 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -256,9 +256,9 @@ def get_axis_dim(self, axis: ptyping.XgridAxis) -> int: assert self.sgrid_metadata.vertical_dimensions is not None return get_cell_count_along_dim(self._ds, self.sgrid_metadata.vertical_dimensions[0]) - def localize( - self, position: dict[ptyping.XgridAxis, tuple[int, float]], dims: list[str] - ) -> dict[str, tuple[int, float]]: + def localise( + self, position: dict[ptyping.XgridAxis, dict[str, np.ndarray]], dims: list[str] + ) -> dict[str, dict[str, np.ndarray]]: """ Uses the grid context (i.e., the staggering of the grid) to convert a position relative to the f-points in the grid to a position relative to the staggered grid the array @@ -271,24 +271,28 @@ def localize( Parameters ---------- position : dict - A mapping of the axis to a tuple of (index, barycentric coordinate) for the - f-points in the grid. + A mapping of the axis to a dict with the arrays ``"index"`` and ``"bcoord"`` (barycentric + coordinate) for the f-points in the grid, as returned by ``search``. dims : list[str] A list of dimension names that the DataArray is defined on. This is used to determine the staggering of the grid and which axis each dimension corresponds to. Returns ------- - dict[str, tuple[int, float]] - A mapping of the dimension names to a tuple of (index, barycentric coordinate) for + dict[str, dict[str, np.ndarray]] + A mapping of the dimension names to a dict with the arrays ``"index"`` and ``"bcoord"`` for the staggered grid the DataArray is defined on. Example ------- - >>> position = {"X": (5, 0.51), "Y": (10, 0.25), "Z": (3, 0.75)} + >>> position = { + ... "X": {"index": np.array([5]), "bcoord": np.array([0.51])}, + ... "Y": {"index": np.array([10]), "bcoord": np.array([0.25])}, + ... "Z": {"index": np.array([3]), "bcoord": np.array([0.75])}, + ... } >>> dims = ["time", "depth", "YC", "XC"] - >>> grid.localize(position, dims) - {'depth': (3, 0.75), 'YC': (9, 0.75), 'XC': (5, 0.01)} + >>> grid.localise(position, dims) + {'depth': {'index': array([3]), 'bcoord': array([0.75])}, 'YC': {'index': array([9]), 'bcoord': array([0.75])}, 'XC': {'index': array([5]), 'bcoord': array([0.01])}} """ dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata) axis_to_var = {dim_to_axis[dim]: dim for dim in dims if dim in dim_to_axis} @@ -297,12 +301,13 @@ def localize( } return { axis_to_var[axis]: _convert_center_pos_to_fpoint( - index=index, - bcoord=bcoord, + index=pos["index"], + bcoord=pos["bcoord"], position=var_positions[axis], f_point_position=self._fpoint_info[axis], ) - for axis, (index, bcoord) in position.items() + for axis, pos in position.items() + if axis in axis_to_var } @property @@ -591,29 +596,29 @@ def assert_valid_depth(da_depth): def _convert_center_pos_to_fpoint( *, - index: int, - bcoord: float, + index: np.ndarray, + bcoord: np.ndarray, position: ptyping.GridPosition, f_point_position: sgrid.Padding, -) -> tuple[int, float]: +) -> dict[str, np.ndarray]: """Converts a physical position relative to the cell edges defined in the grid to be relative to the face center. - This is used to "localize" a position to be relative to the staggered grid at which the field is defined, so that + This is used to "localise" a position to be relative to the staggered grid at which the field is defined, so that it can be easily interpolated. This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm). """ if position != "face": # Data is already defined on the F nodes - return index, bcoord + return {"index": index, "bcoord": bcoord} bcoord = bcoord - 0.5 - if bcoord < 0: - bcoord += 1.0 - index -= 1 + wrapped = bcoord < 0 + bcoord = np.where(wrapped, bcoord + 1.0, bcoord) + index = np.where(wrapped, index - 1, index) # Correct relative to the f-point position # Padding.BOTH was "inner", Padding.LOW was "right" in xgcm vocabulary if f_point_position in (sgrid.Padding.BOTH, sgrid.Padding.LOW): - index += 1 + index = index + 1 - return index, bcoord + return {"index": index, "bcoord": bcoord} diff --git a/src/parcels/_repr_utils.py b/src/parcels/_repr_utils.py index ab5918598..929a21a1d 100644 --- a/src/parcels/_repr_utils.py +++ b/src/parcels/_repr_utils.py @@ -104,7 +104,7 @@ def particleset_repr(pset: ParticleSet) -> str: def particlesetview_repr(pview: Any) -> str: """Return a pretty repr for ParticleSetView""" - time_string = "not_yet_set" if pview.time is None or np.isnan(pview.time) else f"{pview.time:f}" + time_string = "not_yet_set" if pview.t is None or np.isnan(pview.t) else f"{pview.t:f}" out = f"P[{pview.particle_id}]: time={time_string}, z={pview.z:f}, y={pview.y:f}, x={pview.x:f}" vars = [v.name for v in pview._pclass.variables if v.to_write is True and v.name not in ["z", "y", "x", "time"]] for var in vars: diff --git a/tests/test_field.py b/tests/test_field.py index 2272a5cb3..7def3ca8b 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -2,6 +2,7 @@ import numpy as np import pytest +import xarray as xr from re_assert import Matches from parcels import Field, VectorField @@ -82,7 +83,16 @@ def test_field_time_interval(): def test_vectorfield_init_different_time_intervals(): # Tests that a VectorField raises a ValueError if the component fields have different time domains. - ... + dsU = simple_UV_dataset(mesh_str="flat")[["U", "grid"]] + dsV = simple_UV_dataset(mesh_str="flat")[["V"]] + + dsV = dsV.rename({"time": "timeV"}) + dsV["timeV"] = dsV["timeV"] + np.timedelta64(1, "D") + dsV["timeV"].attrs["axis"] = "T" + + ds = xr.merge([dsU, dsV]) + with pytest.raises(ValueError, match="Multiple time coordinates"): + FieldSet.from_sgrid_conventions(ds) def test_field_invalid_interpolator(): @@ -242,15 +252,3 @@ def test_field_eval_out_of_bounds_unstructured(field_name, location, expected): match="Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.", ): np.testing.assert_allclose(field.eval(*location), expected) - - -def test_field_unstructured_grid_creation(): ... - - -def test_field_interpolation(): ... - - -def test_field_interpolation_out_of_spatial_bounds(): ... - - -def test_field_interpolation_out_of_time_bounds(): ... diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 10f7c7045..ea126e664 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -13,11 +13,11 @@ import parcels.tutorial import tests from parcels import Particle, ParticleFile, ParticleSet, Variable, convert -from parcels._core.fieldset import FieldSet, IncompatibleMeshesException, _datetime_to_msg +from parcels._core.fieldset import CalendarError, FieldSet, IncompatibleMeshesException, _datetime_to_msg from parcels._core.mesh import SphericalMesh from parcels._core.model import _default_vector_field_components, _infer_mesh_from_sgrid_data from parcels._core.utils.time import TimeInterval -from parcels._datasets.structured.generated import simple_UV_sigma_dataset +from parcels._datasets.structured.generated import simple_UV_dataset, simple_UV_sigma_dataset from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.structured.generic import datasets_sgrid from parcels._datasets.structured.strategies import reduced_dims_dataset, sgrid_datasets_varying_dims @@ -272,7 +272,16 @@ def test_multi_model_nonoverlapping_time_interval(): def test_fieldset_add_incompatible_calendars(): # tests the adding of fieldsets that have incompatible calendars - ... + ds1 = datasets_structured["ds_2d_left"][["U_A_grid", "V_A_grid", "grid"]] + fieldset = FieldSet.from_sgrid_conventions(ds1) + + ds2 = ds1.copy().rename({"U_A_grid": "U2", "V_A_grid": "V2"}) + # Replace the numpy datetime64 time axis by an array of cftime objects with a different calendar + times = [cftime.DatetimeNoLeap(2000, 1, 1 + i) for i in range(ds2["time"].size)] + ds2["time"] = (ds2["time"].dims, np.array(times, dtype=object), ds2["time"].attrs) + + with pytest.raises(CalendarError): + fieldset + FieldSet.from_sgrid_conventions(ds2) @pytest.mark.parametrize( @@ -295,7 +304,11 @@ def test_datetime_to_msg(input_, expected): def test_fieldset_samegrids_UV(): """Test that if a simple fieldset with U and V is created, that only one grid object is defined.""" - ... + ds = simple_UV_dataset(mesh_str="flat") + fieldset = FieldSet.from_sgrid_conventions(ds) + + assert fieldset.U.grid is fieldset.V.grid + assert len(fieldset.gridset) == 1 def test_fieldset_grid_deduplication(): @@ -305,7 +318,15 @@ def test_fieldset_grid_deduplication(): When grid deduplication is actually implemented, this might need to be refactored into multiple tests (/more might be needed). """ - ... + ds = simple_UV_dataset(mesh_str="flat") + fieldset = FieldSet.from_sgrid_conventions(ds) + + ds2 = simple_UV_dataset(mesh_str="flat").rename({"U": "U2", "V": "V2"}) + fieldset += FieldSet.from_sgrid_conventions(ds2) + assert fieldset.U.grid is fieldset.V.grid + assert fieldset.U2.grid is fieldset.V2.grid + assert fieldset.U.grid is not fieldset.U2.grid + assert fieldset.U.grid is not fieldset.V2.grid def test_fieldset_from_icon(): diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index d1c03c96a..f5b13a865 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -8,7 +8,6 @@ import pyarrow as pa import pyarrow.parquet as pq import pytest -import xarray as xr from re_assert import Matches import parcels.tutorial @@ -156,46 +155,6 @@ def test_write_dtypes_pfile(fieldset, tmp_parquet): assert tab[f"v_{d.__name__}"].type == pa.from_numpy_dtype(d) -@pytest.mark.skip(reason="Pending ParticleFile refactor; see issue #2386") -@pytest.mark.parametrize("dt", [-np.timedelta64(1, "s"), np.timedelta64(1, "s")]) -@pytest.mark.parametrize("maxvar", [2, 4, 10]) -def test_pset_repeated_release_delayed_adding_deleting(fieldset, tmp_parquet, dt, maxvar): - """Tests that if particles are released and deleted based on age that resulting output file is correct.""" - npart = 10 - fieldset.add_context("maxvar", maxvar) - - MyParticle = Particle.add_variable( - [Variable("sample_var", initial=0.0), Variable("v_once", dtype=np.float64, initial=0.0)] - ) - - pset = ParticleSet( - fieldset, - x=np.zeros(npart), - y=np.zeros(npart), - pclass=MyParticle, - t=fieldset.time_interval.left + [np.timedelta64(i + 1, "s") for i in range(npart)], - ) - pfile = ParticleFile(tmp_parquet, outputdt=abs(dt)) - - def IncrLon(particles, fieldset): # pragma: no cover - particles.sample_var += 1.0 - particles.state = np.where( - particles.sample_var > fieldset.maxvar, - StatusCode.Delete, - particles.state, - ) - - for _ in range(npart): - pset.execute(IncrLon, dt=dt, runtime=np.timedelta64(1, "s"), output_file=pfile) - - ds = xr.open_zarr(tmp_parquet) - samplevar = ds["sample_var"][:] - assert samplevar.shape == (npart, min(maxvar, npart + 1)) - # test whether samplevar[:, k] = k - for k in range(samplevar.shape[1]): - assert np.allclose([p for p in samplevar[:, k] if np.isfinite(p)], k + 1) - - def test_file_warnings(fieldset, tmp_parquet): pset = ParticleSet(fieldset, x=[0, 0], y=[0, 0], t=[np.timedelta64(0, "s"), np.timedelta64(1, "s")]) pfile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(2, "s")) diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index 9343307a8..2b37e325a 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -390,12 +390,6 @@ def FieldAccessOutsideTime(particles, fieldset): # pragma: no cover pset.execute(FieldAccessOutsideTime, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(10, "D")) -def test_raise_grid_searching_error(): ... - - -def test_raise_general_error(): ... - - def test_errorinterpolation(fieldset): class NaNInterpolator(VectorInterpolator): # pragma: no cover def interp(self, particle_positions, grid_positions, field): @@ -571,24 +565,23 @@ def test_uxstommelgyre_multiparticle_pset_execute(): ) -@pytest.mark.xfail(reason="Output file not implemented yet") -def test_uxstommelgyre_pset_execute_output(): +def test_uxstommelgyre_pset_execute_output(tmp_parquet): ds = datasets_unstructured["stommel_gyre_delaunay"] fieldset = FieldSet.from_ugrid_conventions(ds) pset = ParticleSet( fieldset, - x=[30.0], - y=[5.0], - z=[50.0], - t=[0.0], + x=30.0, + y=5.0, + z=50.0, + t=np.timedelta64(0, "s"), pclass=Particle, ) output_file = ParticleFile( - name="stommel_uxarray_particles.zarr", # the file name + path=tmp_parquet, # the file name outputdt=np.timedelta64(5, "m"), # the time step of the outputs ) pset.execute( - runtime=np.timedelta64(10, "m"), dt=np.timedelta64(60, "s"), kernels=AdvectionEE, output_file=output_file + runtime=np.timedelta64(10, "m"), dt=np.timedelta64(150, "s"), kernels=AdvectionEE, output_file=output_file ) diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index 465284215..17b2efccd 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -20,7 +20,7 @@ XGrid, _transpose_xfield_data_to_tzyx, ) -from parcels._datasets.structured.generic import X, Y, Z, datasets, datasets_sgrid +from parcels._datasets.structured.generic import X, Y, Z, datasets from tests import utils GridTestCase = namedtuple("GridTestCase", ["ds", "attr", "expected"]) @@ -50,13 +50,10 @@ def assert_equal(actual, expected): assert_allclose(actual, expected) -@pytest.mark.skip( - "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" -) # TODO: Remove or replace @pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) def test_grid_init_param_types(ds): - with pytest.raises(ValueError, match="Invalid value 'invalid'. Valid options are.*"): - XGrid.from_dataset(ds, mesh="invalid") + with pytest.raises(TypeError, match="mesh must be a FlatMesh or SphericalMesh object. Got mesh.*"): + FieldSet.from_sgrid_conventions(ds, mesh="invalid") def test_xgrid_repr(fieldset): @@ -85,20 +82,17 @@ def test_xgrid_mesh_inferred_from_units(ds, mesh): assert isinstance(grid.mesh, mesh) -@pytest.mark.skip( - "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" -) # TODO: Remove or replace @pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) def test_transpose_xfield_data_to_tzyx(ds): da = ds["data_g"] - grid = XGrid.from_dataset(ds) + grid = FieldSet.from_sgrid_conventions(ds).data_g.grid all_combinations = (itertools.combinations(da.dims, n) for n in range(len(da.dims))) all_combinations = itertools.chain(*all_combinations) for subset_dims in all_combinations: isel = {dim: 0 for dim in subset_dims} da_subset = da.isel(isel, drop=True) - da_test = _transpose_xfield_data_to_tzyx(da_subset, grid.xgcm_grid) + da_test = _transpose_xfield_data_to_tzyx(da_subset, grid.sgrid_metadata) utils.assert_valid_field_data(da_test, grid) @@ -109,84 +103,49 @@ def test_xgrid_get_axis_dim(fieldset): assert grid.get_axis_dim("X") == X - 1 -def test_invalid_xgrid_field_array(): - """Stress test initialiser by creating incompatible datasets that test the edge cases""" - ... - - -@pytest.mark.skip( - "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" -) # TODO: Remove or replace def test_invalid_lon_lat(): """Stress test the grid initialiser by creating incompatible datasets that test the edge cases""" ds = datasets["ds_2d_left"].copy() ds["lon"], ds["lat"] = xr.broadcast(ds["YC"], ds["XC"]) + ds["lon"].attrs["units"] = "meters" + ds["lat"].attrs["units"] = "meters" with pytest.raises( ValueError, - match=r".*is defined on the center of the grid, but must be defined on the F points\.", + match=r".*is defined on the faces of the grid, but must be defined on the F nodes\.", ): - XGrid.from_dataset(ds) + FieldSet.from_sgrid_conventions(ds) ds = datasets["ds_2d_left"].copy() ds["lon"], _ = xr.broadcast(ds["YG"], ds["XG"]) + ds["lon"].attrs["units"] = "meters" with pytest.raises( ValueError, match=r".*have different dimensionalities\.", ): - XGrid.from_dataset(ds) + FieldSet.from_sgrid_conventions(ds) ds = datasets["ds_2d_left"].copy() ds["lon"], ds["lat"] = xr.broadcast(ds["YG"], ds["XG"]) ds["lon"], ds["lat"] = ds["lon"].transpose(), ds["lat"].transpose() + ds["lon"].attrs["units"] = "meters" + ds["lat"].attrs["units"] = "meters" with pytest.raises( ValueError, match=r".*must be defined on the X and Y axes and transposed to have dimensions in order of Y, X\.", ): - XGrid.from_dataset(ds) + FieldSet.from_sgrid_conventions(ds) -@pytest.mark.skip( - "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" -) # TODO: Remove or replace def test_invalid_depth(): ds = datasets["ds_2d_left"].copy() ds = ds.reindex({"ZG": ds.ZG[::-1]}) with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"): - XGrid.from_dataset(ds) - - -@pytest.mark.skip( - "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646" -) # TODO: I think we can just rely on the SGRID metadata for this (which already has robust error messaging). How should discrepencies between SGRID and axis attr be handled? -def test_dim_with_duplicate_axis(): - ds = datasets_sgrid["ds_2d_padded_low"].copy() - - # Add an extra Z axis - ds = ds[["data_g", "grid"]] - ds["data_g"] = ds["data_g"].expand_dims("vertical_dimensions_dim2", 1) - - z = ds["vertical_dimensions_dim2"] - z.attrs.update({"axis": "Z"}) - ds["vertical_dimensions_dim2"] = z - - # TODO: Clean up this attribute setting (really, this should be on the source datasets) - lon = ds["lon"] - lat = ds["lat"] - lon.attrs.update({"units": "metres"}) - lat.attrs.update({"units": "metres"}) - ds["lon"] = lon - ds["lat"] = lat - - with pytest.raises(ValueError, match="Two dimensions .*provide values in the axis direction 'Z'."): FieldSet.from_sgrid_conventions(ds) -# TODO restructure: Look into the test below -# remove: eval with timedelta64 time fails; TimeInterval.is_all_time_in_interval expects float (seconds) not timedelta64; time arg type requirement changed -@pytest.mark.skip("remove: see comment above") @pytest.mark.parametrize("ds", [datasets["ds_2d_left"]]) def test_vertical1D_field(ds): ds = ds.drop_vars(set(ds.data_vars) - {"grid"}) @@ -304,7 +263,6 @@ def test_search_1d_array_some_out_of_bounds(search_1d, array, x, expected_xi): np.testing.assert_array_equal(xi, expected_xi) -@pytest.mark.xfail(reason="grid.localize not yet adapted to xi, xsi, yi, eta, zi, zeta position keys") @pytest.mark.parametrize( "ds, da_name, expected", [ @@ -350,11 +308,16 @@ def test_search_1d_array_some_out_of_bounds(search_1d, array, x, expected_xi): ), ], ) -def test_xgrid_localize_zero_position(ds, da_name, expected): - """Test localize function using left and right datasets.""" - grid = XGrid.from_dataset(ds) +def test_xgrid_localise_zero_position(ds, da_name, expected): + """Test localise function using left and right datasets.""" + fieldset = FieldSet.from_sgrid_conventions(ds) + grid = fieldset.models[0].grid + da = ds[da_name] position = grid.search(0, 0, 0) - local_position = grid.localize(position, da.dims) - assert local_position == expected, f"Expected {expected}, got {local_position}" + local_position = grid.localise(position, da.dims) + assert local_position.keys() == expected.keys() + for dim, (index, bcoord) in expected.items(): + np.testing.assert_array_equal(local_position[dim]["index"], [index]) + np.testing.assert_allclose(local_position[dim]["bcoord"], [bcoord])