Changeset 2754


Ignore:
Timestamp:
01/04/10 15:37:05 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added colors in rift plots (also support options: density,colorlevels,...)

Location:
issm/trunk/src/m/classes/public/plot
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/classes/public/plot/plot_quiver.m

    r2649 r2754  
    1010%      plot_quiver(md.x,md.y,md.vx,md.vy,options);
    1111
    12 %keep only non NaN elements
    13 pos=find(~isnan(x) & ~isnan(y) & ~isnan(u) & ~isnan(v));
    14 x=x(pos); y=y(pos);
    15 u=u(pos); v=v(pos);
    16 
    17 %get norm Min and Max
    18 Norm=sqrt(u.^2+v.^2);
    19 Min=min(Norm);
    20 Max=max(Norm);
    21 
    22 %process options: scaling factor?
    23 scalingfactor=getfieldvalue(options,'scaling',0.40);
    24 
    25 %number of colors?
    26 colorlevels=getfieldvalue(options,'colorlevels',30);
    27 if isnumeric(colorlevels),
    28         if isnan(colorlevels),
    29                 numcolors=30;
    30         else
    31                 numcolors=colorlevels;
    32         end
    33         levels=round_ice(linspace(Min,Max,numcolors+1),2);
    34 else
    35         levels=zeros(1,length(colorlevels)+2);
    36         levels(1)=Min;
    37         for i=1:length(colorlevels)
    38                 levels(i+1)=colorlevels{i};
    39         end
    40         levels(end)=Max;
    41         levels=sort(unique(levels));
    42         numcolors=length(levels)-1;
    43 end
    44 
    45 %set the colormap
    46 if numcolors==2;
    47         %blue and red
    48         c=colormap([0 0 1;1 0 0]);
    49 elseif numcolors==3,
    50         %blue yellow and red
    51         c=colormap([0 0 1;1 1 0;1 0 0]);
    52 else
    53         %let jet choose
    54         c=colormap(jet(numcolors));
    55 end
    56 
    57 %Scale data
    58 if strcmpi(getfieldvalue(options,'autoscale','on'),'off'),
    59         delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x);
    60         u=scalingfactor*sqrt(delta)*u./Norm;
    61         v=scalingfactor*sqrt(delta)*v./Norm;
    62 else
    63         delta=((min(x)-max(x))^2+(min(y)-max(y))^2)/numel(x);
    64         u=scalingfactor*sqrt(delta)*u./max(Norm);
    65         v=scalingfactor*sqrt(delta)*v./max(Norm);
    66 end
     12%process fields
     13[quivers,palette]=quiver_process(x,y,u,v,options);
    6714
    6815%loop over the number of colors
    6916hold on
    7017h=[];
    71 for i=1:numcolors
    72         pos=find( (Norm>=levels(i)) & (Norm<=levels(i+1)) );
    73         hprime=quiver(x(pos),y(pos),u(pos),v(pos),'Color',c(i,:),'ShowArrowHead','on','AutoScale','off');
     18for i=1:quivers.numcolors
     19        pos=find(quivers.colorind==i);
     20        hprime=quiver(quivers.x(pos),quivers.y(pos),quivers.u(pos),quivers.v(pos),...
     21                'Color',palette(i,:),'ShowArrowHead','on','AutoScale','off');
    7422        h=[h;hprime];
    7523end
    7624
    7725%take care of colorbar
    78 if  strcmpi(getfieldvalue(options,'colorbar','on'),'on'),
    79 
    80         %build ticks
    81         hcb=colorbar('peer',gca,'location','EastOutside');
    82         ticklabel=cell(1,length(levels));
    83         for i=1:length(levels),
    84                 ticklabel{i}=num2str(round_ice(levels(i),3));
    85         end
    86         tickpos=1:numcolors+1;
    87 
    88         %remove ticks if to many have been created
    89         proportion=round(length(levels)/10);
    90         if proportion>1,
    91                 ticklabel=ticklabel(1:proportion:end);
    92                 tickpos=tickpos(1:proportion:end);
    93         end
    94 
    95         %draw colorbar
    96         set(hcb,'YTickLabel',ticklabel,'YTick',tickpos);
    97         %position
    98         if exist(options,'colorbarpos'),
    99                 set(hcb,'Position',getfieldvalue(options,'colorbarpos'));
    100         end
    101         %fontsize
    102         fontsize=getfieldvalue(options,'fontsize',14);
    103         set(hcb,'FontSize',fontsize);
    104 end
     26quiver_colorbar(quivers,options);

  • issm/trunk/src/m/classes/public/plot/plot_riftrelvel.m

    r2745 r2754  
    77%   See also: PLOTMODEL
    88
     9%some checks
     10if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids),
     11        error('plot_riftvel error message: vx and vy do not have the right size'),
     12end
     13if ~isstruct(md.rifts),
     14        error('plot error message: no rifts available!');
     15end
     16options=addfielddefault(options,'scaling',2);
     17
     18%set as NaN all velocities not on rifts
     19for i=1:size(md.rifts,1),
     20        penaltypairs=md.rifts(i).penaltypairs(:,[1 2]);
     21        u=NaN*ones(md.numberofgrids,1);
     22        v=NaN*ones(md.numberofgrids,1);
     23        u(md.rifts(i).penaltypairs(:,1))=md.vx(penaltypairs(:,1))-md.vx(penaltypairs(:,2));
     24        v(md.rifts(i).penaltypairs(:,1))=md.vy(penaltypairs(:,1))-md.vy(penaltypairs(:,2));
     25end
     26
    927%process data and model
    1028[x y z elements is2d]=processmesh(md,options);
     29[vel isongrid isquiver]=processdata(md,[u v],options);
     30[quivers,palette]=quiver_process(x,y,vel(:,1),vel(:,2),options);
    1131
     32%prepare plot
    1233subplot(width,width,i);
    1334hold on
    … …  
    2243isp2=0;
    2344
    24 if isstruct(md.rifts),
    25         %plot mesh boundaries
    26         for i=1:size(md.segments,1),
    27                 h1=plot(x(md.segments(i,1:2)),y(md.segments(i,1:2)),'b-');
     45%plot mesh boundaries
     46for i=1:size(md.segments,1),
     47        h1=plot(x(md.segments(i,1:2)),y(md.segments(i,1:2)),'b-');
     48end
     49
     50for i=1:size(md.rifts,1),
     51        %get grids on rift
     52        penaltypairs=md.rifts(i).penaltypairs;
     53
     54        segments=md.rifts(i).segments;
     55        for j=1:size(segments,1),
     56                plot(x(segments(j,1:2)),y(segments(j,1:2)),'b-');
    2857        end
    29         for i=1:size(md.rifts,1),
    30                 penaltypairs=md.rifts(i).penaltypairs;
    3158
    32                 segments=md.rifts(i).segments;
    33                 for j=1:size(segments,1),
    34                         plot(x(segments(j,1:2)),y(segments(j,1:2)),'b-');
     59        normal=zeros(2,1);
     60        for j=1:size(penaltypairs,1),
     61                normal(1)=penaltypairs(j,5);
     62                normal(2)=penaltypairs(j,6);
     63
     64                vx1=md.vx(penaltypairs(j,1)); vx2=md.vx(penaltypairs(j,2)); vy1=md.vy(penaltypairs(j,1)); vy2=md.vy(penaltypairs(j,2));
     65                penetration=(vx2-vx1)*normal(1)+(vy2-vy1)*normal(2);
     66                %if penetration is negative, plot in black, positive, plot in red;: ie: if rift is closing, black, if rift is opening, red.
     67                if(penetration>0),
     68                        p2=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'r*');
     69                        isp2=1;
     70                else
     71                        p1=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'k*');
     72                        isp1=1;
    3573                end
     74        end
    3675
    37                 normal=zeros(2,1);
    38                 for j=1:size(penaltypairs,1),
    39                         normal(1)=penaltypairs(j,5);
    40                         normal(2)=penaltypairs(j,6);
     76        %point out the tips
     77        h2=plot(x(md.rifts(i).tips(1)),y(md.rifts(i).tips(1)),'g*');
     78        plot(x(md.rifts(i).tips(2)),y(md.rifts(i).tips(2)),'g*');
     79        segments=md.rifts(i).segments(:,1:2);
     80end
    4181
    42                         vx1=md.vx(penaltypairs(j,1)); vx2=md.vx(penaltypairs(j,2)); vy1=md.vy(penaltypairs(j,1)); vy2=md.vy(penaltypairs(j,2));
    43                         penetration=(vx2-vx1)*normal(1)+(vy2-vy1)*normal(2);
    44                         %if penetration is negative, plot in black, positive, plot in red;: ie: if rift is closing, black, if rift is opening, red.
    45                         if(penetration>0),
    46                                 p2=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'r*');
    47                                 isp2=1;
    48                         else
    49                                 p1=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'k*');
    50                                 isp1=1;
    51                         end
    52                 end
     82%plot rifts vel
     83h3=[];
     84for i=1:quivers.numcolors
     85        pos=find(quivers.colorind==i);
     86        hprime=quiver(quivers.x(pos),quivers.y(pos),quivers.u(pos),quivers.v(pos),...
     87                'Color',palette(i,:),'ShowArrowHead','on','AutoScale','off');
     88        hprime=quiver(quivers.x(pos),quivers.y(pos),-quivers.u(pos),-quivers.v(pos),...
     89                'Color',palette(i,:),'ShowArrowHead','on','AutoScale','off');
     90        h3=[h3;hprime];
     91end
    5392
    54                 %point out the tips
    55                 h2=plot(x(md.rifts(i).tips(1)),y(md.rifts(i).tips(1)),'g*');
    56                 plot(x(md.rifts(i).tips(2)),y(md.rifts(i).tips(2)),'g*');
    57 
    58                 %plot relative velocities
    59                 h=quiver(x(penaltypairs(:,1)),y(penaltypairs(:,1)),md.vx(penaltypairs(:,2))-md.vx(penaltypairs(:,1)),md.vy(penaltypairs(:,2))-md.vy(penaltypairs(:,1)),3);
    60                 %frequency=1:5:size(penaltypairs,1);
    61                 %h=quiver(x(penaltypairs(frequency,1)),y(penaltypairs(frequency,1)),md.vx(penaltypairs(frequency,2))-md.vx(penaltypairs(frequency,1)),md.vy(penaltypairs(frequency,2))-md.vy(penaltypairs(frequency,1)),2);
    62                 set(h,'Color',[0 1 0]);
    63         end
    64         if isp1 & isp2
    65                 legend([h1,h2,p1,p2,h],'mesh boundaries','rift tips',' rifts closing','rifts opening','relative velocities')
    66         elseif isp1
    67                 legend([h1,h2,p1,h],'mesh boundaries','rift tips',' rifts closing','relative velocities')
    68         elseif isp2
    69                 legend([h1,h2,p2,h],'mesh boundaries','rift tips','rifts opening','relative velocities')
    70         else
    71                 legend([h1,h2,h],'mesh boundaries','rift tips','relative velocities')
    72         end
     93%legend
     94if isp1 & isp2
     95        legend([h1,h2,p1,p2],'mesh boundaries','rift tips',' rifts closing','rifts opening')
     96elseif isp1
     97        legend([h1,h2,p1],'mesh boundaries','rift tips',' rifts closing')
     98elseif isp2
     99        legend([h1,h2,p2],'mesh boundaries','rift tips','rifts opening')
    73100else
    74         error('plot error message: no rifts available!');
     101        legend([h1,h2],'mesh boundaries','rift tips')
    75102end
    76103hold off
    77104
    78105%apply options
    79 options=addfielddefault(options,'title','Rift relative velocities');
    80 options=addfielddefault(options,'colorbar',0);
     106quiver_colorbar(quivers,options);
     107options=changefieldvalue(options,'colorbar',2);
     108options=addfielddefault(options,'title','Rift Velocities');
    81109applyoptions(md,[],options);

  • issm/trunk/src/m/classes/public/plot/plot_riftvel.m

    r2695 r2754  
    77%   See also: PLOTMODEL
    88
     9%some checks
     10if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids),
     11        error('plot_riftvel error message: vx and vy do not have the right size'),
     12end
     13if ~isstruct(md.rifts),
     14        error('plot error message: no rifts available!');
     15end
     16options=addfielddefault(options,'scaling',2);
     17
     18%set as NaN all velocities not on rifts
     19for i=1:size(md.rifts,1),
     20        penaltypairs=md.rifts(i).penaltypairs(:,[1 2]);
     21        u=NaN*ones(md.numberofgrids,1);
     22        v=NaN*ones(md.numberofgrids,1);
     23        u(penaltypairs(:))=md.vx(penaltypairs(:));
     24        v(penaltypairs(:))=md.vy(penaltypairs(:));
     25end
     26
    927%process data and model
    1028[x y z elements is2d]=processmesh(md,options);
     29[vel isongrid isquiver]=processdata(md,[u v],options);
     30[quivers,palette]=quiver_process(x,y,vel(:,1),vel(:,2),options);
    1131
     32%prepare plot
    1233subplot(width,width,i);
    1334hold on
    … …  
    2243isp2=0;
    2344
    24 if isstruct(md.rifts),
    25         %plot mesh boundaries
    26         for i=1:size(md.segments,1),
    27                 h1=plot(x(md.segments(i,1:2)),y(md.segments(i,1:2)),'b-');
     45%plot mesh boundaries
     46for i=1:size(md.segments,1),
     47        h1=plot(x(md.segments(i,1:2)),y(md.segments(i,1:2)),'b-');
     48end
     49
     50for i=1:size(md.rifts,1),
     51        %get grids on rift
     52        penaltypairs=md.rifts(i).penaltypairs;
     53
     54        segments=md.rifts(i).segments;
     55        for j=1:size(segments,1),
     56                plot(x(segments(j,1:2)),y(segments(j,1:2)),'b-');
    2857        end
    29         for i=1:size(md.rifts,1),
    30                 penaltypairs=md.rifts(i).penaltypairs;
    3158
    32                 segments=md.rifts(i).segments;
    33                 for j=1:size(segments,1),
    34                         plot(x(segments(j,1:2)),y(segments(j,1:2)),'b-');
     59        normal=zeros(2,1);
     60        for j=1:size(penaltypairs,1),
     61                normal(1)=penaltypairs(j,5);
     62                normal(2)=penaltypairs(j,6);
     63
     64                vx1=md.vx(penaltypairs(j,1)); vx2=md.vx(penaltypairs(j,2)); vy1=md.vy(penaltypairs(j,1)); vy2=md.vy(penaltypairs(j,2));
     65                penetration=(vx2-vx1)*normal(1)+(vy2-vy1)*normal(2);
     66                %if penetration is negative, plot in black, positive, plot in red;: ie: if rift is closing, black, if rift is opening, red.
     67                if(penetration>0),
     68                        p2=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'r*');
     69                        isp2=1;
     70                else
     71                        p1=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'k*');
     72                        isp1=1;
    3573                end
     74        end
    3675
    37                 normal=zeros(2,1);
    38                 for j=1:size(penaltypairs,1),
    39                         normal(1)=penaltypairs(j,5);
    40                         normal(2)=penaltypairs(j,6);
     76        %point out the tips
     77        h2=plot(x(md.rifts(i).tips(1)),y(md.rifts(i).tips(1)),'g*');
     78        plot(x(md.rifts(i).tips(2)),y(md.rifts(i).tips(2)),'g*');
     79        segments=md.rifts(i).segments(:,1:2);
     80end
    4181
    42                         vx1=md.vx(penaltypairs(j,1)); vx2=md.vx(penaltypairs(j,2)); vy1=md.vy(penaltypairs(j,1)); vy2=md.vy(penaltypairs(j,2));
    43                         penetration=(vx2-vx1)*normal(1)+(vy2-vy1)*normal(2);
    44                         %if penetration is negative, plot in black, positive, plot in red;: ie: if rift is closing, black, if rift is opening, red.
    45                         if(penetration>0),
    46                                 p2=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'r*');
    47                                 isp2=1;
    48                         else
    49                                 p1=plot(x(penaltypairs(j,1)) ,y(penaltypairs(j,1)),'k*');
    50                                 isp1=1;
    51                         end
    52                 end
     82%plot rifts vel
     83h3=[];
     84for i=1:quivers.numcolors
     85        pos=find(quivers.colorind==i);
     86        hprime=quiver(quivers.x(pos),quivers.y(pos),quivers.u(pos),quivers.v(pos),...
     87                'Color',palette(i,:),'ShowArrowHead','on','AutoScale','off');
     88        h3=[h3;hprime];
     89end
    5390
    54                 %point out the tips
    55                 h2=plot(x(md.rifts(i).tips(1)),y(md.rifts(i).tips(1)),'g*');
    56                 plot(x(md.rifts(i).tips(2)),y(md.rifts(i).tips(2)),'g*');
    57                 segments=md.rifts(i).segments(:,1:2);
    58                 h3=quiver(x(segments),y(segments),md.vx(segments),md.vy(segments));
    59                 set(h3,'Color',[1 0 0]);
    60         end
    61         if isp1 & isp2
    62                 legend([h1,h2,p1,p2,h3],'mesh boundaries','rift tips',' rifts closing','rifts opening','velocities')
    63         elseif isp1
    64                 legend([h1,h2,p1,h3],'mesh boundaries','rift tips',' rifts closing','velocities')
    65         elseif isp2
    66                 legend([h1,h2,p2,h3],'mesh boundaries','rift tips','rifts opening','velocities')
    67         else
    68                 legend([h1,h2,h3],'mesh boundaries','rift tips','velocities')
    69         end
     91%legend
     92if isp1 & isp2
     93        legend([h1,h2,p1,p2],'mesh boundaries','rift tips',' rifts closing','rifts opening')
     94elseif isp1
     95        legend([h1,h2,p1],'mesh boundaries','rift tips',' rifts closing')
     96elseif isp2
     97        legend([h1,h2,p2],'mesh boundaries','rift tips','rifts opening')
    7098else
    71         error('plot error message: no rifts available!');
     99        legend([h1,h2],'mesh boundaries','rift tips')
    72100end
    73101hold off
    74102
    75103%apply options
     104quiver_colorbar(quivers,options);
     105options=changefieldvalue(options,'colorbar',2);
    76106options=addfielddefault(options,'title','Rift Velocities');
    77 options=addfielddefault(options,'colorbar',0);
    78107applyoptions(md,[],options);
Note: See TracChangeset for help on using the changeset viewer.