Search code examples
pythonnumpyscipyinterpolationsmoothing

scipy interpolation is not smoothing my data


I'm trying to use interpolation in scipy. Here is my code:

from Constants import LOWER_LAT, LOWER_LONG, UPPER_LAT, UPPER_LONG, GRID_RESOLUTION

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
from cmath import sin
from scipy.signal.windows import cosine
from scipy import interpolate
from scipy.interpolate import RectBivariateSpline
import numpy
from numpy import meshgrid

#===============================================================
y_range = GRID_RESOLUTION
delta = (UPPER_LAT - LOWER_LAT)/float(GRID_RESOLUTION)
x_range = int((UPPER_LONG - LOWER_LONG)/delta) + 1
x = numpy.linspace(0,x_range-1,x_range)
y = numpy.linspace(0,y_range-1,y_range)
X,Y = meshgrid(x,y)
Z = numpy.zeros((y.size, x.size))
base_val = 0

# fill values for Z
with open('map.txt','rb') as fp:
    for line in fp:
        parts = line[:-1].split("\t")
        tup = parts[0]
        tup = tup[:-1]
        tup = tup[1:]
        yx = tup.strip().replace(" ","").split(",")
        y_val = int(yx[0])
        x_val = int(yx[1])
        h_val = int(parts[-1])

        for i in range(y_range):
            tx = X[i];
            ty = Y[i];
            tz = Z[i];
            for j in range(x_range):
                if (int(tx[j])==x_val) and (int(ty[j])==y_val):
                    tz[j] = h_val + base_val
Z = numpy.array(Z)

# spline = RectBivariateSpline(y, x, Z)
# Z2 = spline(y, x)
f = interpolate.interp2d(x, y, Z,'cubic')
Z2 = f(x,y)


# Plot here
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, alpha=0.3, cmap='Accent')
ax.set_xlabel('X')
ax.set_xlim(0, 50)
ax.set_ylabel('Y')
ax.set_ylim(0, 50)
ax.set_zlabel('Z')
# ax.set_zlim(0, 1000)
plt.show()

Here are some constants from the top of the above code:

LOWER_LAT = 32.5098
LOWER_LONG = -84.7485
UPPER_LAT = 47.5617
UPPER_LONG = -69.1699
GRID_RESOLUTION = 50

My code creates 1D arrays x, y and then creates the grid with function meshgrid. Values in Z are filled from a text file that you can find here. Each line in the text file has format of (y_value,x_value) z_value. After creating the grid and interpolating the function, I plot it. However, the figure I obtain is the same as the figure I got without interpolation. Concretely, these two lines produce the same figure:

ax.plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.3, cmap='Accent')
ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, alpha=0.3, cmap='Accent')

In the lines above, Z2's values are from the interpolation function, and the Z's values are original values. How can I make the interpolation work? Here is the figure. enter image description here


Solution

  • I think you are confusing smoothing with interpolation.

    In both cases you are fitting a function that yields a continuous approximation of your input data. However, in the case of interpolation the interpolant is constrained to pass exactly through your input points, whereas with smoothing this constraint is relaxed.

    In your example above, you have performed interpolation rather than smoothing. Since you are evaluating your interpolant on the exact same grid of input points as your original data then Z2 is guaranteed to be almost exactly the same as Z. The point of doing interpolation is so that you can evaluate the approximate z-values given a different set of x- and y-values (e.g. a more finely-spaced grid).

    If you want to perform smoothing rather than interpolation, you could try passing a non-zero value as the s= argument to RectBivariateSpline, e.g.:

    spline = RectBivariateSpline(y, x, Z, s=5E7)
    Z2 = spline(y, x)
    
    fig, ax = plt.subplots(1, 2, sharex=True, sharey=True,
                           subplot_kw={'projection':'3d'})
    
    ax[0].plot_surface(X, Y, Z, rstride=1, cstride=1, alpha=0.3, cmap='Accent')
    ax[1].plot_surface(X, Y, Z2, rstride=1, cstride=1, alpha=0.3, cmap='Accent')
    ax[0].set_title('Original')
    ax[1].set_title('Smoothed')
    
    fig.tight_layout()
    plt.show()
    

    enter image description here