interpolation_example.py

import numpy as np
from scipy.interpolate import interp1d

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def interpnjck(q, qi, a, method="linear"):

    # Warning:
    # https://stackoverflow.com/questions/45429831/valueerror-a-value-in-x-new-is-above-the-interpolation-range-what-other-re

    nd = len(q)

    assert len(qi) == nd

    for j in range(nd):
        a = interp1d(q[j], a, axis=j, kind=method,fill_value="extrapolate")(qi[j])

    return a

def main():

    points_in = np.load("points_in.npy")
    points_out = np.load("points_out.npy")

    values = np.load("orbital.npy")

    print values.shape
    print points_in.T.shape

    X, Y, Z = points_in.T
    xu = np.unique(X)
    yu = np.unique(Y)
    zu = np.unique(Z)

    del X
    del Y
    del Z
    del points_in

    print xu.shape, yu.shape, zu.shape

    X, Y, Z = points_out.T
    xuo = np.unique(X)
    yuo = np.unique(Y)
    zuo = np.unique(Z)

    print xuo.shape, yuo.shape, zuo.shape
    print "max", np.max(xu), np.max(xuo)
    print "min", np.min(xu), np.min(xuo)

    plt.plot(xu)
    plt.plot(xuo)
    plt.savefig("test_x")

    new_values = interpnjck([xu, yu, zu], [xuo, yuo, zuo], values)

 if __name__ == '__main__':
    main()