How to deal with out-of-array indexing in the interpolation functions?

[erikvansebille]

With the support of user-provided interpolation methods, we should also think how best to capture/deal with out-of-array indexing in these methods. While OutOfBounds Errors (both space and time) are thrown in search_indices (so before the interpolation function is evaluated) and we can thus be confident that data[ti, zi, yi, xi] is a valid indexing, this is not necessarily true for data[ti+1, zi+1, yi+1, xi+1] or any combination of index and index+1.

The is now captured explicitly in #1895:

Parcels/parcels/_interpolation.py

Lines 142 to 147 in 4502c12

	def _linear_2d(ctx: InterpolationContext2D) -> float:
	ft0 = _interp_on_unit_square(eta=ctx.eta, xsi=ctx.xsi, data=ctx.data[ctx.ti, :, :], yi=ctx.yi, xi=ctx.xi)
	if ctx.tau < EPS or ctx.ti >= ctx.data.shape[0] - 1:
	return ft0
	ft1 = _interp_on_unit_square(eta=ctx.eta, xsi=ctx.xsi, data=ctx.data[ctx.ti + 1, :, :], yi=ctx.yi, xi=ctx.xi)
	return (1 - ctx.tau) * ft0 + ctx.tau * ft1

But it would be nicer to fix this under the hood. Question is: how??

[fluidnumerics-joe]

Maybe we could bake in the concept of "boundary elements" and "interior elements" and not apply interpolation to boundary elements - instead, we'd need "extrapolation" to handle boundaries. Interpolation would then only apply to interior elements.

With unstructured grids, we have the boundary_face_indices that is quite useful. Although in unstructured we only have the face id and the connectivity tables, I suspect some folks may be needing to distinguish between interior and boundary points..

[erikvansebille]

Thanks for this insightful feedback, @fluidnumerics-joe. Could indeed be a good approach to differentiate between boundary elements and interior elements (although the last available time-step would also be a "boundary-in-time" element). How would we implement this in structured grids, @VeckoTheGecko?

[VeckoTheGecko]

I think the suggestion by @fluidnumerics-joe is good, though I'm not entirely sure what exact the implementation would look like (partly since this is closely related to how we handle spatial periodicity, and perhaps even temporal periodicity - though the latter we'll likely implement with xarray).