Interpolation des épaisseurs de retombées volcaniques

**Jean Dumoncel** · 19/12/2011, 20h00

Et y-a-t-il un rapport avec l'interpolation avec le polynôme? As-tu essayé une interpolation 3D?

**LudOL56** · 20/12/2011, 15h22

Et y-a-t-il un rapport avec l'interpolation avec le polynôme? As-tu essayé une interpolation 3D?

Un rapport ? je ne comprend pas ta question magelan...
En fait j'ai essayé d'utiliser la fonction interp2 mais celle ci ne marche pas :/. Mais j'ai pas essayé d'interpolation 3D...

**Jean Dumoncel** · 20/12/2011, 17h24

Par rapport à tes données de départ, pour chaque point tu connais son abscisse, son ordonnée, sa hauteur (épaisseur) et sa distance au cratère, c'est bien cela?

Comment est défini la distance au cratère? c'est par rapport au centre du cratère?

J'ai fait une erreur, j'aurai dû parler d'interpolation 2D. Qu'est-ce qui ne marche pas avec interp2? J'ai vu que tu avais utilisé griddata aussi, cela ne correspondait pas à tes besoins?

**LudOL56** · 20/12/2011, 18h46

Déjà, merci beaucoup magelan de t’intéresser à mon problème. C'est très important pour moi de réussir ce que je cherche à faire.

Effectivement, chaque point de coordonnée X,Y possède une épaisseur mesurée. La distance au cratère a été calculée, comme tu l'as dit, par rapport au conduit éruptif du volcan à l'origine de l’éruption.

En faite la fonction griddata est pas mal mais j'aimerais que mes courbes d'isovaleurs soit lissées. Et dans la documentation de la fonction "contour" présente dans le logiciel, il est marqué qu'en utilisant la fonction interp2 selon le code suivant on pouvait lisser les isovaleurs:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
%Smoothing Contour Data

%Use interp2 to smooth contour lines. Also set the contour label text BackgroundColor to a light yellow and the EdgeColor to light gray.

Z = peaks;
[C,h] = contour(interp2(Z,4));
text_handle = clabel(C,h);
set(text_handle,'BackgroundColor',[1 1 .6],...
    'Edgecolor',[.7 .7 .7])

Voici ici ce que j'ai fait pour utiliser la fonction interp2:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
%figure
%[X,Y]=meshgrid(x,y);
%[XI,YI]=meshgrid(x,y);
%ZI=interp2(X,Y,Z,XI,YI,'linear');
%mesh(x,y,z), hold, mesh(XI,YI,ZI+15)
%hold off

J'obtiens l'erreur suivante:

??? Error using ==> interp2 at 145
When Z is a vector, X and Y must also be vectors.

Error in ==> Prog_sirs at 55
ZI=interp2(X,Y,Z,XI,YI,'linear');

Merci

**Jean Dumoncel** · 20/12/2011, 19h54

Si tes coordonnées sont répartis de façon non régulière, il faut utiliser griddata pour réaliser l'interpolation. Il faut d'abord définir Xi et Yi à l'aide de meshgrid pour définir tes points d'interpolation puis utiliser griddata avec tes vecteurs x y et Z :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

ZI = griddata(x,y,Z,XI,YI)

mais j'ai l'impression que tu l'avais testé au post #5, non?

**LudOL56** · 20/12/2011, 20h08

Envoyé par magelan

Si tes coordonnées sont répartis de façon non régulière, il faut utiliser griddata pour réaliser l'interpolation. Il faut d'abord définir Xi et Yi à l'aide de meshgrid pour définir tes points d'interpolation puis utiliser griddata avec tes vecteurs x y et Z :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

ZI = griddata(x,y,Z,XI,YI)

mais j'ai l'impression que tu l'avais testé au post #5, non?

Ah ok, donc je peux pas utiliser interp2 alors...non ?
Effectivement j'ai utilisé griddata et ça donne ce que j'ai joint

!
Mais ça lisse pas du tout et c'est pas vraiment ce que je veux...

**Jean Dumoncel** · 20/12/2011, 20h32

Si tes points ne sont pas répartis de façon régulière (ie espace constant entre les points) alors non, tu ne peux pas utiliser interp2. Mais griddata fait la même chose pour les grille non régulière.

Peux-tu nous montrer le code associé?

**LudOL56** · 20/12/2011, 20h47

Voici le code associé :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
[XI,YI]=meshgrid(x,y);
ZI=griddata(x,y,z,XI,YI);
figure
%contour3(X,Y,ZI,100)
 
[C,h] = contour(X,Y,ZI);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*25)
colormap cool
hold on
plot(x,y)

