x1, y1 = 1.3, 25000
x2, y2 = 0.2, 50000
x3, y3 = 2.5, 15000
xi = np.array([x1, x2, x3], dtype=np.float)
yi = np.array([y1, y2, y3], dtype=np.float)
xi[xi > xmax] = xmax
xi[xi < xmin] = xmin
yi[yi < ymin] = ymin
xi = (nx - 1) * (xi - xmin) / (xmax - xmin)
yi = (ny - 2) * (yi - ymin) / (ymax - ymin)
yi[yi > 1] = 2.0
z1, z2, z3 = map_coordinates(z, [yi, xi])
for X, Y, Z in zip((x1, x2, x3), (y1, y2, y3), (z1, z2, z3)):
    pass
Scipy interpolation on a numpy array