Numpy scipy 2d interpolation for linear piecewise data

Question

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.

Answer 1

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)

Answer 2

This is a long one: You could write a function to accomplish the same:

import numpy as np
def interpolate(p, points, values):
    xpoint, ypoint = np.unique(points[:,0]), np.unique(points[:,1])
    z = values.reshape(ypoint.size, xpoint.size).T
    x, y = p[:,0], p[:,1]
    xind = np.searchsorted(xpoint, x)
    yind = np.searchsorted(ypoint, y)
    xind_min = np.maximum(xind - 1, 0)
    yind_min = np.maximum(yind - 1, 0)
    Q11 = z[xind_min, yind_min]
    Q12 = z[xind_min, yind]
    Q21 = z[xind, yind_min]
    Q22 = z[xind, yind]
    x1, x2 = xpoint[xind_min], xpoint[xind]
    y1, y2 = ypoint[yind_min], ypoint[yind]
    x_diff = np.where(x2 > x1, x2 - x1, 1)
    xy1 = ((x2 - x)*Q11 + (x - x1)*Q21) / x_diff
    xy2 = ((x2 - x)*Q12 + (x - x1)*Q22) / x_diff
    y_diff = np.where(y2 > y1, y2 - y1, 1)
    res = ((y2 - y) * xy1 + (y - y1) * xy2)/y_diff 
    equalX, equalY = x2 == x, y2 == y
    equalXY = equalX & equalY
    res[equalXY] =  z[xind[equalXY], yind[equalXY]]
    res[equalX & ~equalY] = z.mean(1)[xind[equalX & ~equalY]]
    res[~equalX & equalY] = z.mean(0)[yind[~equalX & equalY]]
    return res


points = np.array([[0, 105],[5000, 105],[0, 135],[5000, 135],[0, 165],[5000, 165]])
values = np.array([[300, 380, 300, 390, 300, 400]]).transpose()
xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])

interpolate(xi, points, values)
array([340.        , 343.33333333, 345.        , 347.5       ,
       350.        ])

Answer 3

i run the code and get the same results as you.

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import numpy as np
from scipy.interpolate import griddata

points = np.array([[0, 105], [5000, 105], [0, 135], [5000, 135], [0, 165], [5000, 165]])
values = np.array([[300, 380, 300, 390, 300, 400]]).transpose()
xi = np.array([[2500, 105],[2500, 125],[2500, 135],[2500, 150],[2500, 165]])

gd = griddata(points, values, xi, method='linear', rescale=True)
print(gd)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the surface
ax.plot_trisurf(points[:,0], points[:,1], values.flatten(), cmap=cm.get_cmap('Blues'), linewidth=0.2)

# Customize the axes
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.xaxis.set_major_locator(LinearLocator(3))
ax.yaxis.set_major_locator(LinearLocator(3))
ax.zaxis.set_major_locator(LinearLocator(5))

plt.show()

the results are this:

[[340.]
 [345.]
 [345.]
 [345.]
 [350.]]

if i change from linear to cubic i get this:

[[340.46073832]
 [343.60340822]
 [345.00000255]
 [347.06944598]
 [349.53926443]]

which is much closer to what you expect... so it looks like an interpolation method issue.

it looks like you are interpolating across this surface:

if you are able to show how you arrive at the expected result then it would be helpful...