chem/wrfout2bluesky_new.m

   1 function wrfout2bluesky(in)
   2 % wrfout2bluesky(in,out)
   3 % extract from wrfout the fire area info for bluesky
   4 % arguments
   5 %   in  wrfout file name, or cell array of file names
   6 %   out bluesky csv file name
   7 % see http://plone.airfire.org/bluesky/framework/comma-separated-value-csv-files
   8
   9 % Jan Mandel, July/August 2013
  10
  11 if ~exist('in','var')
  12     default_in={'wrfout_d05_2012-09-09_00:00:00','wrfout_d05_2012-09-12_00:00:00','wrfout_d05_2012-09-15_00:00:00'};
  13     disp('default wrfout input:')
  14     for i=1:length(default_in)
  15         disp(default_in{i})
  16     end
  17     in=input('enter wrfout file name: ','s');
  18     if isempty(in),
  19         in=default_in;
  20     end
  21 end
  22 if ~iscell(in),
  23     in = {in};
  24 end
  25 w=nc2struct(in{1},{'XLONG','XLAT','FXLONG','FXLAT','NFUEL_CAT','FGIP','FIRE_AREA'},...
  26     {'DX','DY'},1); % get what does not change from the first netcdf file
  27 dx=w.dx;
  28 dy=w.dy;
  29 [m,n]=size(w.xlong);
  30 [fm,fn]=size(w.fxlong);
  31 srx=fm/m;
  32 sry=fn/n;
  33 fdx=dx/srx;
  34 fdy=dy/sry;
  35 fprintf('mesh size %i by %i\n',m,n)
  36 fprintf('fire mesh size %i by %i\n',fm,fn)
  37 fprintf('refinement ratio %g %g\n',srx,sry)
  38 fprintf('fire meshstep %g %g\n',fdx,fdy);
  39 lon=w.fxlong(:,:,1);lon=lon(:);
  40 lat=w.fxlat(:,:,1);lat=lat(:);
  41 min_x=min(lon);
  42 max_x=max(lon);
  43 min_y=min(lat);
  44 max_y=max(lat);
  45 lon_ctr=mean(lon);
  46 lat_ctr=mean(lat);
  47 fprintf('coordinates %g to %g by %g to %g\n',min_x,max_x,min_y,max_y);
  48 fprintf('the center of the domain is at coordinates %g %g\n',lon_ctr,lat_ctr)
  49 %emission_factors=read_namelist_fire_emissions('namelist.fire_emissions_bluesky') % emissions factors table
  50 %species=fields(emission_factors);
  51 species={};
  52 nspecies=length(species);
  53 lat=w.fxlat(:,:,1);
  54 lon=w.fxlong(:,:,1);
  55 fgip=w.fgip;
  56 nfuel_cat=w.nfuel_cat;
  57
  58
  59 loc=open_out('fire_locations.csv');
  60 hou=open_out('fire_hourly.csv');
  61
  62 % header
  63 header_loc='id,hour,latitude,longitude,date_time,area,heat,type,moisture_1hr,moisture_10hr,moisture_100hr,moisture_1000hr,moisture_live';
  64 fprintf(loc,header_loc);
  65 header_hou='fire_id,hour,ignition_date_time,area_fract,flame_profile,smolder_profile,residual_profile';
  66 fprintf(hou,header_hou);
  67 for i=1:nspecies
  68     fprintf(hou,',%s',species{i});
  69 end
  70 fprintf(loc,'\r\n');
  71 fprintf(hou,'\r\n');
  72
  73 % initialize loop over frames
  74 if any(w.fire_area(:)), error('must start from no fire state'),end
  75 area_old=zeros(size(w.fire_area));
  76 mass_old=area_old;
  77 fire_id=0;
  78 for ifile=1:length(in),
  79     fprintf('reading file %s\n',in{ifile});
  80     w=nc2struct(in{ifile},{'Times'},{}); % read netcdf file
  81     t=char(w.times');  % convert times to character matrix one string per row
  82     s=size(t,1);
  83     fprintf('timesteps %i meshstep %g %g\n',s,dx,dy);
  84     for istep=1:s
  85         d=datenum(t(istep,:),'yyyy-mm-dd_HH:MM:SS'); % convert ESMF time string to double
  86         tim=[datestr(d,'yyyymmddHHMM'),'Z'];         % convert to bluesky time format
  87         v = datevec(d);
  88         if(v(5)==0 & v(6) == 0) % whole hour
  89         % if(v(6) == 0) % whole minute
  90             w=nc2struct(in{ifile},{'FIRE_AREA','FUEL_FRAC','FMC_GC'},{},istep); % read step istep
  91             area=w.fire_area*fdx*fdy;               % fire area (m^2)
  92             if any(area(:)< area_old(:) - fdx*fdy*eps(single(1)))
  93                 warning('fire area can only increase')
  94             end
  95             area = max(area,area_old);
  96             ta=sum(area(:));
  97             mass=fgip.*(1-w.fuel_frac)*fdx*fdy;   % fuel mass burned
  98             fprintf('timestep %i/%i %s total fire area %g m^2 fuel mass burned %g kg\n',...
  99                 istep,s,t(istep,:),ta,sum(m(:)))
 100             area_diff = area - area_old;
 101             mass_diff = mass - mass_old;
 102             if(ta>0),
 103                  cc=bwconncomp(area>0); % find connected components, requires image processing toolbox
 104                  % cc.PixelIdxList{1}=[1:length(lon(:))]'; % just take all
 105                 for id=1:length(cc.PixelIdxList)
 106                     sub=cc.PixelIdxList{id};      % subset index list
 107                     fire_id = fire_id+1;
 108                     fire_ctr_lon=mean(lon(sub));  % center longitude of the fire area
 109                     fire_ctr_lat=mean(lat(sub));  % center latitude of the fire area
 110                     a_diff_acres=sum(area_diff(sub))/4046.86; % newly burning area (acres)
 111                     % m_diff_tons=sum(mass_diff(sub))/907.185;  % newly burned fuel mass (tons)
 112                     h_diff_btus=sum(mass_diff(sub))*17.433e+06*4.30e-04; % newly generated heat (BTUs)
 113                     sub_coarse=map_submesh(sub,[m,n],[fm,fn]);
 114                     for m_class=1:5
 115                         fmc=w.fmc_gc(:,:,m_class);
 116                         m_mean(m_class)=mean(fmc(sub_coarse))*100;
 117                     end
 118                     fmt_loc='%i,%g,%08g,%08g,%s,%g,%g,RX,%g,%g,%g,%g,%g';
 119                     arg_loc={fire_id,0,fire_ctr_lat,fire_ctr_lon,tim,...
 120                         a_diff_acres,h_diff_btus,m_mean};
 121                     fprintf(loc,fmt_loc,arg_loc{:});
 122                     fmt_hou='%i,%i,%s,%g,%g,%g,%g';
 123                     arg_hou={fire_id,0,tim,1, 0.9, 0.1, 0};
 124                     fprintf(hou,fmt_hou,arg_hou{:});
 125                     nfuel_cat_sub=nfuel_cat(sub);
 126                     for i=1:nspecies
 127                         factors=emission_factors.(species{i})';
 128                         emiss_diff=m_diff(sub).*factors(nfuel_cat_sub); % emissions = mass burned * emission factor (fuel category), on the subset
 129                         emiss_diff_tons = sum(emiss_diff) /  907.185e3; % total emission of this species in tons, in the subset
 130                         fprintf(hou,',%g',emiss_diff_tons);
 131                     end
 132                     fprintf(loc,'\r\n');
 133                     fprintf(hou,'\r\n');
 134                 end
 135             end
 136             a_old=area;
 137         else
 138             fprintf('skipping %s\n',t(istep,:))
 139         end
 140     end
 141 end
 142 fclose(loc);
 143 fclose(hou);
 144 end
 145
 146 function fid=open_out(out)
 147 fid=fopen(out,'w');
 148 if fid < 0,
 149     error(['Cannot open output file ',out])
 150 end
 151 end