Gridded value at specific coordinate

I’m quite new to Julia.

Suppose you have a 1D gridded data, say, temperature, `[23.5, 20.0, . . . ]`, defined at gridpoints xs = [1.0, 2.0, 5.0, . . . ]. Assume that `length(temperature) == length(xs)` and that `xs` is monotonic.

Now, I want an interpolated value of temperature at x = 3.0, say. I’m aware of the `Interpolations` package and have used it once. So, basically, I can write a program to achieve my goal. But, I will be doing this a lot, so I’d like to know

1. Is there already a library for this?
2. If not, what would be a good design strategy for this library?

If possible, I would like an interface like `temperature[x = 3.0]` to get the (interpolated) value at x = 3.0.

The above example is for a 1D variable `T(x)` but I use up to 5D variables. I would like a few choices for interpolation scheme like “nearest” and “linear”. I would need missing values be handled gracefully and reasonably.

. . . sounds like quite an undertaking. . . .

Could you explain more what you don’t like about `Interpolations`? Say we take the flow of this QuantEcon example

There we have pretty much the interface that you seem to be asking for:

``````# cubic spline interpolation
interp_cubic = CubicSplineInterpolation((xs, ys), A)
@show interp_cubic(3, 2) # exactly log(3 + 2)
@show interp_cubic(3.1, 2.1) # approximately log(3.1 + 2.1);
``````
2 Likes

I illustrate my imaginary use case below. Perhaps using a 2D gridded data is clearer as an example:

``````ds = Dataset("somedatafile.nc", "r") # netCDF file
temp = ds["temperature"]  # read 2D array
xs = ds["x_axis"] # read the x axis
ys = ds["y_axis"] # read the y axis
temp_on_grid  = GriddedData(temp, xs, ys) #<- This is what I want.
plot(ys, temp_on_grid[x = 3.0]) # plot the x=3 section of temp.
``````

So, I want to hide the interpolator behind an operation like “[x = 3.0]”.

The above is a fictitious Julia code. In real Julia, you can easily plot an i-th section:

``````plot(ys, temp[2,:])
``````

I want to extend this convenience to `x` as a continous index, if you like.

I can see that I could define such a type as GriddedData above, but before embarking on such a project, I’d like to know whether there is already such an interface, and if not, what’s the best strategy for it. For example, should I utilize Pair, as `temp( :x => 3.0 )` ? Or perhaps I should override the indexing operation to allow for non-integral indices, like `temp[3.0, :]` Or maybe the constructor should return a function `temp_on_grid(x,y)` which can be used as a regular function `temp_on_grid(3.0, undef)`. . . I don’t know. I’m not that familiar with Julia. I don’t know how to re-define (override) the “” function as in C++.

Moreover, I’d like to eventually extend this to allow for, say, integrals:

`````` # integral from x = 2.5 to 10.0
a = definite_integral(x = (2.5, 10.0), temp_on_grid)
# integral from x = 2.5 to x
f = indefinite_integral(x0 = 2.5, temp_on_grid)
``````

Again interpolation would be automatically carried out.

Sorry- being slow here. I don’t see how to make this work with names (such as `x`) for the dimensions (because the interpolation doesn’t know those). But with dimension indices, something like this works (using undocumented internals, though):

``````using Interpolations

# CubicSpline only workers with range inputs!
x = 1:10;
y = 2:6;
z = x .+ x .* y';
itp = CubicSplineInterpolation((x,y), z);

function dim_slice(itp, x, dim)
innerItp = parent(itp);
ranges = Vector{Any}(collect(innerItp.ranges));
ranges[dim] = x;
return itp(ranges...)
end

@show dim_slice(itp, 3.5, 1) == itp(3.5, y)

@show dim_slice(itp, 2.5, 2) == itp(x, 2.5)

``````
1 Like

I’ve thought about writing up named dimension interpolation for DimensionalData.jl. It would need an `Interp` selector so you could write something like:

``````A[X=Interp(0.7)]
``````

Or maybe we could use e.g. `BSpline` directly instead of `Interp`.

It’s probably only 10 lines of code to add this, but the interpolations.jl dependency is too heavy to include in DimensionalData.jl unfortunately. A glue package could work.

We do something kind of similar in GitHub - cesaraustralia/TimeInterpolatedGeoData.jl but interpolating whole layers with a `sin` interpolator over large lazy-loaded datasets. Maybe a bit too complicated as an example.

1 Like