other/Matlab/impact/burned_by_category.m

   1 function area=burned_by_category(t,burn_time)
   2 % out=burned_by_category(t,burn_time)
   3 % Input:
   4 %    t          structure created from wrfout as follows:
   5 %               t=read_wrfout_tign('wrfout file name')
   6 %               save t, copy t.mat to another computer if needed, load t
   7 %    burn_time  date string in a format understood by Matlab datenum function
   8 %               assumed in UTC time zone
   9 % Output:
  10 %    area(k)  (m^2) area burned in category k until burn_time
  11 %
  12 % Usage:
  13 %    first set search path:
  14 %    clone the wrf-fire git repository
  15 %    cd wrf-fire/other/Matlab
  16 %    startup
  17 %    cd to your working directory
  18 %    load t
  19 %    area=burned_by_category(t,'2011-06-29_00:00:00')
  20
  21
  22 end_times=t.times;
  23 end_datenum = datenum(end_times); % in days from some instant in the past
  24 end_minutes=t.xtime(end); % from simulation start
  25 start_datenum=end_datenum-end_minutes/(24*60);
  26 burn_datenum=datenum(burn_time);
  27 burn_datestr=datestr(burn_datenum,'dd-mmm-yyyy HH:MM');
  28 burn_seconds=(burn_datenum-start_datenum)*24*60*60;
  29 da=t.dx*t.dy/prod(size(t.fxlong)./size(t.xlong));
  30 acre=4046.872609874252; % convert from m^2 to ac
  31
  32 fprintf('Simulation start %s\n',datestr(start_datenum));
  33 fprintf('Simulation end   %s\n',datestr(end_datenum));
  34 fprintf('Burn cut off     %s = %20g from sim start\n',burn_datestr, burn_seconds);
  35 fprintf('Fire mesh cell %g m^2\n',da);
  36
  37 cats = t.nfuel_cat .* (t.tign_g <= burn_seconds);
  38 % now cats(i,i) is cat number if cell burned, 0 if not
  39 num_cats=max(t.nfuel_cat(:));
  40 for i=1:num_cats
  41     count(i)=sum(cats(:)==i);
  42 end
  43
  44 area = count*da;
  45
  46 % pic
  47 cats(cats==0|cats==14)=NaN;
  48 clf,hold off
  49 mesh(t.fxlong,t.fxlat,cats);
  50 view(2)
  51 colorbar
  52 title(['Area burned by ',burn_datestr,' UTC by fuel category'])
  53 xlabel('Longitude')
  54 ylabel('Latitude')
  55
  56 end