chem/wrfout2bluesky.m

   1 function wrfout2bluesky(in,out)
   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
   8 % Jan Mandel, July 2013
   9
  10 if ~exist('in','var')
  11     in=input('enter wrfout file name: ','s');
  12 end
  13 if ~iscell(in),
  14     in = {in};
  15 end
  16 if ~exist('out','var')
  17     out=input('enter bluesky csv file name: ','s');
  18     if isempty(out),
  19         out='bluesky.csv'; % default
  20     end
  21 end
  22 fprintf('reading file %s\n',in{1});
  23 fprintf('writing file %s\n',out);
  24 w=nc2struct(in{1},{'XLONG','XLAT','FXLONG','FXLAT'},{'DX','DY'},1); % get what does not change from the first netcdf file
  25 dx=w.dx;
  26 dy=w.dy;
  27 [m,n]=size(w.xlong);
  28 [fm,fn]=size(w.fxlong);
  29 srx=fm/m;
  30 sry=fn/n;
  31 fdx=dx/srx;
  32 fdy=dy/sry;
  33 fprintf('mesh size %i by %i\n',m,n)
  34 fprintf('fire mesh size %i by %i\n',fm,fn)
  35 fprintf('refinement ratio %g %g\n',srx,sry)
  36 fprintf('fire meshstep %g %g\n',fdx,fdy);
  37 lon=w.fxlong(:,:,1);lon=lon(:);
  38 lat=w.fxlat(:,:,1);lat=lat(:);
  39 min_x=min(lon);
  40 max_x=max(lon);
  41 min_y=min(lat);
  42 max_y=max(lat);
  43 lon_ctr=mean(lon);
  44 lat_ctr=mean(lat);
  45 fprintf('coordinates %g to %g by %g to %g\n',min_x,max_x,min_y,max_y);
  46 fprintf('the center of the domain is at coordinates %g %g\n',lon_ctr,lat_ctr)
  47 fid=fopen(out,'w');
  48 if fid < 0,
  49     error(['Cannot open output file ',out])
  50 end
  51 lat=w.fxlat(:,:,1);
  52 lon=w.fxlong(:,:,1);
  53 w=nc2struct(in{1},{'FIRE_AREA'},{},1); % the first fire area
  54 a_old=w.fire_area;
  55 fire_id=0;
  56 for ifile=1:length(in),
  57     fprintf('reading file %s\n',in{ifile});
  58     w=nc2struct(in{ifile},{'Times'},{}); % read netcdf file
  59     t=char(w.times');  % convert times to character matrix one string per row
  60     s=size(t,1);
  61     fprintf('timesteps %i meshstep %g %g\n',s,dx,dy);
  62     for istep=1:s
  63         w=nc2struct(in{ifile},{'FIRE_AREA'},{},istep); % read step istep
  64         a=w.fire_area;
  65         ta=sum(a(:))*fdx*fdy;
  66         fprintf('timestep %i/%i %s total fire area %g m^2\n',istep,s,t(istep,:),ta)
  67         d=datenum(t(istep,:),'yyyy-mm-dd_HH:MM:SS'); % convert ESMF time string to double
  68         tim=[datestr(d,'yyyymmddHHMM'),'Z'];         % convert to bluesky time format
  69         v = datevec(d);
  70         if(v(5)==0 & v(6) == 0) % whole hour
  71             a_diff = a - a_old;
  72             if any(a_diff(:)< - eps(single(1)))
  73                 warning('fire area can only increase')
  74                 a_diff = max(a_diff,0);
  75             end
  76             if(ta>0),
  77                 cc=bwconncomp(a>0); % find connected components,  image processing toolbox
  78                 for id=1:length(cc.PixelIdxList)
  79                     fire_id = fire_id+1;
  80                     fire_ctr_lon=mean(lon(cc.PixelIdxList{id}));  % center longitude of the fire area
  81                     fire_ctr_lat=mean(lat(cc.PixelIdxList{id}));  % center latitude of the fire area
  82                     fire_acres=fdx*fdy*sum(a_diff(cc.PixelIdxList{id}))/4046.86; % newly burning acres
  83                     fprintf(fid,'%i, %g, %g, %g, %s\n',fire_id,fire_ctr_lon,fire_ctr_lat,fire_acres,tim);
  84                     fprintf('%i, %g, %g, %g, %s\n',fire_id,fire_ctr_lon,fire_ctr_lat,fire_acres,tim);
  85                 end
  86             end
  87             a_old=a;
  88         else
  89             fprintf('not on the hour, skipping %s\n',tim)
  90         end
  91     end
  92 end
  93 fclose(fid);
  94 end
  95