Merge pull request #2341 from reint-fischer/decode_times_issue_v4

Fix decode_times issue by using signed integers for time writing

Author: reint-fischer

docs/examples/tutorial_Argofloats.ipynb

@@ -41,11 +41,9 @@
     "\n",
     "    def SinkingPhase(p):\n",
     "        \"\"\"Phase 0: Sinking with vertical_speed until depth is driftdepth\"\"\"\n",
-    "        p.ddepth += vertical_speed * dt\n",
-    "        p.cycle_phase = np.where(p.depth + p.ddepth >= driftdepth, 1, p.cycle_phase)\n",
-    "        p.ddepth = np.where(\n",
-    "            p.depth + p.ddepth >= driftdepth, driftdepth - p.depth, p.ddepth\n",
-    "        )\n",
+    "        p.dz += vertical_speed * dt\n",
+    "        p.cycle_phase = np.where(p.z + p.dz >= driftdepth, 1, p.cycle_phase)\n",
+    "        p.dz = np.where(p.z + p.dz >= driftdepth, driftdepth - p.z, p.dz)\n",
     "\n",
     "    SinkingPhase(particles[particles.cycle_phase == 0])\n",
     "\n",
@@ -67,11 +65,9 @@
     "\n",
     "    def SecondSinkingPhase(p):\n",
     "        \"\"\"Phase 2: Sinking further to maxdepth\"\"\"\n",
-    "        p.ddepth += vertical_speed * dt\n",
-    "        p.cycle_phase = np.where(p.depth + p.ddepth >= maxdepth, 3, p.cycle_phase)\n",
-    "        p.ddepth = np.where(\n",
-    "            p.depth + p.ddepth >= maxdepth, maxdepth - p.depth, p.ddepth\n",
-    "        )\n",
+    "        p.dz += vertical_speed * dt\n",
+    "        p.cycle_phase = np.where(p.z + p.dz >= maxdepth, 3, p.cycle_phase)\n",
+    "        p.dz = np.where(p.z + p.dz >= maxdepth, maxdepth - p.z, p.dz)\n",
     "\n",
     "    SecondSinkingPhase(particles[particles.cycle_phase == 2])\n",
     "\n",
@@ -81,15 +77,13 @@
     "\n",
     "    def RisingPhase(p):\n",
     "        \"\"\"Phase 3: Rising with vertical_speed until at surface\"\"\"\n",
-    "        p.ddepth -= vertical_speed * dt\n",
-    "        p.temp = fieldset.thetao[p.time, p.depth, p.lat, p.lon]\n",
-    "        p.cycle_phase = np.where(\n",
-    "            p.depth + p.ddepth <= fieldset.mindepth, 4, p.cycle_phase\n",
-    "        )\n",
-    "        p.ddepth = np.where(\n",
-    "            p.depth + p.ddepth <= fieldset.mindepth,\n",
-    "            fieldset.mindepth - p.depth,\n",
-    "            p.ddepth,\n",
+    "        p.dz -= vertical_speed * dt\n",
+    "        p.temp = fieldset.thetao[p.time, p.z, p.lat, p.lon]\n",
+    "        p.cycle_phase = np.where(p.z + p.dz <= fieldset.mindepth, 4, p.cycle_phase)\n",
+    "        p.dz = np.where(\n",
+    "            p.z + p.dz <= fieldset.mindepth,\n",
+    "            fieldset.mindepth - p.z,\n",
+    "            p.dz,\n",
     "        )\n",
     "\n",
     "    RisingPhase(particles[particles.cycle_phase == 3])\n",
@@ -163,7 +157,7 @@
     "    pclass=ArgoParticle,\n",
     "    lon=[32],\n",
     "    lat=[-31],\n",
-    "    depth=[fieldset.mindepth],\n",
+    "    z=[fieldset.mindepth],\n",
     ")\n",
     "\n",
     "# combine Argo vertical movement kernel with built-in Advection kernel\n",
@@ -209,9 +203,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds_out = xr.open_zarr(\n",
-    "    output_file.store, decode_times=False\n",
-    ")  # TODO fix without using decode_times=False\n",
+    "ds_out = xr.open_zarr(output_file.store)\n",
     "x = ds_out[\"lon\"][:].squeeze()\n",
     "y = ds_out[\"lat\"][:].squeeze()\n",
     "z = ds_out[\"z\"][:].squeeze()\n",

docs/examples/tutorial_sampling.ipynb

@@ -96,7 +96,7 @@
     "time = np.repeat(\n",
     "    ds.time.values[0], npart\n",
     ")  # release all particles at the start time of the fieldset\n",
-    "depth = np.repeat(ds.depth.values[0], npart)\n",
+    "z = np.repeat(ds.depth.values[0], npart)\n",
     "\n",
     "# Plot temperature field and initial particle locations\n",
     "plt.figure()\n",
@@ -135,7 +135,7 @@
     "\n",
     "def SampleT(particles, fieldset):\n",
     "    particles.temperature = fieldset.thetao[\n",
-    "        particles.time, particles.depth, particles.lat, particles.lon\n",
+    "        particles.time, particles.z, particles.lat, particles.lon\n",
     "    ]"
    ]
   },
