detection/barker2.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 % /home/akochans/NASA_WSU/wrf-fire/WRFV3/test/em_barker_moist/wrfoutputfiles_live_0.25
   8 % then in Matlab
   9 % f='wrfout_d05_2012-09-15_00:00:00';
  10 % t=nc2struct(f,{'Times'},{});  n=size(t.times,2);  w=nc2struct(f,{'TIGN_G','FXLONG','FXLAT','UNIT_FXLAT','UNIT_FXLONG'},{},n);
  11 % save ~/w.mat w
  12 %
  13 % to create c.mat
  14 % c=nc2struct(f,{'NFUEL_CAT'},{},1);
  15 % save ~/c.mat c
  16 %
  17 % to create s.mat:
  18 % s=read_wrfout_sel({'wrfout_d05_2012-09-09_00:00:00','wrfout_d05_2012-09-12_00:00:00','wrfout_d05_2012-09-15_00:00:00'},{'FGRNHFX'});
  19 % save ~/s.mat s
  20 %
  21 % fuels.m is created by WRF-SFIRE at the beginning of the run
  22
  23 % ****** REQUIRES Matlab 2013a - will not run in earlier versions *******
  24
  25 conus = input('enter 0 for viirs, 1 for modis> ');
  26 if conus==0,
  27         v=read_fire_kml('conus_viirs.kml');
  28         detection='VIIRS';
  29 elseif conus==1,
  30         v=read_fire_kml('conus_modis.kml');
  31         detection='MODIS';
  32 else
  33         error('need kml file')
  34 end
  35
  36 load w
  37 load s
  38 load c
  39 fuels
  40
  41 % establish boundaries from simulations
  42
  43 min_lat = min(w.fxlat(:))
  44 max_lat = max(w.fxlat(:))
  45 min_lon = min(w.fxlong(:))
  46 max_lon = max(w.fxlong(:))
  47 min_tign= min(w.tign_g(:))
  48
  49 default_bounds{1}=[min_lon,max_lon,min_lat,max_lat];
  50 default_bounds{2}=[-119.5, -119.0, 47.95, 48.15];
  51 for i=1:length(default_bounds),fprintf('default bounds %i: %8.5f %8.5f %8.5f %8.5f\n',i,default_bounds{i});end
  52
  53 bounds=input('enter bounds [min_lon,max_lon,min_lat,max_lat] or number of bounds above> ');
  54 if length(bounds)==1, bounds=default_bounds{bounds}; end
  55 [ii,jj]=find(w.fxlong>=bounds(1) & w.fxlong<=bounds(2) & w.fxlat >=bounds(3) & w.fxlat <=bounds(4));
  56 ispan=min(ii):max(ii);
  57 jspan=min(jj):max(jj);
  58
  59 % restrict data
  60 w.fxlat=w.fxlat(ispan,jspan);
  61 w.fxlong=w.fxlong(ispan,jspan);
  62 w.tign_g=w.tign_g(ispan,jspan);
  63 c.nfuel_cat=c.nfuel_cat(ispan,jspan);
  64
  65 min_lat = min(w.fxlat(:))
  66 max_lat = max(w.fxlat(:))
  67 min_lon = min(w.fxlong(:))
  68 max_lon = max(w.fxlong(:))
  69 min_tign= min(w.tign_g(:))
  70
  71 % rebase time on the largest tign_g = the time of the last frame, in days
  72
  73 last_time=datenum(char(w.times)');
  74 max_tign_g=max(w.tign_g(:));
  75
  76 tim_all = v.tim - last_time;
  77 tign= (w.tign_g - max_tign_g)/(24*60*60);
  78 min_tign= min(tign(:)); % initial ignition time
  79 tign(tign==0)=NaN;
  80
  81 % select fire detection within the domain and time
  82 bii=(v.lon > min_lon & v.lon < max_lon & v.lat > min_lat & v.lat < max_lat);
  83
  84 tim_in = tim_all(bii);
  85 u_in = unique(tim_in);
  86 fprintf('detection times from first ignition\n')
  87 for i=1:length(u_in)
  88     fprintf('%8.5f days %s UTC %3i %s detections\n',u_in(i)-min_tign,...
  89     datestr(u_in(i)+last_time),sum(tim_in==u_in(i)),detection);
  90 end
  91 detection_bounds=input('enter detection bounds as [upper,lower]: ')
  92 bi = bii & detection_bounds(1)  + min_tign <= tim_all ...
  93          & tim_all <= detection_bounds(2)  + min_tign;
  94 % now detection selected in time and space
  95 lon=v.lon(bi);
  96 lat=v.lat(bi);
  97 res=v.res(bi);
  98 tim=tim_all(bi);
  99 tim_ref = mean(tim);
 100
 101 fprintf('%i detections selected\n',sum(bi)),
 102 fprintf('mean detection time %g days from ignition %s UTC\n',...
 103     tim_ref-min_tign,datestr(tim_ref+last_time));
 104 fprintf('days from ignition  min %8.5f max %8.5f\n',min(tim)-min_tign,max(tim)-min_tign);
 105 fprintf('longitude           min %8.5f max %8.5f\n',min(lon),max(lon));
 106 fprintf('latitude            min %8.5f max %8.5f\n',min(lat),max(lat));
 107
 108 % detection selected in time and space
 109 lon=v.lon(bi);
 110 lat=v.lat(bi);
 111 res=v.res(bi);
 112 tim=tim_all(bi);
 113
 114 % plot satellite detection points
 115
 116 % plot3(lon,lat,tim,'ko'),
 117 % m=500; n=500;
 118 [m,n]=size(w.fxlong);
 119 for j=1:length(fuel), W(j)=fuel(j).weight;end
 120 for j=length(fuel)+1:max(c.nfuel_cat(:)),W(j)=NaN;end
 121 T = zeros(m,n);
 122 for j=1:n, for i=1:m
 123         T(i,j)=W(c.nfuel_cat(i,j));
 124 end,end
 125
 126 while 1
 127 tscale = 1000; % multiply fuel.weight by this to get detection time scale
 128 a=input('enter [tscale dir add1p add1m add2p add2m] (1,rad,s/m,s/m)\n');
 129 % magic match to u(1) and viirs seems [1000 3*pi/4 -0.0002  -0.0001 0 0.00005]
 130 tscale=a(1)
 131 if tscale<=0, break, end
 132 dir=a(2)
 133 add1p=a(3)
 134 add1m=a(4)
 135 add2p=a(5)
 136 add2m=a(6)
 137
 138 v1=[cos(dir),sin(dir)]; % main direction
 139 v2=[cos(dir+pi/2),sin(dir+pi/2)]; % secondary direction
 140
 141 % find ignition point
 142 [i_ign,j_ign]=find(w.tign_g == min(w.tign_g(:)));
 143
 144 % vector field (x,y) - (x_ign,y_ign)
 145 VDx=(w.fxlong-w.fxlong(i_ign,j_ign))*w.unit_fxlong;
 146 VDy=(w.fxlat-w.fxlat(i_ign,j_ign))*w.unit_fxlat;
 147
 148 p1 = VDx*v1(1)+VDy*v1(2); % projections on main direction
 149 p2 = VDx*v2(1)+VDy*v2(2); % projections on secondary direction
 150 [theta,rho]=cart2pol(p1,p2); % in polar coordinates, theta between [-pi,pi]
 151 thetas=pi*[-3/2,-1,-1/2,0,1/2,1,3/2];
 152 deltas=[add2p add1m add2m add1p add2p add1m add2m];
 153 delta = interp1(thetas,deltas,theta,'pchip').*rho;
 154
 155 tign_mod = tign + delta;
 156
 157 % probability of detection, assuming selected times are close
 158 pmap = p_map(tscale*T/(24*60*60),tign_mod - tim_ref);
 159
 160 [mm,nn]=size(w.fxlong);
 161 mi=1:ceil(mm/m):mm;
 162 ni=1:ceil(nn/n):nn;
 163 mesh_fxlong=w.fxlong(mi,ni);
 164 mesh_fxlat=w.fxlat(mi,ni);
 165 mesh_fgrnhfx=s.fgrnhfx(mi,ni,:);
 166 mesh_pmap = pmap(mi,ni);
 167 mesh_lon = mesh_fxlong(mi,ni);
 168 mesh_lat = mesh_fxlat(mi,ni);
 169 mesh_tign= tign(mi,ni);
 170
 171 % draw the fireline
 172 figure(1)
 173 hold off, clf
 174 %contour(mesh_fxlong,mesh_fxlat,mesh_tign,[tim_ref,tim_ref]);
 175
 176 % add the detection probability map
 177 %hold on
 178 h=pcolor(mesh_fxlong,mesh_fxlat,pmap);
 179 set(h,'EdgeAlpha',0,'FaceAlpha',0.5);
 180 colorbar
 181 cc=colormap; cc(1,:)=1; colormap(cc);
 182
 183 hold on
 184 plot(w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign),'xk')
 185 fprintf('first ignition at %g %g, marked by black x\n',w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign))
 186 drawnow
 187 pause(1)
 188
 189 % plot detection squares
 190
 191 C=0.5*ones(1,length(res));
 192 rlon=0.5*res/w.unit_fxlong;
 193 rlat=0.5*res/w.unit_fxlat;
 194 X=[lon-rlon,lon+rlon,lon+rlon,lon-rlon]';
 195 Y=[lat-rlat,lat-rlat,lat+rlat,lat+rlat]';
 196 hold on
 197 hh=fill(X,Y,C);
 198 hold off
 199
 200 grid,drawnow
 201
 202 % hold on, (mesh_lon,mesh_lat,mesh_tim2),grid on
 203 assim_string='';if any(delta(:)),assim_string='assimilation',end
 204 daspect([w.unit_fxlat,w.unit_fxlong,1])                % axis aspect ratio length to length
 205 title(sprintf('Barker Canyon fire %s %s %s UTC',...
 206     assim_string, detection, datestr(tim_ref+last_time)))
 207 ylabel('latitude')
 208 xlabel('longitude')
 209 drawnow
 210
 211 % evaluate likelihood
 212 ndetect=length(res);
 213 likelihood=zeros(1,ndetect);
 214 mask=cell(1,ndetect);
 215 for i=1:ndetect
 216     mask{i}=lon(i)-rlon(i) <= w.fxlong & w.fxlong <= lon(i)+rlon(i) & ...
 217             lat(i)-rlat(i) <= w.fxlat  & w.fxlat  <= lat(i)+rlat(i);
 218     likelihood(i)=ssum(pmap(mask{i}))/ssum(mask{i});
 219 end
 220 likelihood,
 221 total_likelihood=sum(likelihood) % should be really a product if they are independeny
 222                                  % but the sum seems to work little better
 223 end