detection/patch3.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=2;
  36 fig3d=0;
  37
  38 % convert tign_g to datenum
  39 w.time=datenum(char(w.times)');
  40 red.tign=(red.tign_g - max(red.tign_g(:)))/(24*60*60) + w.time;
  41 min_tign=min(red.tign(:));
  42 max_tign=max(red.tign(:));
  43
  44 red.tign(find(red.tign==max_tign))=NaN; % squash the top
  45 if fig3d>0,
  46     figure(fig3d); clf
  47     hold off
  48     h=surf(red.fxlong,red.fxlat,red.tign-min_tign);
  49     xlabel('longitude'),ylabel('latitude'),zlabel('days')
  50     set(h,'EdgeAlpha',0,'FaceAlpha',0.5); % show faces only
  51     drawnow
  52 end
  53
  54 cmap=cmapmod14;
  55 cmap2=cmap;
  56 cmap2(1:7,:)=NaN;
  57 plot_all_level2=true;
  58
  59 figure(figmap);clf
  60 iframe=1;
  61 for i=1:length(r.x),
  62     x=r.x{i};
  63     [x.xlon,x.xlat]=meshgrid(x.lon,x.lat);
  64     if any(x.data(:)>6), % some data present - not all absent, water, or cloud
  65             figure(figmap);clf
  66             showmod14(x)
  67             hold on
  68             contour(red.fxlong,red.fxlat,red.tign,[x.time x.time],'-r');
  69             fprintf('image time            %s\n',datestr(x.time));
  70             fprintf('min wind field time   %s\n',datestr(ss.min_time));
  71             fprintf('max wind field time   %s\n',datestr(ss.max_time));
  72             if x.time>=ss.min_time && x.time <= ss.max_time,
  73                 step=interp1(ss.time,ss.num,x.time);
  74                 step0=floor(step);
  75                 if step0 < ss.steps,
  76                     step1 = step0+1;
  77                 else
  78                     step1=step0;
  79                     step0=step1-1;
  80                 end
  81                 w0=step1-step;
  82                 w1=step-step0;
  83                 uu=w0*ss.uh(:,:,step0)+w1*ss.uh(:,:,step1);
  84                 vv=w0*ss.vh(:,:,step0)+w1*ss.vh(:,:,step1);
  85                 fprintf('wind interpolated to %s from\n',datestr(x.time))
  86                 fprintf('step %i %s weight %8.3f\n',step0,datestr(ss.time(step0)),w0)
  87                 fprintf('step %i %s weight %8.3f\n',step1,datestr(ss.time(step1)),w1)
  88                 fprintf('wind interpolated to %s from\n',datestr(x.time))
  89                 sc=0.006;quiver(w.xlong,w.xlat,sc*uu,sc*vv,0);
  90             end
  91             hold off
  92             axis(red.axis)
  93             drawnow
  94             M(i)=getframe(gcf);
  95     end
  96     if any(x.data(:)>7) && fig3d>0,
  97             figure(fig3d)
  98             x.C2=cmap2(x.data+1,:); % translate data to RGB colormap, NaN=no detection
  99             x.C2=reshape(x.C2,[size(x.data),size(cmap2,2)]);
 100             hold on
 101             h2=surf(x.xlon,x.xlat,(v.time-min_tign)*ones(size(x.data)),x.C2);
 102             set(h2,'EdgeAlpha',0,'FaceAlpha',1); % show faces only
 103             hold off
 104             drawnow
 105     end
 106 end
 107
 108 return
 109
 110 tim_in = tim_all(bii);
 111 u_in = unique(tim_in);
 112 fprintf('detection times from first ignition\n')
 113 for i=1:length(u_in)
 114     fprintf('%8.5f days %s UTC %3i %s detections\n',u_in(i)-min_tign,...
 115     datestr(u_in(i)+w.time),sum(tim_in==u_in(i)),detection);
 116 end
 117 detection_bounds=input('enter detection bounds as [upper,lower]: ')
 118 bi = bii & detection_bounds(1)  + min_tign <= tim_all ...
 119          & tim_all <= detection_bounds(2)  + min_tign;
 120 % now detection selected in time and space
 121 lon=v.lon(bi);
 122 lat=v.lat(bi);
 123 res=v.res(bi);
 124 tim=tim_all(bi);
 125 tim_ref = mean(tim);
 126
 127 fprintf('%i detections selected\n',sum(bi)),
 128 fprintf('mean detection time %g days from ignition %s UTC\n',...
 129     tim_ref-min_tign,datestr(tim_ref+w.time));
 130 fprintf('days from ignition  min %8.5f max %8.5f\n',min(tim)-min_tign,max(tim)-min_tign);
 131 fprintf('longitude           min %8.5f max %8.5f\n',min(lon),max(lon));
 132 fprintf('latitude            min %8.5f max %8.5f\n',min(lat),max(lat));
 133
 134 % detection selected in time and space
 135 lon=v.lon(bi);
 136 lat=v.lat(bi);
 137 res=v.res(bi);
 138 tim=tim_all(bi);
 139
 140 % plot satellite detection points
 141
 142 % plot3(lon,lat,tim,'ko'),
 143 % m=500; n=500;
 144 [m,n]=size(w.fxlong);
 145 for j=1:length(fuel), W(j)=fuel(j).weight;end
 146 for j=length(fuel)+1:max(c.nfuel_cat(:)),W(j)=NaN;end
 147 T = zeros(m,n);
 148 for j=1:n, for i=1:m
 149         T(i,j)=W(c.nfuel_cat(i,j));
 150 end,end
 151
 152 while 1
 153 tscale = 1000; % multiply fuel.weight by this to get detection time scale
 154 a=input('enter [tscale dir add1p add1m add2p add2m] (1,rad,s/m,s/m)\n');
 155 % magic match to u(1) and viirs seems [1000 3*pi/4 -0.0002  -0.0001 0 0.00005]
 156 tscale=a(1)
 157 if tscale<=0, break, end
 158 dir=a(2)
 159 add1p=a(3)
 160 add1m=a(4)
 161 add2p=a(5)
 162 add2m=a(6)
 163
 164 v1=[cos(dir),sin(dir)]; % main direction
 165 v2=[cos(dir+pi/2),sin(dir+pi/2)]; % secondary direction
 166
 167 % find ignition point
 168 [i_ign,j_ign]=find(w.tign_g == min(w.tign_g(:)));
 169
 170 % vector field (x,y) - (x_ign,y_ign)
 171 VDx=(w.fxlong-w.fxlong(i_ign,j_ign))*w.unit_fxlong;
 172 VDy=(w.fxlat-w.fxlat(i_ign,j_ign))*w.unit_fxlat;
 173
 174 p1 = VDx*v1(1)+VDy*v1(2); % projections on main direction
 175 p2 = VDx*v2(1)+VDy*v2(2); % projections on secondary direction
 176 [theta,rho]=cart2pol(p1,p2); % in polar coordinates, theta between [-pi,pi]
 177 thetas=pi*[-3/2,-1,-1/2,0,1/2,1,3/2];
 178 deltas=[add2p add1m add2m add1p add2p add1m add2m];
 179 delta = interp1(thetas,deltas,theta,'pchip').*rho;
 180
 181 tign_mod = tign + delta;
 182
 183 % probability of detection, assuming selected times are close
 184 pmap = p_map(tscale*T/(24*60*60),tign_mod - tim_ref);
 185
 186 [mm,nn]=size(w.fxlong);
 187 mi=1:ceil(mm/m):mm;
 188 ni=1:ceil(nn/n):nn;
 189 mesh_fxlong=w.fxlong(mi,ni);
 190 mesh_fxlat=w.fxlat(mi,ni);
 191 mesh_fgrnhfx=s.fgrnhfx(mi,ni,:);
 192 mesh_pmap = pmap(mi,ni);
 193 mesh_lon = mesh_fxlong(mi,ni);
 194 mesh_lat = mesh_fxlat(mi,ni);
 195 mesh_tign= tign(mi,ni);
 196
 197 % draw the fireline
 198 figure(1)
 199 hold off, clf
 200 %contour(mesh_fxlong,mesh_fxlat,mesh_tign,[tim_ref,tim_ref]);
 201
 202 % add the detection probability map
 203 %hold on
 204 h=pcolor(mesh_fxlong,mesh_fxlat,pmap);
 205 set(h,'EdgeAlpha',0,'FaceAlpha',0.5);
 206 colorbar
 207 cc=colormap; cc(1,:)=1; colormap(cc);
 208
 209 hold on
 210 plot(w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign),'xk')
 211 fprintf('first ignition at %g %g, marked by black x\n',w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign))
 212 drawnow
 213 pause(1)
 214
 215 % plot detection squares
 216
 217 C=0.5*ones(1,length(res));
 218 rlon=0.5*res/w.unit_fxlong;
 219 rlat=0.5*res/w.unit_fxlat;
 220 X=[lon-rlon,lon+rlon,lon+rlon,lon-rlon]';
 221 Y=[lat-rlat,lat-rlat,lat+rlat,lat+rlat]';
 222 hold on
 223 hh=fill(X,Y,C);
 224 hold off
 225
 226 grid,drawnow
 227
 228 % hold on, (mesh_lon,mesh_lat,mesh_tim2),grid on
 229 assim_string='';if any(delta(:)),assim_string='assimilation',end
 230 daspect([w.unit_fxlat,w.unit_fxlong,1])                % axis aspect ratio length to length
 231 title(sprintf('Barker Canyon fire %s %s %s UTC',...
 232     assim_string, detection, datestr(tim_ref+w.time)))
 233 ylabel('latitude')
 234 xlabel('longitude')
 235 drawnow
 236
 237 % evaluate likelihood
 238 ndetect=length(res);
 239 likelihood=zeros(1,ndetect);
 240 mask=cell(1,ndetect);
 241 for i=1:ndetect
 242     mask{i}=lon(i)-rlon(i) <= w.fxlong & w.fxlong <= lon(i)+rlon(i) & ...
 243             lat(i)-rlat(i) <= w.fxlat  & w.fxlat  <= lat(i)+rlat(i);
 244     likelihood(i)=ssum(pmap(mask{i}))/ssum(mask{i});
 245 end
 246 likelihood,
 247 total_likelihood=sum(likelihood) % should be really a product if they are independeny
 248                                  % but the sum seems to work little better
 249 end