detection/patch4.m

   1 % to create conus.kml:
   2 % download http://firemapper.sc.egov.usda.gov/data_viirs/kml/conus_hist/conus_20120914.kmz
   3 % and gunzip
   4 %
   5 % to create w.mat:
   6 % run Adam's simulation, currently results in
   7 %
   8 % /share_home/akochans/WRF341F/wrf-fire/WRFV3/test/em_utfire_1d_med_4km_200m
   9 % then in Matlab
  10 % arrays needed only once
  11 % f='wrfout_d01_2013-08-20_00:00:00';
  12 % t=nc2struct(f,{'Times'},{});  n=size(t.times,2)
  13 % w=nc2struct(f,{'Times','TIGN_G','FXLONG','FXLAT','UNIT_FXLAT','UNIT_FXLONG','XLONG','XLAT','NFUEL_CAT'},{},n);
  14 % save ~/w.mat w
  15 %
  16 % array at fire resolution every saved timestep
  17 % to create s.mat:
  18 % a=dir('wrfout_d01*');
  19 % s=read_wrfout_sel({a.name},{'FGRNHFX',Times});
  20 % save ~/s.mat s
  21 %
  22 % arrays at atm resolution every saved timestep
  23 % to create ss.mat
  24 % a=dir('wrfout_d01*')
  25 % s=read_wrfout_sel({a.name},{'Times','UAH','VAH'})
  26 % save ss s
  27 %
  28 % fuels.m is created by WRF-SFIRE at the beginning of the run
  29
  30 % ****** REQUIRES Matlab 2013a - will not run in earlier versions *******
  31
  32 % run patch_load first
  33
  34 % figures
  35 figmap=1;
  36 printing=0;
  37 historical='all';
  38 %historical='previous';
  39
  40 % convert tign_g to datenum
  41 w.time=datenum(char(w.times)');
  42 red.tign=(red.tign_g - max(red.tign_g(:)))/(24*60*60) + w.time;
  43 min_tign=min(red.tign(:));
  44 max_tign=max(red.tign(:));
  45
  46 cmap2=cmap;
  47 cmap2(1:7,:)=NaN;
  48 [cmap,imax]=cmapmod14;
  49 granules{2}=[];
  50 figure(figmap);clf
  51 clear M
  52 for step=2:length(ss.time)  % over WRF frames
  53     figure(figmap);clf;hold off
  54     granules{1}=find(r.time <= ss.time(step));
  55     switch historical
  56         case 'all'
  57             % take all satellite overpasses from the beginning of time
  58             granules{2}=granules{1};
  59         case 'last'
  60             % last time step only
  61             granules{2}=find(r.time <= ss.time(step) & r.time > ss.time(step-1));
  62         case 'previous'
  63             % the previous granules persist until replaced
  64             new_det=find(r.time <= ss.time(step) & r.time > ss.time(step-1));
  65             if new_det,
  66                 granules{2}=new_det;
  67             end
  68         otherwise
  69             historical
  70             error('unknown parameter historical')
  71     end
  72     for ipass=1:2,
  73         det=granules{ipass};
  74         ndet=length(det);
  75         if printing>1,
  76             if det, fprintf('pass %i using satellite granules ',ipass),disp(granules{1}), end
  77             for idet=1:ndet
  78                 fprintf('%i %s\n',det(idet),datestr(r.time(det(idet))))
  79             end
  80         end
  81         % detections to now
  82         for idet=1:ndet,
  83             x=r.x{det(idet)}; % load fire detection image
  84             age=t-r.time(idet); % age of detection in days
  85             offset = min(imax,floor(4*age)); % offset in colormap for detection age
  86             dd=x.data(:)>6;  % indices of detected
  87             x.data(dd)=x.data(dd)+3*offset;   % transition to yellow
  88             if printing>1,
  89                 fprintf('step %i pass %i granule %i detections %i\n',step,ipass,idet,sum(dd))
  90             end
  91             if ipass==1, % build up the background
  92                showmod14(x)
  93             elseif any(dd), % all except fire transparent
  94                x.data(~dd)=0;
  95                showmod14(x)
  96             end
  97             hold on
  98         end % idet
  99     end % ipass
 100     u=ss.uh(:,:,step);
 101     v=ss.vh(:,:,step);
 102     maxw=max(sqrt(u(:).*u(:)+v(:).*v(:)));
 103     fprintf('step %i max windspeed %g granules %i\n',step,maxw,ndet)
 104     sc=0.006;quiver(w.xlong,w.xlat,sc*u,sc*v,0); % wind
 105     hold on
 106     t=ss.time(step);
 107     contour(red.fxlong,red.fxlat,red.tign,[t t],'-k'); % fireline
 108     title(datestr(t));
 109     axis(red.axis)
 110     avg_lat=0.5*(red.min_lat+red.max_lat);
 111     daspect([1,cos(avg_lat*pi/180),1]);
 112     hold off
 113     drawnow
 114     pause(0.1)
 115     M(step-1)=getframe(gcf);
 116     title(datestr(t));
 117 end % step
 118
 119 mov2mpeg(M,'M')