Add checking that Fields have expected number of dimensions when calling an interpolation method

[erikvansebille]

What version of Parcels are you running?

v4alpha

Is your feature request related to a problem?

Due to a bug in convert.nemo_to_sgrid(), the MOi-benchmark resulted in a DataSet with dimensions (depth: 1, time: 3, x: 4322, y_center: 3059, depth_center: 50, x_center: 4322, y: 3059) and a variable U with dimensions (depth, time, depth_center, y_center, x) (so five dimensions, instead of four)

This then breaks e.g. the interpolator code that expects exactly four dimensions (and in order time, depth, lat, lon) at

Parcels/src/parcels/interpolators/_xinterpolators.py

Line 186 in 2012123

tdim, zdim, ydim, xdim = U.shape[0], U.shape[1], U.shape[2], U.shape[3]

Describe the solution you'd like

It would be good to put some guarding in place so that Fields have the expected number of dimensions when passed to an interpolator. Perhaps this can be done with typing?

We could expect users to write interpolators for other numbers of dimensions; so I wouldn't hard-code a check on four dimensions. But a bit more guarding would be useful

Also, perhaps the code in _xinterpolators.py above could be a bit more robust against the order of the dimensions?

Describe alternatives you've considered

No response

Additional context

No response

[VeckoTheGecko]

I think that this touches on another problem as well - which assumptions can those writing interpolators make about the fieldsets that they're working work? This mentions the dimensionality, but there's other things (riffing off the top of my head):

• Dimensionality in the case of 2D etc fields
• Types of the underlying arrays (dask, numpy, cupy etc)
• Units

Since you've written a bunch of interpolators already, what do you think about this @erikvansebille ? I don't think we have to have all the answers now, but if we don't - perhaps its worth setting up a tracking issue so that we can work towards this

[erikvansebille]

Yes, all of these are good points. And a fourth one is land-value (see also #2451).

As for Dimensionality, I thought that we always assumed fields to be 4D (time, depth, lat, lon) in Parcels? So writers of interpolators can be confident that the order is that way?

As for the array-type; perhaps we can declare a generalised array operation, that then links to the same array package that the rest of Parcels?

As for units, that is up to the user - until we have something like support for pint; but that is for a later version...)

[VeckoTheGecko]

Just made #2627 to continue discussion