detection/patch2.m

   1 function M=patch2
   2 % to create w.mat:
   3 % run Adam's simulation, currently results in
   4 %
   5 % /share_home/akochans/WRF341F/wrf-fire/WRFV3/test/em_utfire_1d_med_4km_200m
   6 % then in Matlab
   7 % arrays needed only once
   8 % f='wrfout_d01_2013-08-20_00:00:00';
   9 % t=nc2struct(f,{'Times'},{});  n=size(t.times,2)
  10 % w=nc2struct(f,{'Times','TIGN_G','FXLONG','FXLAT','UNIT_FXLAT','UNIT_FXLONG','XLONG','XLAT','NFUEL_CAT'},{'DX','DY'},n);
  11 % save ~/w.mat w
  12 %
  13 % array at fire resolution every saved timestep
  14 % to create s.mat:
  15 % a=dir('wrfout_d01*');
  16 % s=read_wrfout_sel({a.name},{'FGRNHFX',Times});
  17 % save ~/s.mat s
  18 %
  19 % arrays at atm resolution every saved timestep
  20 % to create ss.mat
  21 % a=dir('wrfout_d01*')
  22 % s=read_wrfout_sel({a.name},{'Times','UAH','VAH'})
  23 % save ss s
  24 %
  25 % fuels.m is created by WRF-SFIRE at the beginning of the run
  26
  27 % ****** REQUIRES Matlab 2013a - will not run in earlier versions *******
  28
  29 % figures
  30 figmap=2;
  31 fig3d=0;
  32
  33 load w
  34
  35 % wind
  36 ss=load('ss');ss=ss.s;
  37 ss.steps=size(ss.times,2);
  38 ss.time=zeros(1,ss.steps);
  39 for i=1:ss.steps,
  40     ss.time(i)=datenum(char(ss.times(:,i)'));
  41 end
  42 ss.num=1:ss.steps;
  43 ss.min_time=min(ss.time);
  44 ss.max_time=max(ss.time);
  45 % interpolate surface wind to center of the grid
  46 ss.uh=0.5*(ss.uah(1:end-1,:,:)+ss.uah(2:end,:,:));
  47 ss.vh=0.5*(ss.vah(:,1:end-1,:)+ss.vah(:,2:end,:));
  48 % load s
  49 % load c
  50 fuels
  51
  52 % establish boundaries from simulations
  53
  54 sim.min_lat = min(w.fxlat(:))
  55 sim.max_lat = max(w.fxlat(:));
  56 sim.min_lon = min(w.fxlong(:));
  57 sim.max_lon = max(w.fxlong(:));
  58 sim.min_tign= min(w.tign_g(:));
  59 sim.max_tign= max(w.tign_g(:));
  60 act.x=find(w.tign_g(:)<sim.max_tign);
  61 act.min_lat = min(w.fxlat(act.x));
  62 act.max_lat = max(w.fxlat(act.x));
  63 act.min_lon = min(w.fxlong(act.x));
  64 act.max_lon = max(w.fxlong(act.x));
  65 margin=0.5;
  66 min_lon=max(sim.min_lon,act.min_lon-margin*(act.max_lon-act.min_lon));
  67 min_lat=max(sim.min_lat,act.min_lat-margin*(act.max_lat-act.min_lat));
  68 max_lon=min(sim.max_lon,act.max_lon+margin*(act.max_lon-act.min_lon));
  69 max_lat=min(sim.max_lat,act.max_lat+margin*(act.max_lat-act.min_lat));
  70
  71 default_bounds{1}=[min_lon,max_lon,min_lat,max_lat];
  72 default_bounds{2}=[sim.min_lon,sim.max_lon,sim.min_lat,sim.max_lat];
  73 for i=1:length(default_bounds),fprintf('default bounds %i: %8.5f %8.5f %8.5f %8.5f\n',i,default_bounds{i});end
  74
  75 bounds=input('enter bounds [min_lon,max_lon,min_lat,max_lat] or number of bounds above (1)> ');
  76 if isempty(bounds),bounds=1;end
  77 if length(bounds)==1,
  78     bounds=default_bounds{bounds};
  79 end
  80 [ii,jj]=find(w.fxlong>=bounds(1) & w.fxlong<=bounds(2) & w.fxlat >=bounds(3) & w.fxlat <=bounds(4));
  81 ispan=min(ii):max(ii);
  82 jspan=min(jj):max(jj);
  83 if isempty(ispan) | isempty(jspan), error('selection empty'),end
  84
  85 % restrict data for display
  86
  87 red.fxlat=w.fxlat(ispan,jspan);
  88 red.fxlong=w.fxlong(ispan,jspan);
  89 red.tign_g=w.tign_g(ispan,jspan);
  90 red.nfuel_cat=w.nfuel_cat(ispan,jspan);
  91
  92 red.min_lat = min(red.fxlat(:))
  93 red.max_lat = max(red.fxlat(:))
  94 red.min_lon = min(red.fxlong(:))
  95 red.max_lon = max(red.fxlong(:))
  96
  97 % convert tign_g to datenum
  98 w.time=datenum(char(w.times)');
  99 red.tign=(red.tign_g - max(red.tign_g(:)))/(24*60*60) + w.time;
 100 min_tign=min(red.tign(:));
 101 max_tign=max(red.tign(:));
 102
 103
 104 red.tign(find(red.tign==max_tign))=NaN; % squash the top
 105 if fig3d>0,
 106     figure(fig3d); clf
 107     hold off
 108     h=surf(red.fxlong,red.fxlat,red.tign-min_tign);
 109     xlabel('longitude'),ylabel('latitude'),zlabel('days')
 110     set(h,'EdgeAlpha',0,'FaceAlpha',0.5); % show faces only
 111     drawnow
 112 end
 113
 114 prefix='TIFs/';    % the level2 files processed by geotiff2mat.py
 115 p=sort_rsac_files(prefix);
 116 d=p.file;
 117
 118 cmap=cmapmod14;
 119 cmap2=cmap;
 120 cmap2(1:7,:)=NaN;
 121 plot_all_level2=true;
 122
 123 figure(figmap);clf
 124 iframe=1;
 125 for i=1:length(d),
 126     file=d{i};
 127     v=readmod14(prefix,file);
 128     % select fire detection within the domain
 129     xj=find(v.lon > red.min_lon & v.lon < red.max_lon);
 130     xi=find(v.lat > red.min_lat & v.lat < red.max_lat);
 131     ax=[red.min_lon red.max_lon red.min_lat red.max_lat];
 132     if ~isempty(xi) & ~isempty(xj)
 133         x=[];
 134         x.axis=[red.min_lon,red.max_lon,red.min_lat,red.max_lat];
 135         x.file=v.file;
 136         x.time=v.time;
 137         x.data=v.data(xi,xj);    % subset data
 138         x.lon=v.lon(xj);
 139         x.lat=v.lat(xi);
 140         [x.xlon,x.xlat]=meshgrid(x.lon,x.lat);
 141         if any(x.data(:)>6), % some data present - not all absent, water, or cloud
 142             figure(figmap);clf
 143             showmod14(x)
 144             hold on
 145             contour(red.fxlong,red.fxlat,red.tign,[v.time v.time],'-k');
 146             fprintf('image time            %s\n',datestr(x.time));
 147             fprintf('min wind field time   %s\n',datestr(ss.min_time));
 148             fprintf('max wind field time   %s\n',datestr(ss.max_time));
 149             if x.time>=ss.min_time && x.time <= ss.max_time,
 150                 step=interp1(ss.time,ss.num,x.time);
 151                 step0=floor(step);
 152                 if step0 < ss.steps,
 153                     step1 = step0+1;
 154                 else
 155                     step1=step0;
 156                     step0=step1-1;
 157                 end
 158                 w0=step1-step;
 159                 w1=step-step0;
 160                 uu=w0*ss.uh(:,:,step0)+w1*ss.uh(:,:,step1);
 161                 vv=w0*ss.vh(:,:,step0)+w1*ss.vh(:,:,step1);
 162                 fprintf('wind interpolated to %s from\n',datestr(x.time))
 163                 fprintf('step %i %s weight %8.3f\n',step0,datestr(ss.time(step0)),w0)
 164                 fprintf('step %i %s weight %8.3f\n',step1,datestr(ss.time(step1)),w1)
 165                 fprintf('wind interpolated to %s from\n',datestr(x.time))
 166                 sc=0.006;quiver(w.xlong,w.xlat,sc*uu,sc*vv,0);
 167             end
 168             hold off
 169             axis(ax)
 170             drawnow
 171             M(iframe)=getframe(gcf);
 172             iframe=iframe+1;
 173             print(figmap,'-dpng',['fig',v.timestr],'-r1600');
 174         end
 175         if any(x.data(:)>7) && fig3d>0,
 176             figure(fig3d)
 177             x.C2=cmap2(x.data+1,:); % translate data to RGB colormap, NaN=no detection
 178             x.C2=reshape(x.C2,[size(x.data),size(cmap2,2)]);
 179             hold on
 180             h2=surf(x.xlon,x.xlat,(v.time-min_tign)*ones(size(x.data)),x.C2);
 181             set(h2,'EdgeAlpha',0,'FaceAlpha',1); % show faces only
 182             hold off
 183             drawnow
 184         end
 185     end
 186 end
 187
 188 end