Hello,
I would like to do a grid search by assining zeros to each grid point first. And then I would like to calculate required values for each grid point and finally I am required to plot the grid with relevant values. here my grid represent error values for each grid point. So I want to plot this grid using some colour scheme method.
Attached here is the type of graph I am looking for.
Could someone please help me with this?
Below is the code I wrote for this. However I am not sure how much of assigning zeros, for loop and then plotting is correct.
I tried to use plot and quiver fuctions but different length x and y axis appeared to be an error in that case. For now i used surf function, though i am not sure if that is correct. But it gives me a z axis too, which is not what i want.
Please let me know if someone could help me with this.
% Obserevd Z location km/Myr ## 176.7 W, 19.2 S, 64 mm/yr (km/Myr) velocity;
Zlon=-176.324;
Zlat=-18.395;
Zobsw=64;
R=6371; %km
% evaluate error at each grid point
AClon =-177.5:0.01:-176.7;
AClat =-20.5:0.01:-19;
[X,Y]=meshgrid(AClon,AClat);
AC=zeros(151,81);
for i=[1:151]
for j=[1:81]
arclen = distance(Zlon,Zlat,AClon(j),AClat(i)); % arc distant between 1.Z location and 2.each grid point.
w =(Zobsw/(R*sin(theta*pi/180)))*180/pi; % deg/Myr
vl = plate_vel(Zlon,Zlat,AClon(j),AClat(i),w); % i found this plate-vel funtion online. please use any set of data as an example here.
misfit=(sum((Zobsw-vl).^2)./(numel(vl)));
AC(i,j)=misfit;
end
end
surf(X,Y,AC)

18 commentaires

