Nearest Neighbor Interpolation

[simone-silvestri]

This PR introduces an NN interpolation routine to interpolate data from a TRG to a latitude-longitude grid. It is quite crude at the moment and not reliable near the fold, but it looks like it performs ok for the rest of the domain.

I am now using it only for visualization and to perform zonal averages, I am not even sure we want to use this algorithm for actual interpolation in Oceananigans.

Anywas, since I am using it, I ll leave it here so we can decide what to do about it

[glwagner]

What's wrong with the interpolation we currently have?

[simone-silvestri]

It uses a sorted binary search on a 1D array to find the index corresponding to the closest point. In this case, the coordinates are 2D and one of them is not sorted (the longitude)

[glwagner]

It uses a sorted binary search on a 1D array to find the index corresponding to the closest point. In this case, the coordinates are 2D and one of them is not sorted (the longitude)

Certainly this bug should be documented in an issue...

[simone-silvestri]

I don't think we ever tackled the problem of interpolation from an OrthogonalSphericalShellGrid, however we probably should throw an instructive error when trying to interpolate from a grid with 2D coordinates.

[glwagner]

Just show us the code that has the problem

[navidcy]

Suggested change

	ifelse(λ₀ > 180, λ + 360, λ - 360), λ)
	ifelse(λ₀ > 180, λ + 360, λ - 360), λ)

[simone-silvestri]

Just show us the code that has the problem

I wouldn't necessarily say that the code has a problem, just that it does not support interpolation on a "semi-unstructured" grid. We don't have an interpolation method for an OrthogonalSphericalShellGrid in the interpolate.jl file, just a linear binary search that is perfect for 1D coordinates. We haven't yet implemented a method for OrthogonalSphericalShellGrids, so since I needed it, I thought I would implement the simplest version. I am not sure this simple interpolation can cut it, but if we want, we can add it to Oceananigans as a first pass for interpolation from an OrthogonalSphericalShellGrid

[simone-silvestri]

Probably, this algorithm wouldn't work for a general OrthogonalSphericalGrid though because, to be efficient, it assumes that the φ coordinate is monotonic in j, which is true for the tripolar grid, but it is not the case in general, for example, the latitude on the north face of a conformal cubed sphere is not monotonic in j. A general algorithm might need a full 2D unsorted search or a coordinate transformation. We can discard this PR and focus on developing a general regridding method in Oceananigans. I just thought it might be useful to leave this algorithm here since it seems easy enough.

[glwagner]

show the code that doesn't work

[simone-silvestri]

https://github.com/CliMA/Oceananigans.jl/blob/93cf3e9fee3c9c479641f10f99207c83b6538984/src/Fields/interpolate.jl#L65-L129

[glwagner]

No I mean can you produce a code snippet that will fail

like grid = , and then field = , and then interpolate!(field1, field2) and then shows there is either an error or the result is wrong

[simone-silvestri]

Ah ok in that sense!

using Oceananigans
using OrthogonalSphericalShellGrids
using Oceananigans.Fields: interpolate!

trg = TripolarGrid(size = (10, 10, 10), z = (0, 1))
llg = LatitudeLongitudeGrid(size = (10, 10, 10), latitude = (-75, 75), longitude = (0, 360), z = (0, 1))

ctrg = CenterField(trg)
cllg = CenterField(llg)

interpolate!(cllg, ctrg)

[glwagner]

hallelujah