perimeter_new/read_file_perimeter.m

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