diff --git a/pyproject.toml b/pyproject.toml index ed96ce498..35d63ce3e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,6 +68,7 @@ markers = [ # can be skipped by doing `pytest -m "not slow"` etc. filterwarnings = [ "error:.*removed in a future release of Parcels.*:DeprecationWarning", # Have Parcels DeprecationWarnings fail CI (prevents deprecated items being used in internal code) + "error:::parcels.*", ] [tool.ruff] diff --git a/src/parcels/convert.py b/src/parcels/convert.py index b9eeffa98..c50be82a6 100644 --- a/src/parcels/convert.py +++ b/src/parcels/convert.py @@ -14,7 +14,6 @@ import enum import typing -import warnings from typing import cast import numpy as np @@ -136,7 +135,7 @@ def _maybe_bring_other_depths_to_depth(ds: xr.Dataset): ds[var] = ds[var].rename({old_depth: target}) if "depth" not in ds.dims: - warnings.warn("No depth dimension found in your dataset. Assuming no depth (i.e., surface data).", stacklevel=1) + logger.info("No depth dimension found in your dataset. Assuming no depth (i.e., surface data).", stacklevel=1) ds = ds.expand_dims({"depth": [0]}) ds["depth"] = xr.DataArray([0], dims=["depth"]) return ds @@ -200,11 +199,14 @@ def _set_axis_attrs(ds: xr.Dataset, dim_axis: dict[str, XgcmAxisDirection]): def _ds_rename_using_standard_names(ds: xr.Dataset | ux.UxDataset, name_dict: dict[str, str]) -> xr.Dataset: for standard_name, rename_to in name_dict.items(): - name = ds.cf[standard_name].name - ds = ds.rename({name: rename_to}) - logger.info( - f"cf_xarray found variable {name!r} with CF standard name {standard_name!r} in dataset, renamed it to {rename_to!r} for Parcels simulation." - ) + if standard_name in ds: + ds = ds.rename({standard_name: rename_to}) + else: + name = ds.cf[standard_name].name + ds = ds.rename({name: rename_to}) + logger.info( + f"cf_xarray found variable {name!r} with CF standard name {standard_name!r} in dataset, renamed it to {rename_to!r} for Parcels simulation." + ) return ds @@ -419,7 +421,7 @@ def mitgcm_to_sgrid(*, fields: dict[str, xr.Dataset | xr.DataArray], coords: xr. coords = _pick_expected_coords(coords, _MITGCM_EXPECTED_COORDS) - ds = xr.merge(list(fields.values()) + [coords]) + ds = xr.merge(list(fields.values()) + [coords], compat="override") ds.attrs.clear() # Clear global attributes from the merging ds = _maybe_rename_variables(ds, _MITGCM_VARNAMES_MAPPING) diff --git a/tests/test_advection.py b/tests/test_advection.py index 0a033f213..e98b26014 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -94,7 +94,7 @@ def test_advection_zonal_periodic(): halo = ds.isel(XG=0) halo.lon.values = ds.lon.values[1] + 1 halo.XG.values = ds.XG.values[1] + 2 - ds = xr.concat([ds, halo], dim="XG") + ds = xr.concat([ds, halo], dim="XG", data_vars="all") fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") @@ -286,6 +286,8 @@ def test_moving_eddy(kernel, rtol): if kernel == AdvectionRK45: fieldset.add_context("RK45_tol", rtol) + fieldset.add_context("RK45_min_dt", 1) + fieldset.add_context("RK45_max_dt", 24 * 60 * 60) pset = ParticleSet( fieldset, pclass=DEFAULT_PARTICLES[kernel], x=start_lon, y=start_lat, z=start_z, t=np.timedelta64(0, "s") @@ -325,6 +327,7 @@ def test_decaying_moving_eddy(kernel, rtol): if kernel == AdvectionRK45: fieldset.add_context("RK45_tol", rtol) fieldset.add_context("RK45_min_dt", 10 * 60) + fieldset.add_context("RK45_max_dt", 24 * 60 * 60) pset = ParticleSet(fieldset, pclass=DEFAULT_PARTICLES[kernel], x=start_lon, y=start_lat, t=np.timedelta64(0, "s")) pset.execute(kernel, dt=dt, endtime=endtime) @@ -372,6 +375,8 @@ def test_stommelgyre_fieldset(kernel, rtol, grid_type): if kernel == AdvectionRK45: fieldset.add_context("RK45_tol", rtol) + fieldset.add_context("RK45_min_dt", 1) + fieldset.add_context("RK45_max_dt", 24 * 60 * 60) def UpdateP(particles, fieldset): # pragma: no cover particles.p = fieldset.P[particles.t, particles.z, particles.y, particles.x] @@ -407,6 +412,8 @@ def test_peninsula_fieldset(kernel, rtol, grid_type): if kernel == AdvectionRK45: fieldset.add_context("RK45_tol", rtol) + fieldset.add_context("RK45_min_dt", 1) + fieldset.add_context("RK45_max_dt", 24 * 60 * 60) def UpdateP(particles, fieldset): # pragma: no cover particles.p = fieldset.P[particles.t, particles.z, particles.y, particles.x] diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index ba78da679..00a84c676 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -41,7 +41,7 @@ def field(): temporal_data = np.array([spatial_data, spatial_data + 10, spatial_data + 20]) # each t is +10 from the previous ds = xr.Dataset( - {"U": (["time", "depth", "lat", "lon"], temporal_data)}, + {"P": (["time", "depth", "lat", "lon"], temporal_data)}, coords={ "time": (["time"], [np.timedelta64(t, "s") for t in [0, 2, 4]], {"axis": "T"}), "depth": (["depth"], [0, 1, 2, 3], {"axis": "Z"}), @@ -64,7 +64,7 @@ def field(): vertical_dimensions=(sgrid.FaceNodePadding("ZC", "depth", sgrid.Padding.HIGH),), ), ) - field = FieldSet.from_sgrid_conventions(ds, mesh="flat").U + field = FieldSet.from_sgrid_conventions(ds, mesh="flat").P assert isinstance(field.interp_method, XLinear) return field @@ -192,8 +192,9 @@ def test_interpolation_mesh_type(mesh, npart=10): time = 0.0 u_expected = 1.0 if mesh == "flat" else 1.0 / (1852 * 60 * np.cos(np.radians(lat))) - assert fieldset.U.eval(time, 0, lat, 0) == 1.0 - assert fieldset.V[time, 0, lat, 0] == 0.0 + with pytest.warns(RuntimeWarning, match="Sampling of velocities should normally be done"): + assert fieldset.U.eval(time, 0, lat, 0) == 1.0 + assert fieldset.V[time, 0, lat, 0] == 0.0 u, v = fieldset.UV[time, 0, lat, 0] assert np.isclose(u, u_expected, atol=1e-7) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index d3e8bb811..3cf3e6e80 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -215,13 +215,8 @@ def test_write_timebackward(fieldset, tmp_parquet): df = pd.read_parquet(tmp_parquet) assert df["particle_id"].dtype == "int64" - assert bool( - df.groupby("particle_id") - .apply( - lambda x: (np.diff(x["t"]) < 0).all() # for each particle - set True if it has decreasing time - ) - .all() # ensure for all particles - ) + dt_per_particle = df.groupby("particle_id")["t"].diff().dropna() + assert (dt_per_particle < 0).all() @pytest.mark.xfail @@ -293,7 +288,13 @@ def IncreaseAge(particles, fieldset): # pragma: no cover pset = ParticleSet(fieldset, pclass=AgeParticle, x=npart * [0], y=npart * [0], t=time) ofile = ParticleFile(tmp_parquet, outputdt=outputdt) - pset.execute(IncreaseAge, runtime=np.timedelta64(npart * 2, "s"), dt=np.timedelta64(1, "s"), output_file=ofile) + if outputdt > np.timedelta64(1, "s"): + warning_ctx = pytest.warns(ParticleSetWarning, match="Some of the particles have a start time difference.*") + else: + warning_ctx = does_not_raise() + + with warning_ctx: + pset.execute(IncreaseAge, runtime=np.timedelta64(npart * 2, "s"), dt=np.timedelta64(1, "s"), output_file=ofile) df = parcels.read_particlefile(tmp_parquet) diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index 9b347fe55..244da9643 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -21,7 +21,7 @@ from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels.interpolators import Ux_Velocity, UxConstantFaceConstantZC -from parcels.interpolators._base import ScalarInterpolator +from parcels.interpolators._base import VectorInterpolator from parcels.kernels import AdvectionEE, AdvectionRK2, AdvectionRK4, AdvectionRK4_3D, AdvectionRK45 from tests.common_kernels import DoNothing from tests.utils import DEFAULT_PARTICLES @@ -49,12 +49,14 @@ def zonal_flow_fieldset() -> FieldSet: def test_pset_execute_invalid_arguments(fieldset, fieldset_no_time_interval): - for dt in [np.timedelta64(0, "s"), np.timedelta64(None)]: - with pytest.raises( - ValueError, - match="dt must be a non-zero datetime.timedelta or np.timedelta64 object, got .*", - ): - ParticleSet(fieldset, x=[0.2], y=[5.0], pclass=Particle).execute(AdvectionRK4, dt=dt) + with pytest.raises(RuntimeWarning, match="invalid value encountered in cast.*"): + ParticleSet(fieldset, x=[0.2], y=[5.0], pclass=Particle).execute(AdvectionRK4, dt=np.timedelta64(None)) + + with pytest.raises( + ValueError, + match="dt must be a non-zero datetime.timedelta or np.timedelta64 object, got .*", + ): + ParticleSet(fieldset, x=[0.2], y=[5.0], pclass=Particle).execute(AdvectionRK4, dt=np.timedelta64(0, "s")) with pytest.raises( ValueError, @@ -130,15 +132,28 @@ def test_particleset_endtime_type(fieldset, endtime, expectation): pset.execute(endtime=endtime, dt=np.timedelta64(10, "m"), kernels=DoNothing) +def test_sampleUonly(fieldset): + + def SampleU(particles, fieldset): # pragma: no cover + _ = fieldset.U[particles] + + pset = ParticleSet(fieldset, x=[0.2], y=[5.0]) + with pytest.raises( + RuntimeWarning, + match="Sampling of velocities should normally be done using fieldset.UV or fieldset.UVW object; tread carefully", + ): + pset.execute(SampleU, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) + + def test_particleset_run_to_endtime(fieldset): starttime = fieldset.time_interval.left endtime = fieldset.time_interval.right - def SampleU(particles, fieldset): # pragma: no cover - _ = fieldset.U[particles] + def SampleUV(particles, fieldset): # pragma: no cover + _, _ = fieldset.UV[particles] pset = ParticleSet(fieldset, x=[0.2], y=[5.0], t=[starttime]) - pset.execute(SampleU, endtime=endtime, dt=np.timedelta64(1, "D")) + pset.execute(SampleUV, endtime=endtime, dt=np.timedelta64(1, "D")) assert np.timedelta64(int(pset[0].t), "s") + fieldset.time_interval.left == endtime @@ -149,6 +164,11 @@ def test_particleset_run_RK_to_endtime_fwd_bwd(fieldset, kernel, dt): starttime = fieldset.time_interval.left endtime = fieldset.time_interval.right + if kernel == AdvectionRK45: + fieldset.add_context("RK45_tol", 10) + fieldset.add_context("RK45_min_dt", 1) + fieldset.add_context("RK45_max_dt", 24 * 60 * 60) + # Setting zero velocities to avoid OutofBoundsErrors fieldset.U.data[:] = 0.0 fieldset.V.data[:] = 0.0 @@ -168,11 +188,11 @@ def test_particleset_interpolate_on_domainedge(zonal_flow_fieldset): MyParticle = Particle.add_variable(Variable("var")) - def SampleU(particles, fieldset): # pragma: no cover - particles.var = fieldset.U[particles] + def SampleUV(particles, fieldset): # pragma: no cover + particles.var, _ = fieldset.UV[particles] pset = ParticleSet(fieldset, pclass=MyParticle, x=fieldset.U.grid.lon[-1], y=fieldset.U.grid.lat[-1]) - pset.execute(SampleU, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) + pset.execute(SampleUV, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) np.testing.assert_equal(pset[0].var, 1) @@ -180,7 +200,7 @@ def test_particleset_interpolate_outside_domainedge(zonal_flow_fieldset): fieldset = zonal_flow_fieldset def SampleU(particles, fieldset): # pragma: no cover - particles.dx = fieldset.U[particles] + particles.dx, _ = fieldset.UV[particles] dlat = 1e-3 pset = ParticleSet(fieldset, x=fieldset.U.grid.lon[-1], y=fieldset.U.grid.lat[-1] + dlat) @@ -301,7 +321,7 @@ def test_some_particles_throw_outofbounds(zonal_flow_fieldset): def test_delete_on_all_errors(fieldset): def MoveRight(particles, fieldset): # pragma: no cover particles.dx += 1 - fieldset.U[particles.t, particles.z, particles.y, particles.x, particles] + fieldset.UV[particles.t, particles.z, particles.y, particles.x, particles] def DeleteAllErrorParticles(particles, fieldset): # pragma: no cover particles[particles.state > 20].state = StatusCode.Delete @@ -316,7 +336,7 @@ def test_some_particles_throw_outoftime(fieldset): pset = ParticleSet(fieldset, x=np.zeros_like(time), y=np.zeros_like(time), t=time) def FieldAccessOutsideTime(particles, fieldset): # pragma: no cover - fieldset.U[particles.t + 400 * 86400, particles.z, particles.y, particles.x, particles] + fieldset.UV[particles.t + 400 * 86400, particles.z, particles.y, particles.x, particles] with pytest.raises(OutsideTimeInterval): pset.execute(FieldAccessOutsideTime, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(10, "D")) @@ -329,17 +349,18 @@ def test_raise_general_error(): ... def test_errorinterpolation(fieldset): - class NaNInterpolator(ScalarInterpolator): # pragma: no cover + class NaNInterpolator(VectorInterpolator): # pragma: no cover def interp(self, particle_positions, grid_positions, field): - return np.nan * np.zeros_like(particle_positions["x"]) + nanvals = np.nan * np.zeros_like(particle_positions["x"]) + return nanvals, nanvals, nanvals - def SampleU(particles, fieldset): # pragma: no cover - fieldset.U[particles.t, particles.z, particles.y, particles.x, particles] + def SampleUV(particles, fieldset): # pragma: no cover + fieldset.UV[particles.t, particles.z, particles.y, particles.x, particles] - fieldset.U.interp_method = NaNInterpolator() + fieldset.UV.interp_method = NaNInterpolator() pset = ParticleSet(fieldset, x=[0, 2], y=[0, 0]) with pytest.raises(FieldInterpolationError): - pset.execute(SampleU, runtime=np.timedelta64(2, "s"), dt=np.timedelta64(1, "s")) + pset.execute(SampleUV, runtime=np.timedelta64(2, "s"), dt=np.timedelta64(1, "s")) def test_execution_check_stopallexecution(fieldset): @@ -357,7 +378,7 @@ def test_execution_recover_out_of_bounds(fieldset): npart = 2 def MoveRight(particles, fieldset): # pragma: no cover - fieldset.U[particles.t, particles.z, particles.y, particles.x + 0.1, particles] + fieldset.UV[particles.t, particles.z, particles.y, particles.x + 0.1, particles] particles.dx += 0.1 def MoveLeft(particles, fieldset): # pragma: no cover