2015-11-28 2 views
2

J'essaye d'utiliser l'interpolation dans scipy. Voici mon code:L'interpolation scipy ne lisse pas mes données

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() 

Voici quelques constantes du haut du code ci-dessus:

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

Mon code crée des tableaux 1D x, y et crée ensuite la grille avec la fonction meshgrid. Les valeurs dans Z sont remplies à partir d'un fichier texte que vous pouvez trouver here. Chaque ligne du fichier texte a le format (y_value,x_value) z_value. Après avoir créé la grille et interpolé la fonction, je la trace. Cependant, le chiffre que j'obtiens est le même que celui que j'ai obtenu sans interpolation. Concrètement, ces deux lignes produisent le même chiffre:

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') 

Dans les lignes ci-dessus, les valeurs de Z2 sont de la fonction d'interpolation, et les valeurs du Z sont des valeurs d'origine. Comment puis-je faire fonctionner l'interpolation? Voici la figure. enter image description here

Répondre

2

Je pense que vous confondez smoothing avec interpolation.

Dans les deux cas, vous utilisez une fonction qui fournit une approximation continue de vos données d'entrée. Cependant, dans le cas de l'interpolation, l'interpolant est contraint de passer exactement par vos points d'entrée, alors qu'avec le lissage, cette contrainte est relâchée.

Dans votre exemple ci-dessus, vous avez effectué par interpolation au lieu de lisser. Puisque vous évaluez votre interpolant sur exactement la même grille de points d'entrée que vos données d'origine, alors Z2 est garanti être presque exactement le même que Z. Le point de l'interpolation est de sorte que vous pouvez évaluer les valeurs z approximatives données différents ensemble de valeurs x et y (par exemple une grille plus finement espacées).

Si vous souhaitez effectuer lissage plutôt que l'interpolation, vous pouvez essayer de faire passer une valeur non nulle comme argument s= à RectBivariateSpline, par exemple:

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

+0

Merci, exactement ce que Je veux et en effet je confond les concepts. Pouvez-vous expliquer la signification de '5E7' ou me donner une référence pour cela? – lenhhoxung

+1

C'est juste la notation scientifique pour 50000000, c'est-à-dire 5 x 10^7 –