perimeter_new/plot_ros.m

   1 function plot_ros(long,lat,ros)
   2 % plot_ros(long,lat,ros)
   3 % in:
   4 %   long    east-west coordinates (m), size (m,n)
   5 %   lat     south-north coordinates (m) size (m,n)
   6 %   ros     rate of spread in 8 directions,from read_ros_from_wrfout
   7 %           size (m,n,3,3)
   8
   9 % resample
  10 step=20;
  11
  12 long_r=long(1:step:end,1:step:end);
  13 lat_r=lat(1:step:end,1:step:end);
  14 ros_r=ros(1:step:end,1:step:end,:,:);
  15 m_r=size(ros_r,1);
  16 n_r=size(ros_r,2);
  17
  18 d=zeros(8,2);
  19 u_r=zeros(m_r,n_r);
  20 v_r=zeros(m_r,n_r);
  21 for i=1:m_r
  22     for j=1:n_r
  23         k=0;
  24         for a=1:3
  25             for b=1:3
  26                 u=a-2;v=b-2;s=norm([u,v]);
  27                 if s>0,
  28                     k=k+1;
  29                     u=u/s;v=v/s;
  30                     d(k,:)=[ros_r(i,j,a,b)*u,ros_r(i,j,a,b)*v];
  31                 end
  32             end
  33         end
  34     % now d has the ends of ros vectors as rows
  35     center=mean(d);
  36     C=cov(d);
  37     [V,D]=eig(C);
  38     % components of the vector of max spread
  39     [lambda1,i1]=max(diag(D));   % max eigenvalue
  40     s1=sqrt(lambda1);            % singular value
  41     ax=V(:,i1);                 % the eigenvector
  42     vertex=(center+s1*ax'*sign(center*ax));  % vertex of the ellipse away from the center
  43     u_r(i,j)=vertex(1);
  44     v_r(i,j)=vertex(2);
  45     end
  46 end
  47
  48 % scale
  49 dx=long_r(2:end,:)-long_r(1:end-1,:);dx=mean(dx(:));
  50 dy=lat_r (:,2:end)-lat_r (:,1:end-1);dy=mean(dy(:));
  51 maxros=max(sqrt((u_r(:).^2+v_r(:).^2)));
  52 scaleros=min(dx,dy)/6;
  53
  54 % plot
  55 clf
  56 quiver(long_r,lat_r,scaleros*u_r,scaleros*v_r,'k')
  57 %hold on
  58 %pcolor(long,lat,mr); shading flat
  59 %hold off
  60 % drawnow
  61 %colorbar
  62 end