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.