Interpolation d'une surface avec des splines

Version imprimable

1 pièce(s) jointe(s)

Bonjour

J'essaye d'interpoler une surface par des splines. J'utilise les fonction bisplrep et bisplev mais j'ai une erreur.

Voici mon code :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
#!/usr/bin/python
# -*- coding=utf-8 -*-
 
import scipy as sp
 
nx = 29
ny = 34
 
def spline(data):
 
    #from scipy.interpolate import BivariateSpline
    from scipy.interpolate import bisplrep, bisplev
 
    # data
    x, y, z = data.transpose()
 
    m = nx
    print("s = [{0}:{1}]".format(m-sp.sqrt(2.*m), m + sp.sqrt(2.*m)))
 
    # create mgrid type object
    gx = x.reshape( nx, ny)
    gy = y.reshape( nx, ny)
    gz = z.reshape( nx, ny)
 
    # create spline object
    tck = bisplrep( gx, gy, gz)
 
    # new grid
    xn, yn = sp.mgrid[x.min():x.max():50j,y.min():y.max():50j]
    zspline = bisplev( xn[:,0], yn[0,:], tck)
 
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_wireframe( gx, gy, gz)
    ax.plot_wireframe( xn, yn, zspline)
    plt.show()
 
if __name__ == "__main__":
    # load data
    data = sp.loadtxt("resultats.dat", usecols = (0,1,3))
    spline(data)

Et voici l'erreur que j'obtient :

Code:

1
2
3
4
5
6
7
8
9
 
Traceback (most recent call last):
  File "interp2D.py", line 56, in <module>
    spline(data)
  File "interp2D.py", line 26, in spline
    tck = bisplrep( gx, gy, gz)
  File "/usr/lib/python2.7/dist-packages/scipy/interpolate/fitpack.py", line 782, in bisplrep
    tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)
TypeError: integer argument expected, got float

J'ai essayé d'exécuter le code proposé sur cette page : doc scipy bisplrep et ça fonctionne très bien ...

Une idée ?

J'ai mis le fichier resultats.dat en pièce jointe sauf que j'ai du l'appeler resultats.txt pour qu'il soit accepter.

Merci

28/06/2011, 08h45
VinsS
Salut,

As-tu essayé
Code:

1 2 tck = bisplrep( int(gx), int(gy), int(gz))
comme semble le suggérer le message d'erreur ?

Bonjour

Non, mais le problème c'est que gx, gy et gz ne sont pas des entiers ! Il s'agit de flottants donc même si ça marche ça n'aura aucun sens par rapport à mon application.

Cependant j'ai fais le test et j'ai la même erreur. Visiblement bisplrep est une surcouche d'une fonction dans fitpack qui pose problème.

Ce qui me perturbe, c'est que le code suivant fonctionne et je ne vois pas la différence :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 import numpy as np from scipy import interpolate import matplotlib.pyplot as plt # Define function over sparse 20x20 grid x,y = np.mgrid[-1:1:20j,-1:1:20j] z = (x+y)*np.exp(-6.0*(x*x+y*y)) # Interpolate function over new 70x70 grid xnew,ynew = np.mgrid[-1:1:70j,-1:1:70j] tck = interpolate.bisplrep(x,y,z,s=0) znew = interpolate.bisplev(xnew[:,0],ynew[0,:],tck) plt.figure() plt.pcolor(xnew,ynew,znew) plt.colorbar() plt.title("Interpolated function.") plt.show()
Merci

28/06/2011, 13h21
Alexis.M

Je pense que ça viens du fait que bisplrep attend des tableaux a une dimension.

D'après la doc :

Code:

1
2
3
4
5
6
7
8
9
10
scipy.interpolate.bisplrep(x, y, z, w=None, xb=None, xe=None, yb=None, ye=None, kx=3, ky=3, task=0, s=None, eps=1e-16, tx=None, ty=None, full_output=0, nxest=None, nyest=None, quiet=1)
 
    Find a bivariate B-spline representation of a surface.
 
    Given a set of data points (x[i], y[i], z[i]) representing a surface z=f(x,y), compute a B-spline representation of the surface. Based on the routine SURFIT from FITPACK.
    Parameters :	
 
    x, y, z : ndarray
 
        Rank-1 arrays of data points.