Du coup as tu une idée pour lisser les contours des isovaleurs ?
Normalement les isopaques dans ce cas devraient avoir une forme elliptiques...

**Jean Dumoncel** · 20/12/2011, 21h15

Si c'est possible, peux-tu joindre un fichier contenant les données que tu traites dans ton prochain message pour que l'on puisse faire des tests?

As-tu essayé de creer des grilles x et y plus fine :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

[Xi Yi] = meshgrid(xmin:pasx:xmax,ymin:pasy:ymax)

ou pasx et pasy sont les écarts entre 2 valeurs?

**LudOL56** · 20/12/2011, 21h20

Non je n'ai pas essayé d'affiner mes grilles...je vais tenter ça !
Je te mets mon fichier au format .txt en pièce jointe

.

**Jean Dumoncel** · 20/12/2011, 21h45

En affinant un peu la grille avec tes données :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
load Ech_sirs.txt
 
plot3(Ech_sirs(:,1),Ech_sirs(:,2),Ech_sirs(:,3),'*')
x = Ech_sirs(:,1);
y = Ech_sirs(:,2);
z = Ech_sirs(:,3);
 
 
[XI,YI]=meshgrid(min(x):100:max(x),min(y):100:max(y));
ZI=griddata(x,y,z,XI,YI,'cubic');
 
figure,surf(XI,YI,ZI,'EdgeColor','none')
hold on, plot3(x,y,z,'*')
 
figure,contour(XI,YI,ZI,30)
 
figure,contour3(XI,YI,ZI,30)

Ca correspond à ce que tu souhaites?

**LudOL56** · 20/12/2011, 21h56

Ouai ça commence à ressembler à ce que je cherche ! C'est même ça d'ailleurs, mais peut être que je reviendrais vers toi pour d'autres conseils

!

Encore merci magelan !

**LudOL56** · 21/12/2011, 13h11

Bonjour,

C'est encore moi, je voudrais créer une boite de dialogue qui permettrai de demander à l'utilisateur quelle méthode d'interpolation utiliser: cubic, linear, v4 ou nearest ?

Pour l'instant j'ai fait ça mais je ne sais plus comment faire pour integrer ça dans les lignes de codes qui permettent de tracer les cartes...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
prompt={'cubic, nearest, v4 ou linear ?'};
dlg_title='Méthode d''interpolation';
num_lines=1;
answer = inputdlg(prompt,dlg_title,num_lines);

En faite ce que je voudrais c'est que, quand l'utilisateur tape "cubic" dans la boite de dialogue, le programme n'execute uniquement les carte interpolées réalisées avec la méthode cubic...

Merci

**Jean Dumoncel** · 21/12/2011, 16h09

Tu peux utiliser la structure switch...case :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
prompt={'cubic, nearest, v4 ou linear ?'};
dlg_title='Méthode d''interpolation';
num_lines=1;
answer = inputdlg(prompt,dlg_title,num_lines);
switch answer{1}
    case 'cubic'
        ...
    case 'nearest'
        ...
    case 'v4'
        ...
    case 'linear'
        ...
    otherwise
        ...
end

**LudOL56** · 21/12/2011, 17h39

Ça marche très bien, merci

!
Sais-tu si c'est possible de faire une coupe à travers le graphe avec la surface interpolée afin de montrer la décroissance des épaisseurs avec la distance ?

Parce que j'aurais des calculs à faire dessus pour à terme remonter au débit de l'éruption...

Merci

**Jean Dumoncel** · 21/12/2011, 19h50

Une coupe dans quelle direction?

tes lignes et tes colonnes de ZI sont autant de coupes sur les axes x et y déjà.

**LudOL56** · 21/12/2011, 22h42

Okay donc juste avec un plot je devrais y arriver...merci !

**LudOL56** · 11/01/2012, 17h36

Bonjour,

C'est encore moi, je n'arrive pas à tracer une coupe comme je le veux (pièce jointe).
En faite j'aimerais tracer une coupe avec en x la longitude, en y l'épaisseur et le tout à une latitude donnée. Comment puis-je faire ?

Merci !

**Jean Dumoncel** · 13/01/2012, 13h50

Comme dit précédemment, il te suffit de selectionner les données correspondantes dans tes tableaux, par exemple :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
figure,contour3(XI,YI,ZI,30)
hold on,
coupe = 300;
plot3(XI(coupe,:),YI(coupe,:),ZI(coupe,:))

et pour la tracer en 2D :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

plot(XI(coupe,:),ZI(coupe,:))

Interpolation des épaisseurs de retombées volcaniques

MATLAB

Discussions similaires

Partager

Partager