detection/detect_fit.m

   1 function analysis=detect_fit
   2 % from a copy of barker2
   3
   4 disp('input data')
   5     % to create conus.kml:
   6     % download http://firemapper.sc.egov.usda.gov/data_viirs/kml/conus_hist/conus_20120914.kmz
   7     % and gunzip
   8     %
   9     % to create w.mat:
  10     % run Adam's simulation, currently results in
  11     % /home/akochans/NASA_WSU/wrf-fire/WRFV3/test/em_barker_moist/wrfoutputfiles_live_0.25
  12     % then in Matlab
  13     % f='wrfout_d05_2012-09-15_00:00:00';
  14     % t=nc2struct(f,{'Times'},{'DX','DY'});  n=size(t.times,2);  w=nc2struct(f,{'TIGN_G','FXLONG','FXLAT','UNIT_FXLAT','UNIT_FXLONG'},{},n);
  15     % save ~/w.mat w
  16     %
  17     % to create c.mat
  18     % c=nc2struct(f,{'NFUEL_CAT'},{},1);
  19     % save ~/c.mat c
  20     %
  21     % to create s.mat:
  22     % 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'});
  23     % save ~/s.mat 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
  30     conus = input_num('0 for viirs, 1 for modis',0);
  31     if conus==0,
  32             v=read_fire_kml('conus_viirs.kml');
  33             detection='VIIRS';
  34     elseif conus==1,
  35             v=read_fire_kml('conus_modis.kml');
  36             detection='MODIS';
  37     else
  38             error('need kml file')
  39     end
  40
  41     a=load('w');w=a.w;
  42     if ~isfield('dx',w),
  43         w.dx=444.44;
  44         w.dy=444.44;
  45         warning('fixing up w for old w.mat file from Barker fire')
  46     end
  47
  48     a=load('s');s=a.s;
  49     a=load('c');c=a.c;
  50     fuel.weight=0; % just to let Matlab know what fuel is going to be at compile time
  51     fuels
  52
  53
  54 disp('subset and process inputs')
  55
  56     % establish boundaries from simulations
  57
  58     min_lat = min(w.fxlat(:))
  59     max_lat = max(w.fxlat(:))
  60     min_lon = min(w.fxlong(:))
  61     max_lon = max(w.fxlong(:))
  62     min_tign= min(w.tign_g(:))
  63
  64     default_bounds{1}=[min_lon,max_lon,min_lat,max_lat];
  65     default_bounds{2}=[-119.5, -119.0, 47.95, 48.15];
  66     display_bounds=default_bounds{2};
  67     for i=1:length(default_bounds),fprintf('default bounds %i: %8.5f %8.5f %8.5f %8.5f\n',i,default_bounds{i});end
  68
  69     bounds=input_num('bounds [min_lon,max_lon,min_lat,max_lat] or number of bounds above',1);
  70     if length(bounds)==1, bounds=default_bounds{bounds}; end
  71     [ii,jj]=find(w.fxlong>=bounds(1) & w.fxlong<=bounds(2) & w.fxlat >=bounds(3) & w.fxlat <=bounds(4));
  72     ispan=min(ii):max(ii);
  73     jspan=min(jj):max(jj);
  74
  75     % restrict data
  76     w.fxlat=w.fxlat(ispan,jspan);
  77     w.fxlong=w.fxlong(ispan,jspan);
  78     w.tign_g=w.tign_g(ispan,jspan);
  79     c.nfuel_cat=c.nfuel_cat(ispan,jspan);
  80
  81     min_lat = min(w.fxlat(:))
  82     max_lat = max(w.fxlat(:))
  83     min_lon = min(w.fxlong(:))
  84     max_lon = max(w.fxlong(:))
  85     min_lon = display_bounds(1);
  86     max_lon = display_bounds(2);
  87     min_lat = display_bounds(3);
  88     max_lat = display_bounds(4);
  89
  90     min_tign= min(w.tign_g(:))
  91
  92     % rebase time on the largest tign_g = the time of the last frame, in days
  93
  94     last_time=datenum(char(w.times)');
  95     max_tign_g=max(w.tign_g(:));
  96
  97     tim_all = v.tim - last_time;
  98     tign= (w.tign_g - max_tign_g)/(24*60*60);  % now tign is in days
  99     min_tign= min(tign(:)); % initial ignition time
 100     tign_disp=tign;
 101     tign_disp(tign==0)=NaN;      % for display
 102
 103     % select fire detection within the domain and time
 104     bii=(v.lon > min_lon & v.lon < max_lon & v.lat > min_lat & v.lat < max_lat);
 105
 106     tim_in = tim_all(bii);
 107     u_in = unique(tim_in);
 108     fprintf('detection times from first ignition\n')
 109     for i=1:length(u_in)
 110         detection_freq(i)=sum(tim_in==u_in(i));
 111         fprintf('%8.5f days %s UTC %3i %s detections\n',u_in(i)-min_tign,...
 112         datestr(u_in(i)+last_time),detection_freq(i),detection);
 113     end
 114     [max_freq,i]=max(detection_freq);
 115     tol=0.01;
 116     detection_bounds=input_num('detection bounds as [upper,lower]',...
 117         [u_in(i)-min_tign-tol,u_in(i)-min_tign+tol]);
 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     detection_days_from_ignition=tim_ref-min_tign;
 129     detection_datestr=datestr(tim_ref+last_time);
 130     fprintf('mean detection time %g days from ignition %s UTC\n',...
 131         detection_days_from_ignition,detection_datestr);
 132     fprintf('days from ignition  min %8.5f max %8.5f\n',min(tim)-min_tign,max(tim)-min_tign);
 133     fprintf('longitude           min %8.5f max %8.5f\n',min(lon),max(lon));
 134     fprintf('latitude            min %8.5f max %8.5f\n',min(lat),max(lat));
 135
 136     % detection selected in time and space
 137     lon=v.lon(bi);
 138     lat=v.lat(bi);
 139     res=v.res(bi);
 140     tim=tim_all(bi);
 141
 142     % set up reduced resolution plots
 143     [m,n]=size(w.fxlong);
 144     m_plot=m; n_plot=n;
 145
 146     m1=map_index(display_bounds(1),bounds(1),bounds(2),m);
 147     m2=map_index(display_bounds(2),bounds(1),bounds(2),m);
 148     n1=map_index(display_bounds(3),bounds(3),bounds(4),n);
 149     n2=map_index(display_bounds(4),bounds(3),bounds(4),n);
 150     mi=m1:ceil((m2-m1+1)/m_plot):m2; % reduced index vectors
 151     ni=n1:ceil((n2-n1+1)/n_plot):n2;
 152     mesh_fxlong=w.fxlong(mi,ni);
 153     mesh_fxlat=w.fxlat(mi,ni);
 154     [mesh_m,mesh_n]=size(mesh_fxlat);
 155
 156     % find ignition point
 157     [i_ign,j_ign]=find(w.tign_g == min(w.tign_g(:)));
 158     if length(i_ign)~=1,error('assuming single ignition point here'),end
 159
 160     % set up constraint on ignition point being the same
 161     Constr_ign = zeros(m,n); Constr_ign(i_ign,j_ign)=1;
 162
 163     detection_mask=zeros(m,n);
 164     detection_time=tim_ref*ones(m,n);
 165
 166     % resolution diameter in longitude/latitude units
 167     rlon=0.5*res/w.unit_fxlong;
 168     rlat=0.5*res/w.unit_fxlat;
 169
 170     lon1=lon-rlon;
 171     lon2=lon+rlon;
 172     lat1=lat-rlat;
 173     lat2=lat+rlat;
 174     for i=1:length(lon),
 175         square = w.fxlong>=lon1(i) & w.fxlong<=lon2(i) & ...
 176                  w.fxlat >=lat1(i) & w.fxlat <=lat2(i);
 177         detection_mask(square)=1;
 178     end
 179     C=0.5*ones(1,length(res));
 180     X=[lon-rlon,lon+rlon,lon+rlon,lon-rlon]';
 181     Y=[lat-rlat,lat-rlat,lat+rlat,lat+rlat]';
 182     plotstate(1,detection_mask,['Fire detection at ',detection_datestr],[])
 183     % add ignition point
 184     hold on, plot(w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign),'xw'); hold off
 185     % legend('first ignition at %g %g',w.fxlong(i_ign,j_ign),w.fxlat(i_ign,j_ign))
 186
 187     fuelweight(length(fuel)+1:max(c.nfuel_cat(:)))=NaN;
 188     for j=1:length(fuel),
 189         fuelweight(j)=fuel(j).weight;
 190     end
 191     W = zeros(m,n);
 192     for j=1:n, for i=1:m
 193            W(i,j)=fuelweight(c.nfuel_cat(i,j));
 194     end,end
 195
 196     plotstate(2,W,'Fuel weight',[])
 197
 198 disp('optimization loop')
 199 h =zeros(m,n); % initial increment
 200 plotstate(3,tign,'Forecast fire arrival time',detection_time(1));
 201 for istep=1:5
 202
 203     % can change the objective function here
 204     alpha=input_num('penalty coefficient alpha, <0 to end',1e-2);
 205     if alpha<0, break, end
 206     % TC = W/(900*24); % time constant = fuel gone in one hour
 207     TC = 1/24;  % detection time constants in hours
 208     stretch=input_num('Tmin,Tmax,Tneg,Tpos',[0.5,10,5,10]);
 209     nodetw=input_num('no fire detection weight',0.1);
 210     power=input_num('negative laplacian power',0.51);
 211
 212     psi = detection_mask - nodetw*(1-detection_mask);
 213
 214     [Js,search]=objective(tign,h);
 215     search = -search/big(search); % initial search direction
 216
 217     plotstate(4,search,'Search direction',0);
 218     h=zeros(m,n); % initial increment
 219     stepsize=0;
 220     % initial estimate of stepsize
 221     last_stepsize = 3;
 222     for i=2:100 % crude manual line search
 223         s=input_num('step size',last_stepsize);
 224         stepsize(i)=s;
 225         last_stepsize=s;
 226         plotstate(5,tign+h+last_stepsize*search,'Line search',detection_time(1));
 227         [Js(i),delta]=objective(tign,h+last_stepsize*search,'noplot');
 228         c=input_num('try another step size: 0/1',1)
 229         if c==0, break, end
 230     end
 231     h = h + last_stepsize*search;
 232     plotstate(6,tign+h,sprintf('Analysis descent iteration %i',istep),detection_time(1));
 233 end
 234 disp('converting analysis fire arrival time from days with zero at the end of the fire to original scale')
 235 analysis=max_tign_g+(24*60*60)*(tign+h);
 236 disp('input the analysis as tign in WRF-SFIRE with fire_perimeter_time=detection time')
 237
 238     function [J,delta]=objective(tign,h,noplot)
 239     % compute objective function and optionally ascent direction
 240     T=tign+h;
 241     [f0,f1]=like1(psi,detection_time-T,TC*stretch);
 242     F = f1;             % forcing
 243     % objective function and preconditioned gradient
 244     Ah = poisson_fft2(h,[w.dx,w.dy],1);
 245     J = alpha*0.5*(h(:)'*Ah(:)) - ssum(psi.*f0)/(m*n);
 246     fprintf('Objective function J=%g\n',J);
 247     gradJ = alpha*Ah + F;
 248     if ~exist('noplot','var'),
 249         plotstate(7,f0,'Detection likelihood',0.5,'-w');
 250         plotstate(8,f1,'Detection likelihood derivative',0);
 251         plotstate(9,F,'Forcing',0);
 252         plotstate(10,gradJ,'gradient of J',0);
 253     end
 254     delta = solve_saddle(Constr_ign,h,F,0,@(u) poisson_fft2(u,[w.dx,w.dy],-power)/alpha);
 255     % plotstate(11,delta,'Preconditioned gradient',0);
 256     fprintf('norm(grad(J))=%g norm(delta)=%g\n',norm(gradJ,'fro'),norm(delta,'fro'))
 257     end
 258
 259     function plotstate(fig,T,s,c,linespec)
 260             fprintf('Figure %i %s\n',fig,s)
 261             plotmap(fig,mesh_fxlong,mesh_fxlat,T(mi,ni),' ');
 262             hold on
 263             hh=fill(X,Y,C,'EdgeAlpha',1,'FaceAlpha',0);
 264             if ~exist('c','var') || isempty(c) || isnan(c),
 265                 title(s);
 266             else
 267                 title(sprintf('%s, contour=%g',s,c(1)))
 268                 if ~exist('linespec','var'),
 269                     linespec='-k';
 270                 end
 271                 contour(mesh_fxlong,mesh_fxlat,T(mi,ni),[c c],linespec)
 272             end
 273             hold off
 274             ratio=[w.unit_fxlat,w.unit_fxlong];
 275             xlabel longtitude
 276             ylabel latitude
 277             ratio=[ratio/norm(ratio),1];
 278             daspect(ratio)
 279             axis tight
 280             drawnow
 281     end
 282
 283 end % detect_fit
 284
 285 function i=map_index(x,a,b,n)
 286 % find image of x under linear map [a,b] -> [1,m]
 287 % and round to integer
 288 i=round(1+(n-1)*(x-a)/(b-a));
 289 end