@@ -154,7 +154,7 @@
    "outputs": [],
    "source": [
     "pset = parcels.ParticleSet(\n",
-    "    fieldset=fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=time, depth=depth\n",
+    "    fieldset=fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=time, z=z\n",
     ")\n",
     "\n",
     "output_file = parcels.ParticleFile(\"SampleTemp.zarr\", outputdt=timedelta(hours=1))\n",
@@ -181,7 +181,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Particle_data = xr.open_zarr(\"SampleTemp.zarr\", decode_times=False)\n",
+    "Particle_data = xr.open_zarr(\"SampleTemp.zarr\")\n",
     "\n",
     "plt.figure()\n",
     "ax = plt.axes()\n",
@@ -228,7 +228,7 @@
     "\n",
     "\n",
     "pset = parcels.ParticleSet(\n",
-    "    fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, depth=depth\n",
+    "    fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z\n",
     ")\n",
     "\n",
     "pset.execute(\n",
@@ -257,7 +257,7 @@
     "\n",
     "\n",
     "pset = parcels.ParticleSet(\n",
-    "    fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, depth=depth\n",
+    "    fieldset=fieldset, pclass=parcels.Particle, lon=lon, lat=lat, time=time, z=z\n",
     ")\n",
     "\n",
     "pset.execute(SampleVel_correct, runtime=timedelta(hours=30), dt=timedelta(minutes=5))"
@@ -331,7 +331,7 @@
     "    lon=lon,\n",
     "    lat=lat,\n",
     "    time=time,\n",
-    "    depth=depth,\n",
+    "    z=z,\n",
     ")"
    ]
   },
@@ -406,7 +406,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.11"
+   "version": "3.11.0"
   },
   "pycharm": {
    "stem_cell": {

src/parcels/_core/particlefile.py

@@ -32,7 +32,7 @@
     np.dtype(np.int8): np.iinfo(np.int8).max,
     np.dtype(np.int16): np.iinfo(np.int16).max,
     np.dtype(np.int32): np.iinfo(np.int32).max,
-    np.dtype(np.int64): np.iinfo(np.int64).max,
+    np.dtype(np.int64): np.iinfo(np.int64).min,
     np.dtype(np.uint8): np.iinfo(np.uint8).max,
     np.dtype(np.uint16): np.iinfo(np.uint16).max,
     np.dtype(np.uint32): np.iinfo(np.uint32).max,
@@ -335,7 +335,7 @@ def _maybe_convert_time_dtype(dtype: np.dtype | _SAME_AS_FIELDSET_TIME_INTERVAL)
     """Convert the dtype of time to float64 if it is not already."""
     if dtype is _SAME_AS_FIELDSET_TIME_INTERVAL.VALUE:
         return np.dtype(
-            np.uint64
+            np.int64
         )  #! We need to have here some proper mechanism for converting particle data to the data that is to be output to zarr (namely the time needs to be converted to float seconds by subtracting the time_interval.left)
     return dtype
 

tests/test_particlefile.py

@@ -30,14 +30,14 @@ def fieldset() -> FieldSet:  # TODO v4: Move into a `conftest.py` file and remov
     )
 
 
-@pytest.mark.skip
 def test_metadata(fieldset, tmp_zarrfile):
     pset = ParticleSet(fieldset, pclass=Particle, lon=0, lat=0)
 
-    pset.execute(DoNothing, runtime=1, output_file=ParticleFile(tmp_zarrfile, outputdt=np.timedelta64(1, "s")))
+    ofile = ParticleFile(tmp_zarrfile, outputdt=np.timedelta64(1, "s"))
+    pset.execute(DoNothing, runtime=np.timedelta64(1, "s"), dt=np.timedelta64(1, "s"), output_file=ofile)
 
-    ds = xr.open_zarr(tmp_zarrfile, decode_cf=False)  # TODO v4: Fix metadata and re-enable decode_cf
-    assert ds.attrs["parcels_kernels"].lower() == "ParticleDoNothing".lower()
+    ds = xr.open_zarr(tmp_zarrfile)
+    assert ds.attrs["parcels_kernels"].lower() == "DoNothing".lower()
 
 
 def test_pfile_array_write_zarr_memorystore(fieldset):
@@ -76,11 +76,8 @@ def test_pfile_array_remove_particles(fieldset, tmp_zarrfile):
     pset._data["time"][:] = new_time
     pset._data["time_nextloop"][:] = new_time
     pfile.write(pset, new_time)
-    ds = xr.open_zarr(tmp_zarrfile, decode_cf=False)
+    ds = xr.open_zarr(tmp_zarrfile)
     timearr = ds["time"][:]
-    pytest.skip(
-        "TODO v4: Set decode_cf=True, which will mean that missing values get decoded to NaT rather than fill value"
-    )
     assert (np.isnat(timearr[3, 1])) and (np.isfinite(timearr[3, 0]))
 
 
@@ -129,7 +126,7 @@ def Update_lon(particles, fieldset):  # pragma: no cover
         output_file=ofile,
     )
 
-    ds = xr.open_zarr(tmp_zarrfile, decode_cf=False)  # TODO v4: Fix metadata and re-enable decode_cf
+    ds = xr.open_zarr(tmp_zarrfile)
     lons = ds["lon"][:]
     assert isinstance(lons.values[0, 0], np.float64)