perimeter_new/read_file_perimeter_jm.m

   1 function [fxlong,fxlat,fire_perimeter,timestep_end]=read_file_perimeter(wrfout,wrfout_fire,time,input_type,input_file)
   2 % Volodymyr Kondratenko           April 3 2012
   3
   4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   5
   6 % Input: wrfout : String with the name of the wrfout file,
   7 %                  It is needed for reading the latitude and longtitude
   8 %                  coordinates of the mesh and also UNIT_FXLONG and
   9 %                  UNIT_FXLAT variables
  10 %        time    : time step index in the wrfout (last index)
  11 %                  example FXLAT(:,:,time)
  12 % Output:
  13 %        p      structure with various fiels
  14 %        long = FXLONG, longtitude coordinates of the mesh converted to (m)
  15 %        lat = FXLAT, latitude coordinates of the mesh converted to (m)
  16 %        fire_area = fire map,[0,1] array, where 0- not
  17 %                  burning area, >0 burning area, 1 - area that was burnt
  18 %                  If input_type=1, then fire_area - set of ordered points
  19 %                  of the boundary 1st=last;
  20 %                  bound(i,1)-horisontal; bound(i,1)-vertical coordinate
  21
  22 datafile=sprintf('data_%s_%i_%i',wrfout,time,input_type);
  23 global saved_data  % 0 = read from original files and store in matlab files, 1=read saved data
  24 disp(['read_file_perimeter time=',num2str(time),' input_type=',num2str(input_type)]);
  25
  26 if (input_type==0)
  27     error('not supported')
  28     if saved_data
  29         w=load(datafile);
  30         p=w.p; q=w.q;
  31     else
  32         p=nc2struct(wrfout,{'UNIT_FXLONG','UNIT_FXLAT','FXLONG','FXLAT','Times','ITIMESTEP'},{'DT'},time);
  33         q=nc2struct(wrfout_fire,{'FIRE_AREA'},{},time);
  34         save(datafile,'p','q')
  35     end
  36     fire_area=q.fire_area;
  37     fxlong=p.fxlong*p.unit_fxlong;
  38     fxlat=p.fxlat*p.unit_fxlat;
  39 elseif (input_type==1)
  40     error('not supported')
  41     if saved_data
  42         w=load(datafile);
  43         p=w.p;
  44     else
  45         p=nc2struct(wrfout,{'UNIT_FXLONG','UNIT_FXLAT','FXLONG','FXLAT'},{},time);
  46         save(datafile,'p')
  47     end
  48     fxlong=p.fxlong*p.unit_fxlong;
  49     fxlat=p.fxlat*p.unit_fxlat;
  50     fire_area=read_perim_from_file(input_file,fxlong,fxlat,p.unit_fxlong,p.unit_fxlat);
  51
  52 % Test
  53     fid = fopen('output_fire_area.txt', 'w');
  54     dlmwrite('output_fire_area.txt', fire_area, 'delimiter', '\t','precision', '%.4f');
  55     fclose(fid);
  56 elseif (input_type==2)
  57     if saved_data
  58         disp(['loading from ',datafile])
  59         w=load(datafile);
  60         p=w.p;
  61     else
  62         p=nc2struct(wrfout,{'UNIT_FXLONG','UNIT_FXLAT','FXLONG','FXLAT','ITIMESTEP'},{'DT','DX','DY'},time);
  63         disp(['storing to ',datafile])
  64         save(datafile,'p')
  65     end
  66     disp(['reading fire_perimeter_big from ',input_file,' frame ',num2str(time)])
  67     fire_perimeter_big=dlmread(input_file);
  68     fxlong=p.fxlong*p.unit_fxlong;
  69     fxlat=p.fxlat*p.unit_fxlat;
  70     fire_area=zeros(size(fire_perimeter_big(:,1:2)));
  71     fire_perimeter(:,1)=fire_perimeter_big(:,1)*p.unit_fxlong;
  72     fire_perimeter(:,2)=fire_perimeter_big(:,2)*p.unit_fxlat;
  73     timestep_end=p.itimestep*p.dt;
  74     %fire_area=fire_area_big(:,1:2);
  75 end
  76
  77 end
  78
  79
  80 % Input                 data : String - data, that contains the name of the Text file.
  81 %                  First 2 columns - coordinates of all the
  82 %                  points on the boundary (lon,lat).
  83 %                  1rt row - time_now (second number is not needed, is set to 0);
  84 %                  2nd row - size of the mesh;
  85 %                  3rd row - coordinates of ignition point;
  86 %                  All next rows - coordinates of all the
  87 %                  points on the boundary (lon,lat). %
  88 % Output           bound - set of ordered points of the boundary 1st=last
  89 %                  bound(i,1)-horisontal; bound(i,1)-vertical coordinate
  90
  91
  92 function IN=read_perim_from_file(data,long,lat,unit_long,unit_lat)
  93 fid = fopen(data);
  94 bound = fscanf(fid,'%9g %*1s %8g %*3s',[2 inf]);
  95 % bound - set of ordered points of the boundary 1st=last
  96 %         bound(i,1)-horisontal; bound(i,1)-vertical coordinate
  97
  98 bound = bound';
  99 fclose(fid);
 100
 101 bound(:,1)=bound(:,1)*unit_long;
 102 bound(:,2)=bound(:,2)*unit_lat;
 103
 104 xv=bound(:,1);
 105 yv=bound(:,2);
 106 xv=xv*100000;
 107 yv=yv*100000;
 108 lat1=lat*100000;
 109 long1=long*100000;
 110 [IN,ON] = inpolygon(long1,lat1,xv,yv);
 111
 112 max(max(IN(:,:)))
 113
 114 end
 115