parcels.FieldSet.from_croco() has issues with Fields that have multiple times

[andrew-s28]

Parcels version

3.1.1

Description

There seems to be an issue with the parcels.FieldSet.from_croco() function when using it with a file with multiple times. This first occurred when trying to run Parcels on regional CROCO fields that I am using from a collaborator (note I am not a CROCO expert by any means). Of course, this could have just been something funky with my model files, so I investigated a bit further, and can report what I think is an issue with the way Parcels is handling the Cs_w parameter when loading fields with multiple times.

The example runs just fine, but this example uses a file with a single time and allow_time_extrapolation=True. I then used the example file and created a new file with 10 times, each with the exact same data as the original file, which one would think would provide the same results. However, this results in the error shown below. This was all done in a fresh environment with v3.1.1 of Parcels, following the basic installation in the docs. Worth noting that a similar error occurs in the Field.py file if using multiple times and multiple files (e.g., one file per time), but it is a different error, and that deferred_load=True creates a slightly different error as well, but I believe they are all related. Happy to provide minimal example for these too (but it's very similar to below).

I have taken some stabs at fixing it locally. One attempt that seemed to make progress was adding an additional condition in parcels.Field.computeTimeChunk() that attempts to account for the fact that Cs_w only has dimensions of depth and reshaping it to match other fields:

# field.py
# lines 1898-1899
if len(buffer_data.shape) == 1 and g.zdim > 1:
    buffer_data = lib.reshape(buffer_data, sum(((1,), buffer_data.shape, (1, 1,)), ()))

but it just seemed to result in another error cropping up down the road, so I got to the point that I thought it was worth some discussion here.

One question before diving any deeper is: has anyone ran this function with files containing multiple times with success?

Code sample

import os

import numpy as np
import parcels
import xarray as xr

# create and save new dataset with multiple times
example_dataset_folder = parcels.download_example_dataset("CROCOidealized_data")
file = os.path.join(example_dataset_folder, "CROCO_idealized.nc")
ds = xr.open_dataset(file)

time = ds.time.values
time = [time[0]*x + time[0] for x in np.arange(10)]
ds_list = []
for t in time:
    ds_copy = ds.copy()
    ds_copy["time"] = [t]
    ds_list.append(ds_copy)
xr.concat(ds_list, dim="time").to_netcdf("./CROCO_idealized-extended.nc")

# rerun example with new dataset
file = "./CROCO_idealized-extended.nc"

variables = {"U": "u", "V": "v", "W": "w", "H": "h", "Zeta": "zeta", "Cs_w": "Cs_w"}

lon_rho = "x_rho"  # Note, this would be "lon_rho" for a dataset on a spherical grid
lat_rho = "y_rho"  # Note ,this would be "lat_rho" for a dataset on a spherical grid

dimensions = {
    "U": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "V": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "W": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "H": {"lon": lon_rho, "lat": lat_rho},
    "Zeta": {"lon": lon_rho, "lat": lat_rho, "time": "time"},
    "Cs_w": {"depth": "s_w"},
}

fieldset = parcels.FieldSet.from_croco(
    file,
    variables,
    dimensions,
    hc=xr.open_dataset(
        file
    ).hc.values,
    allow_time_extrapolation=False,  # same error happens with True or False
    mesh="flat", 
)

X, Z, T = np.meshgrid(
    [40e3, 80e3, 120e3],
    [100, -10, -130, -250, -400, -850, -1400, -1550],
    [19200],   # added time to release particles as well
)

Y = np.ones(X.size) * 10e3

def DeleteParticle(particle, fieldset, time):
    if particle.state >= 50:
        particle.delete()


pset = parcels.ParticleSet(
    fieldset=fieldset, pclass=parcels.JITParticle, lon=X, lat=Y, depth=Z, time=T
)

outputfile = pset.ParticleFile(name="croco_particles3D.zarr", outputdt=5000)

pset.execute(
    [parcels.AdvectionRK4_3D, DeleteParticle],
    runtime=5e4,
    dt=100,
    output_file=outputfile,
)

Traceback (most recent call last):
  File "~/code/test.py", line 38, in <module>
    fieldset = parcels.FieldSet.from_croco(
        file,
    ...<6 lines>...
        mesh="flat",
    )
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/fieldset.py", line 779, in from_croco
    fieldset = cls.from_netcdf(
        filenames,
    ...<9 lines>...
        **kwargs,
    )
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/fieldset.py", line 528, in from_netcdf
    fields[var] = Field.from_netcdf(
                  ~~~~~~~~~~~~~~~~~^
        paths,
        ^^^^^^
    ...<12 lines>...
        **kwargs,
        ^^^^^^^^^
    )
    ^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/field.py", line 760, in from_netcdf
    data_list.append(buffer_data.reshape(sum(((len(tslice), 1), buffer_data.shape[1:]), ())))
                              ~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: cannot reshape array of size 31110 into shape (1,1,51,61)

[VeckoTheGecko]

Thanks for the minimal example @andrew-s28 !!

Stepping through the script, the error is being encountered at the reading in of the "H" field. @erikvansebille any immediate ideas as to why that would be the case? I'm not as familiar with this part of the code

[erikvansebille]

Hi @andrew-s28, thanks indeed for providing this nice minimal example!

I dug around, and it seems the error you're getting above is because the h field now also has gotten a time dimensions; so that error disappears when you update the dimensions["H"] = {"lon": lon_rho, "lat": lat_rho, "time": "time"},

Another error then pops up though, about Cs_w also being the wrong shape; I expect again it now has a time dimensions. That would not be normal, right? I understood that Cs_w was time-invariant in CROCO?

Could you update your code so that the CROCO_idealized-extended.nc file only extends the variables that do vary with time (u, v, w and zeta?)?

[andrew-s28]

Thanks for your reply @VeckoTheGecko and @erikvansebille.

Indeed, both H and Cs_w should not have dimensions of time. Fixing the example I gave with xr.concat(ds_list, dim="time", data_vars="minimal").to_netcdf("./CROCO_idealized-extended.nc") (adding the data_vars="minimal" kwarg) does indeed make the run successful, in the case of a single input file. Somehow in troubleshooting my original use case, I added a second bug along the way. Sorry about that oversight - that's what I get for trying to figure this out late on a Friday.

Back to that original use case though, I mentioned in the issue that when using multiple files instead (since the model I am using is split over ~40 files, each with 72 timestamps), I was getting a different error. See the code block below for this minimal example and error outputs (I include three error outputs, one for JITParticle, one for ScipyParticle, and one with deferred_load=False, since they are different but seem related to a common dimensions issue). I think/hope I've done my due diligence this time around and this isn't just another issue with the file creation, but of course let me know if I'm missing something obvious.

Happy to close this issue as completed and open another one or move to the discussions instead if this feels like a separate topic, but this was the original motivation for the issue and is closely related.

Updated Code Sample

Code

import glob
import os

import numpy as np
import parcels
import xarray as xr

# create and save new dataset with multiple times
example_dataset_folder = parcels.download_example_dataset("CROCOidealized_data")
file = os.path.join(example_dataset_folder, "CROCO_idealized.nc")
ds = xr.open_dataset(file)

time = ds.time.values
time = [time[0]*x + time[0] for x in np.arange(10)]
ds_list = []
for t in time:
    ds_copy = ds.copy()
    ds_copy["time"] = [t]
    ds_copy.to_netcdf('CROCO_idealized-'+str(t)+'.nc')  # this time making 10 separate files, each with 1 timestep

# rerun example with new dataset
files = glob.glob("./CROCO_idealized-*.0.nc")

variables = {"U": "u", "V": "v", "W": "w", "H": "h", "Zeta": "zeta", "Cs_w": "Cs_w"}

lon_rho = "x_rho"  # Note, this would be "lon_rho" for a dataset on a spherical grid
lat_rho = "y_rho"  # Note ,this would be "lat_rho" for a dataset on a spherical grid

dimensions = {
    "U": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "V": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "W": {"lon": lon_rho, "lat": lat_rho, "depth": "s_w", "time": "time"},
    "H": {"lon": lon_rho, "lat": lat_rho},
    "Zeta": {"lon": lon_rho, "lat": lat_rho, "time": "time"},
    "Cs_w": {"depth": "s_w"},
}

fieldset = parcels.FieldSet.from_croco(
    files,
    variables,
    dimensions,
    hc=xr.open_dataset(
        files[0]
    ).hc.values,
    allow_time_extrapolation=False,  # same error happens with True or False
    mesh="flat",
    timestamps=[[t] for t in time],  # requires manually setting timestamps or another error pops up
)

X, Z, T = np.meshgrid(
    [40e3, 80e3, 120e3],
    [100, -10, -130, -250, -400, -850, -1400, -1550],
    [19200],   # added time to release particles as well
)

Y = np.ones(X.size) * 10e3

def DeleteParticle(particle, fieldset, time):
    if particle.state >= 50:
        particle.delete()


pset = parcels.ParticleSet(
    fieldset=fieldset, pclass=parcels.JITParticle, lon=X, lat=Y, depth=Z, time=T
)

outputfile = pset.ParticleFile(name="croco_particles3D.zarr", outputdt=5000)

pset.execute(
    [parcels.AdvectionRK4_3D, DeleteParticle],
    runtime=5e4,
    dt=100,
    output_file=outputfile,
)

Error Output

JITParticle

Traceback (most recent call last):
  File "/mnt/d/code/parcels/test.py", line 69, in <module>
    pset.execute(
    ~~~~~~~~~~~~^
        [parcels.AdvectionRK4_3D, DeleteParticle],
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    ...<2 lines>...
        output_file=outputfile,
        ^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/particleset.py", line 1046, in execute
    self._kernel = self.Kernel(pyfunc, delete_cfiles=delete_cfiles)
                   ~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/particleset.py", line 889, in Kernel
    return Kernel.from_list(
           ~~~~~~~~~~~~~~~~^
        self.fieldset,
        ^^^^^^^^^^^^^^
    ...<3 lines>...
        delete_cfiles=delete_cfiles,
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/kernel.py", line 517, in from_list
    pyfunc_list[0] = cls(fieldset, ptype, pyfunc_list[0], *args, **kwargs)
                     ~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/kernel.py", line 254, in __init__
    kernel_ccode = kernelgen.generate(deepcopy(self.py_ast), self.funcvars)
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/compilation/codegenerator.py", line 445, in generate
    self.visit(py_ast)
    ~~~~~~~~~~^^^^^^^^
  File "~/miniconda3/envs/parcels3/lib/python3.13/ast.py", line 428, in visit
    return visitor(node)
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/compilation/codegenerator.py", line 491, in visit_FunctionDef
    self.visit(stmt)
    ~~~~~~~~~~^^^^^^
  File "~/miniconda3/envs/parcels3/lib/python3.13/ast.py", line 428, in visit
    return visitor(node)
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/compilation/codegenerator.py", line 862, in visit_VectorFieldEvalNode
    Cs_w = [float(self.fieldset.Cs_w.data[0][zi][0][0]) for zi in range(self.fieldset.Cs_w.data.shape[1])]
                                                                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'DeferredArray' object has no attribute 'shape'

ScipyParticle

Traceback (most recent call last):
  File "/mnt/d/code/parcels/test.py", line 69, in <module>
    pset.execute(
    ~~~~~~~~~~~~^
        [parcels.AdvectionRK4_3D, DeleteParticle],
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    ...<2 lines>...
        output_file=outputfile,
        ^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/particleset.py", line 1161, in execute
    next_input = self.fieldset.computeTimeChunk(time, dt)
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/fieldset.py", line 1585, in
 computeTimeChunk
    data = f.computeTimeChunk(data, tind)
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/field.py", line 1906, in computeTimeChunk
    data = self._data_concatenate(data, buffer_data, tindex)
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/field.py", line 1848, in _data_concatenate
    data = lib.concatenate([data_to_concat, data[tindex + 1 :, :]], axis=0)
ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 1 dimension(s) and the array at index 1 has 4 dimension(s)

deferred_load=False

Traceback (most recent call last):
  File "/mnt/d/code/parcels/test.py", line 38, in <module>
    fieldset = parcels.FieldSet.from_croco(
        files,
    ...<8 lines>...
        timestamps=[[t] for t in time],  # requires manually setting timestamps or another error pops up
    )
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/fieldset.py", line 779, in from_croco
    fieldset = cls.from_netcdf(
        filenames,
    ...<9 lines>...
        **kwargs,
    )
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/fieldset.py", line 528, in from_netcdf
    fields[var] = Field.from_netcdf(
                  ~~~~~~~~~~~~~~~~~^
        paths,
        ^^^^^^
    ...<12 lines>...
        **kwargs,
        ^^^^^^^^^
    )
    ^
  File "~/miniconda3/envs/parcels3/lib/python3.13/site-packages/parcels/field.py", line 748, in from_netcdf
    assert buffer_data.shape[0] == grid.tdim, errormessage
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError: Field u expecting a data shape of [tdim=10, zdim=10, ydim=51, xdim=61] but got shape (1, 10, 51, 61). Flag transpose=True could help to reorder the data.

[erikvansebille]

Thanks for this updated code, @andrew-s28. I have dug around in the code a bit, and this patch below seems to at least fix the errors

diff --git a/parcels/field.py b/parcels/field.py
index 67ac56d4..58f233c7 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -682,14 +682,17 @@ class Field:
             raise RuntimeError("Multiple files given but no time dimension specified")
 
         if grid is None:
-            # Concatenate time variable to determine overall dimension
-            # across multiple files
-            time, time_origin, timeslices, dataFiles = cls._collect_timeslices(
-                timestamps, data_filenames, _grid_fb_class, dimensions, indices, netcdf_engine
-            )
-            grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
-            grid.timeslices = timeslices
-            kwargs["dataFiles"] = dataFiles
+            if "time" in dimensions:
+                time, time_origin, timeslices, dataFiles = cls._collect_timeslices(
+                    timestamps, data_filenames, _grid_fb_class, dimensions, indices, netcdf_engine
+                )
+                grid = Grid.create_grid(lon, lat, depth, time, time_origin=time_origin, mesh=mesh)
+                grid.timeslices = timeslices
+                kwargs["dataFiles"] = dataFiles
+            else:
+                grid = Grid.create_grid(lon, lat, depth, np.array([0.]), time_origin=TimeConverter(0.), mesh=mesh)
+                grid.timeslices = [[0]]
+                data_filenames = [data_filenames[0]]
         elif grid is not None and ("dataFiles" not in kwargs or kwargs["dataFiles"] is None):
             # ==== means: the field has a shared grid, but may have different data files, so we need to collect the
             # ==== correct file time series again.

Could you check if this patch also fixes your simulation? If so, I'll clean it up further (it's a quite hacky workaround now) and submit to main

[andrew-s28]

This does seem to be the fix - thanks for helping me get this one resolved. Perhaps including a warning of some kind if time is not in dimensions and the variable isn't H or Cs_w (or some other condition that makes sense) would be good, since if I'm understanding correctly, I think(?) what this fix does is force the single time in the field to be reused for all times if a time dimension is not present. Especially considering the previous behavior was to error.

PS: the development of the Parcels CROCO support came at a perfect time for my work, and it didn't look like an easy lift, so thanks in general for developing this! Excited to make use of it.

[erikvansebille]

Thanks @andrew-s28, for reporting back that the patch above worked.

The multi-file handling in Parcles is a bit (very?) clunky and not so flexible, which is one of the reasons we're planning to move over to full xarray-based FieldSet in Parcels v4.

In the meantime, I guess that the patch is good enough for now, so I created PR #1835 to merge this into main

And thanks for your nice comments about Parcles Croco support. Do note that this is still a bit experimental, so let us know whenever you run into issues/strange behaviour!

[andrew-s28]

Sorry to ping this issue again @erikvansebille / @VeckoTheGecko, but I noticed a small bug with the fix in #1835.

The condition in line 688 that raises an error is equivalent to the negation of the condition in the fix introduced in line 694, meaning the else statement in line 701 will never trigger since an error is raised:

https://github.com/OceanParcels/Parcels/blob/00959f6d57bd0baa617d1cb97bf4571cfa7236ff/parcels/field.py#L688-L704

I would propose to fix this by turning the error in line 689 into a warning (FieldSetWarning?), allowing the code to continue but maintaining a user notification (since the condition it was throwing an error for is now explicitly handled). I would be happy to submit a PR if this fix sounds reasonable.

[erikvansebille]

Thanks for pointing this out, @andrew-s28. I agree that in many cases the RuntimeError will be triggered before the else statement of line 701. This is code that will very likely be removed in our v4 development, so I think that in the meantime changing the RuntimeError to a warning so that you can keep using the code is a good idea