# Gridded value at specific coordinate

**URL:** https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427
**Category:** General Usage
**Created:** [May 19, 2021, 9:59am UTC](https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427 "2021-05-19T09:59:02Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)
#### Post date: [May 19, 2021, 9:59am UTC](https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427/1 "2021-05-19T09:59:02Z")

</div>

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

---

<div class="post-metadata">

### Author: ![hendri54](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hendri54/32/9621_2.png) [@hendri54](https://discourse.julialang.org/u/hendri54)
#### Post date: [May 19, 2021, 12:39pm UTC](https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427/2 "2021-05-19T12:39:29Z")

</div>

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

> **[General Purpose Packages](https://julia.quantecon.org/more_julia/general_packages.html#Interpolation)**
>
> This website presents a set of lectures on quantitative economic modeling, designed and written by Jesse Perla, Thomas J. Sargent and John Stachurski. The language instruction is Julia.

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

```julia
# 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);

```

---

<div class="post-metadata">

### Author: ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)
#### Post date: [May 19, 2021, 1:16pm UTC](https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427/3 "2021-05-19T13:16:42Z")

</div>

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

```julia
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:

```julia
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:

```julia
 # 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.

---

<div class="post-metadata">

### Author: ![hendri54](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hendri54/32/9621_2.png) [@hendri54](https://discourse.julialang.org/u/hendri54)
#### Post date: [May 22, 2021, 12:08pm UTC](https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427/4 "2021-05-22T12:08:00Z")

</div>

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):

```julia
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)

```

---

<div class="post-metadata">

### Author: ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [May 22, 2021, 12:20pm UTC](https://discourse.julialang.org/t/gridded-value-at-specific-coordinate/61427/5 "2021-05-22T12:20:10Z")

</div>

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

```julia
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](https://github.com/cesaraustralia/TimeInterpolatedGeoData.jl) but interpolating whole layers with a `sin` interpolator over large lazy-loaded datasets. Maybe a bit too complicated as an example.