Pourtant si je reprend le code qui fonctionne:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
 
# Define function over sparse 20x20 grid
x,y = np.mgrid[-1:1:20j,-1:1:20j]
z = (x+y)*np.exp(-6.0*(x*x+y*y))
 
print("x.shape = {0}".format(x.shape))
print("y.shape = {0}".format(y.shape))
print("z.shape = {0}".format(z.shape))
 
# Interpolate function over new 70x70 grid
xnew,ynew = np.mgrid[-1:1:70j,-1:1:70j]
tck = interpolate.bisplrep(x,y,z,s=0)
znew = interpolate.bisplev(xnew[:,0],ynew[0,:],tck)
 
plt.figure()
plt.pcolor(xnew,ynew,znew)
plt.colorbar()
plt.title("Interpolated function.")
plt.show()

En sortie j'ai :

Code:

1
2
3
x.shape = (20, 20)
y.shape = (20, 20)
z.shape = (20, 20)

Donc les entrées de bisplrep ne sont pas des tableaux 1D ... J'ai quand même fais le test et j'ai toujours la même erreur.

28/06/2011, 14h08
Alexis.M

Citation:

Code:

zspline = bisplev( xn[:,0], y[0,:], tck)

Ne serait-ce pas plutôt :

Code:

zspline = bisplev( xn[:,0], yn[0,:], tck)
28/06/2011, 14h28
Ladgalen

Si très juste :roll:

Mais l'erreur est avant, au niveau de bisplrep et pas bisplev donc l'erreur est toujours là :cry:

J'ai corrigé le code sur mon premier post.

En fait pour moi ça fonctionne avec un

"/usr/lib/python2.6/dist-packages/scipy/interpolate/fitpack.py:763: DeprecationWarning: integer argument expected, got float"

j'utilise la version 0.7.0 de scipy

Au passage, si tu retires "s=0" dans le code qui marche, tu a la même erreur.

As-tu essayer d'utiliser SmoothBivariateSpline ?

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#!/usr/bin/python
# -*- coding=utf-8 -*-
 
import scipy as sp
 
nx = 29
ny = 34
 
def spline(data):
 
    from scipy.interpolate import SmoothBivariateSpline
 
    # data
    x, y, z = data.transpose()
 
 
    m = nx
    print("s = [{0}:{1}]".format(m-sp.sqrt(2.*m), m + sp.sqrt(2.*m)))
 
    # create mgrid type object
    gx = x.reshape( nx, ny)
    gy = y.reshape( nx, ny)
    gz = z.reshape( nx, ny)
 
 
    # create spline object
    spl = SmoothBivariateSpline(x,y,z)
 
    # new grid
    xn, yn = sp.mgrid[x.min():x.max():50j,y.min():y.max():50j]
    zspline = spl( xn[:,0], yn[0,:])
 
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot_wireframe( gx, gy, gz)
    ax.plot_wireframe( xn, yn, zspline)
    plt.show()
 
if __name__ == "__main__":
    # load data
    data = sp.loadtxt("resultats.dat", usecols = (0,1,3))
    spline(data)

28/06/2011, 17h51
Ladgalen

Nickel ça marche ! :ccool:

C'est un peu comme la version POO de l'interpolation avec UnivariateSpline qu'ils présentent dans la doc. J'avais essayé avec la classe BivariateSpline dont ils parlent mais ça ne marchait pas. Visiblement le bon nom pour cette classe est celle dont tu parles :

http://docs.scipy.org/doc/scipy/refe...ivariatespline

J'ai aussi trouvé une classe RectBivariateSpline:

http://docs.scipy.org/doc/scipy/refe...ivariatespline

Cependant je ne vois pas encore la différence entre la Smoth et la Rect (en dehors de la syntaxe ) ... Visiblement le paramètre s est définit dans les deux cas mais par défaut il vaut 0 dans un cas et None dans l'autre. J'arrive à lui donner une valeur sans que ça plante ...

Merci en tout cas
28/06/2011, 21h18
Alexis.M

Manifestement Rect est fait pour lorsque tes données sont définies sur une grille alors que Smooth les positions peuvent être quelconques.
29/06/2011, 00h22
Ladgalen

Merci pour tout