surf(X,Y,AC) uses AC for the color data and surface height. If you want a flat surface (e.g., all ZData are 0), you can do this:
surf(X,Y,zeros(size(X)),AC);
Then set the axes view to be from directly above the surface (looking down on the X-Y plane):
set(gca(),'View',[0 90]);
I'm not sure if this helps you get what you're after though.
Anitha Limann
Anitha Limann le 18 Déc 2021
thank you for your kind response. using your code i get this kind of an output. But i want to plot error calculations on each x,y intersection point (which are location cordinates) instead of squares like this. would you be able to help me with that?
Voss
Voss le 18 Déc 2021
Modifié(e) : Voss le 18 Déc 2021
I'm not sure what you mean by "plot error calculations". Do you want a dot or a circle, say, at each (x,y) point, with the color of the dot or circle corresponding to the error value?
If so, you can try this (to get dots):
scatter3(X(:),Y(:),AC(:),[],AC(:),'Marker','.');
set(gca(),'View',[0 90]); % set to X-Y view
Anitha Limann
Anitha Limann le 18 Déc 2021
Exactly. But instead of dots or circles with a color corresponding to the error value i want to asign different colours to those dots/circles corresponding to an error range. For an example blue for error values less than 5; red for error values 5-10; green for error values higher than 10 likewise;
Could you please help me with this?
Thank you
Put this after your for loops:
AC_group = zeros(size(AC));
% AC_group(AC < 5) = 0; % it is not necessary to specify this one
AC_group(AC >= 5 & AC <= 10) = 1;
AC_group(AC > 10) = 2;
scatter3(X(:),Y(:),AC_group(:),[],AC_group(:),'Marker','.');
set(gca(),'View',[0 90]);
set(gcf(),'Colormap',[0 0 1; 1 0 0; 0 1 0]); % color is [blue; red; green] <-> AC_group is [0; 1; 2]
Thanks again for your kind help.
AC_group(AC >= min(AC) & AC <= median(AC)) = 1;
AC_group(AC > median(AC)) = 2;
I tried with median value to group data but it seems like groups are not enough to represent any error variaotions based on location. I found the type of graph I was looking for in a book but not sure how to plot something like that. If you could please take a look, attached here is the type of graph I need
AC_median = median(AC(:));
AC_group = zeros(size(AC));
AC_group(AC <= AC_median) = 1;
AC_group(AC > AC_median) = 2;
scatter3(X(:),Y(:),AC_group(:),[],AC_group(:),'Marker','.');
set(gca(),'View',[0 90]);
set(gcf(),'Colormap',[0 0 1; 1 0 0]); % color is [blue; red] <-> AC_group is [1; 2]
Anitha Limann
Anitha Limann le 19 Déc 2021
Thank you. Please let me know if you know any method to obtain a graph that looks like in the attached image. That is how my error graph should looks like, but i am not sure how to plot like that.
That image looks more like a surface than a scatter plot. And with no thresholding into groups, and with a greyscale color map:
surf(X,Y,zeros(size(X)),AC,'EdgeColor','none');
set(gca(),'View',[0 90]);
colormap('gray');
See if that gives you something similar.
Anitha Limann
Anitha Limann le 19 Déc 2021
I guess something wrong with my code. When i use surf i get a map similar to very first map i attached in this thread. IT might be something wrong with my grid search section.
Anyway Thank you very much for your help. If you see any error somewhere please let me know. I am not sure how to show error changes with the location (x,y).
Voss
Voss le 19 Déc 2021
You can post the latest version of your code and I'll take a look at it when I can.
sure. Thank you. posted below is the code.
clear;clc;
% Obserevd Z location km/Myr ## 176.7 W, 19.2 S, 64 mm/yr (km/Myr) velocity;
Zlon=-176.324;
Zlat=-18.395;
Zomega=64;
R=6371; %km
% evaluate error at each grid point
AClon =-176.8:0.01:-176.5;
AClat =-19.5:0.01:-19;
[X,Y]=meshgrid(AClon,AClat);
AC=zeros(51,31);
for i=[1:51]
for j=[1:31]
theta = arcdist(Zlon,Zlat,AClon(j),AClat(i)); % arc distant between 1.Z location and 2.each grid point.
AComega =(Zomega/(R*sin(theta*pi/180)))*180/pi; % deg/Myr
[pvel,azi] = plate_motion(Zlon,Zlat,AClon(j),AClat(i),AComega);
misfit=(sum((Zomega-pvel).^2));
AC(i,j)=misfit;
end
end
a1 = median(AC(:));
a2 = median(min(AC(:)):a1);
a3 = median(a1:max(AC(:)));
AC_group = zeros(size(AC));
AC_group(AC <=a2)=1; % 4 groups. not sure if i need to make groups at all.
AC_group(AC >a2 & AC <=a1)=2;
AC_group(AC >a1 & AC <=a3)=3;
AC_group(AC >a3)=4;
% Both scatter and surf were used, since i am not sure how to the required type of a map.
figure(1)
scatter3(X(:),Y(:),AC_group(:),[],AC_group(:),'Marker','.');
set(gca(),'View',[0 90]);
set(gcf(),'Colormap',winter); % color is [blue; red] <-> AC_group is [1; 2]
figure(2)
surf(X,Y,AC);
% surf(X,Y,zeros(size(X)),AC,'EdgeColor','none');
% set(gca(),'View',[0 90]);
% colormap('winter');
Try replacing the last few lines, from figure(2) to the end, with this:
figure(2)
% surf(X,Y,AC);
surf(X,Y,zeros(size(X)),AC,'EdgeColor','none','FaceColor','interp');
set(gca(),'View',[0 90]);
set(gcf(),'Colormap',flip(gray(),1));
Does it look similar enough to the reference image? If not, what is different?
Anitha Limann
Anitha Limann le 19 Déc 2021
It does look quite similar.
I am trying to plot this graph to find the location (x,y) where I can add an euler pole (geophysical context). the Z location has an observed plate velocity value. So i conducted the grid search to back calculate plate velocity value at Z location correspondent to each grid point. Thus whereever (x,y) gives me the minimum errors would be my preffered locations.
This is why I first tried with grouping, however it does not seems to be working because, grouping criteria are not good enough I guess.
The correct locations should roughly fall along the upper part of the right diagonal line.
That is what makes me hesitate quite. It might be my error calculations are at fault. I will recheck them.
Thank you very much for your kind help with the graph.
Voss
Voss le 19 Déc 2021
Modifié(e) : Voss le 19 Déc 2021
No problem. Since the question regarding making an image like the reference image appears to be resolved, I'll make an answer with the relevant plotting code, which you can accept so that others may find it.
If you have additional questions about why the data/calculations might be incorrect, please post them as new questions.
Anitha Limann
Anitha Limann le 20 Déc 2021
Thank you so much for your hardworks
Voss
Voss le 20 Déc 2021
No problem! Don't hesitate to ask another question if you suspect your calculations are off.
Anitha Limann
Anitha Limann le 20 Déc 2021
Sure will do. Thank you!

Connectez-vous pour commenter.

 Réponse acceptée

Voss
Voss le 19 Déc 2021
surf(X,Y,zeros(size(X)),AC,'EdgeColor','none','FaceColor','interp');
set(gca(),'View',[0 90]);
set(gcf(),'Colormap',flip(gray(),1));

Plus de réponses (0)

Catégories

En savoir plus sur Graphics Performance dans Centre d'aide et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by