Gridded value at specific coordinate

General Usage
1

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. . . .

2

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
3

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.

4

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
5

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