I have the points:
points = np.array([[0, 105],[5000, 105],[0, 135],[5000, 135],[0, 165],[5000, 165]])
and values
values = np.array([[300, 380, 300, 390, 300, 400]]).transpose()
with the inputs I'm trying to interpolate for
xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])
with expected result for bilinear interpolation (ref: https://en.wikipedia.org/wiki/Bilinear_interpolation)
[340, 343.3333, 345, 347.5, 350]
My working for the second example using bilinear interpolation
x1=2500, y1=105 giving z1=340
x2=2500, y2=135 giving z2=345
Hence for x3=2500, y3=125 gives z3=343.3333
however, with
gd = griddata(points, values, xi, method='linear', rescale=True)
I'm getting the result
[340, 345, 345, 345, 350]
I must be missing something simple here, but have gotten nowhere trying multiple different approaches.
This can be done using scipy.interpolate.interpn
if you provide the data correctly. It expects you to provide the points as a list of individual x and y values (for the 2D case) that define the grid. The values are then defined in the format that corresponds to the grid, which is the result of np.meshgrid(x, y, indexing="ij")
for the 2D case. Be sure to provide x and y in strictly ascending or descending order or interpn
will throw an error.
import numpy as np
from scipy.interpolate import interpn
x = np.array([0, 5000])
y = np.array([105, 135, 165])
# data corresponds to (x, y) given by np.meshgrid(x, y, indexing="ij")
values = np.array([[300, 300, 300],
[380, 390, 400]])
xi = np.array([[2500, 105],
[2500, 125],
[2500, 135],
[2500, 150],
[2500, 165]])
interpolated = interpn((x, y), values, xi) # array([340., 343.33333333, 345., 347.5, 350.])
Here it is written as a function, though it doesn't have all the functionality necessary for general use, namely checking the inputs, ensuring proper datatypes, etc. I also haven't tested it beyond the given inputs, so I might have overlooked something.
import numpy as np
from scipy.interpolate import interpn
# moved one of the points and values to show it works with both unsorted
# x and y values
points = np.array([[0, 105],[5000, 105],[5000, 135],[0, 165],[5000, 165],[0, 135]])
values = np.array([[300, 380, 390, 300, 400, 300]]).T
xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])
def bilinear_interpolation(points, values, xi):
p = points.copy()
v = values.copy()
# sorting the points to ensure ascending order
ysorted = np.argsort(p[:,1])
p = p[ysorted]
v = v[ysorted]
xsorted = np.argsort(p[:,0])
p = p[xsorted]
v = v[xsorted]
x = np.unique(p[:,0])
y = np.unique(p[:,1])
v = v.reshape(x.size, y.size)
return interpn((x, y), v, xi)
interpolated = bilinear_interpolation(points, values, xi)