Numpy scipy 2d interpolation for linear piecewise data

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.

Solution

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)