ignition/ignition7.m

   1 function A=ignition7(data,wrf_out)
   2
   3 % The function Ignition plots and showes how the fire was propagating,
   4 %           given ignition point, boundary region and time_now
   5 % Input:    String - data, that contains the name of the Text file.
   6 %           First 2 columns - coordinates of all the
   7 %           points on the boundary (lon,lat).
   8 %           1rt row - time_now (second number is not needed, is set to 0);
   9 %           2nd row - size of the mesh;
  10 %           3rd row - coordinates of ignition point;
  11 %           All next rows - coordinates of all the
  12 %           points on the boundary (lon,lat).
  13 %
  14 % Output:   Matrix of time of ignition
  15 %
  16 %
  17 % Volodymyr Kondratenko           February 17 2011
  18 %--------------------------------------
  19 % Command line: what to do to call the function
  20 %--------------------------------------
  21 %
  22 % This code is done for the real case for the witch_fire
  23 % put the file and data.txt in the run folder
  24 %
  25 %  addpath ../../other/Matlab/util1_jan
  26 %  addpath ../../other/Matlab/netcdf
  27   addpath('../../other/Matlab/util1_jan');
  28   addpath('../../other/Matlab/netcdf');
  29   %  data='data.txt';
  30 %  wrf_out='wrfout_d01_2007-10-21_12:00:00_real'; // or if you use another wrf_out, that you have
  31 %  B=ignition6(data,wrf_out);
  32
  33
  34 % Code:
  35 % Getting Latitude
  36 % 1st row - # of latitude rows
  37 % 2nd row - # of longtitude
  38 % Rest coordinates of the grid points
  39 format long
  40 %var=ncload(wrf_out);
  41 unit_long=ncread(wrf_out,'UNIT_FXLONG');
  42 unit_lat=ncread(wrf_out,'UNIT_FXLAT');
  43 unit_long=unit_long(1);
  44 unit_lat=unit_lat(1);
  45 long=ncread(wrf_out,'FXLONG');
  46 lat=ncread(wrf_out,'FXLAT');
  47 long=long*unit_long;
  48 lat=lat*unit_lat;
  49 mat_size=size(long);
  50 grid_1=mat_size(1);
  51 grid_2=mat_size(2);
  52
  53
  54 fid = fopen(data);
  55 data = fscanf(fid,'%g %g',[2 inf]); % It has two rows now.
  56 data = data';
  57 fclose(fid)
  58 data_size=size(data);
  59
  60 time_now=data(1,1);  % time of ignition on the boundary
  61 mesh_size=data(2,:); % size of the matrix
  62 ign_pnt=data(3,:);   % coordinates of the ignition point
  63 ign_pnt(1)=ign_pnt(1)*unit_long;
  64 ign_pnt(2)=ign_pnt(2)*unit_lat;
  65 bound=data(4:data_size(1),:);
  66 bound(:,1)=bound(:,1)*unit_long;
  67 bound(:,2)=bound(:,2)*unit_lat;
  68
  69 % bound - set of ordered points of the boundary 1st=last
  70 % bound(i,1)-horisontal; bound(i,1)-vertical coordinate
  71 bnd_size=size(bound);
  72
  73 B=zeros(mesh_size);        % "time of ignition matrix" of the nodes
  74                            % Ignition point corresponds to B()=0
  75 C=zeros(mesh_size);        % Matrix of distances from ignition point to all points
  76
  77
  78 %  IN - matrix that shows, whether the point is inside (IN(x,y)>0) the burning region
  79 %  or outside (IN(x,y)<0)
  80 %  ON - matrix that, shows whether the point is on the boundary or not
  81 %  Both matrices evaluated using "polygon"
  82 xv=bound(:,1);
  83 yv=bound(:,2);
  84 xv=xv*100000;
  85 yv=yv*100000;
  86 lat1=lat*100000;
  87 long1=long*100000;
  88 [IN,ON] = inpolygon(long1,lat1,xv,yv);
  89
  90 %xv(5)
  91 %yv(5)
  92 %lat1(5,5)
  93 %long1(5,5)
  94
  95 % Calculation of the matrix of distances
  96 % || ( (coord_xf(i,j) - x)*unit_xf , (coord_yf(i,j) - y)*unit_yf ) ||
  97 % unit_fxlat=pi2/(360.*reradius)
  98 % unit_fxlong=cos(lat_ctr*pi2/360.)*unit_fxlat
  99 % lat_ctr=config_flags%cen_lat
 100 %reradius,    & ! 1/earth radiusw
 101 %                 pi2   ! 2*pi
 102
 103 %for i=1:grid_1
 104 %     for j=1:grid_2
 105 %          C(i,j)=sqrt((unit_long(1,1,1)*(long(i)-ign_pnt(1)))^2+(unit_lat(1,1,1)*(lat(i)-ign_pnt(2)))^2);
 106 %     end
 107 %end
 108
 109 %--------------------------------------------------------------------%
 110 % Algorythm
 111 %--------------------------------------------------------------------%
 112
 113 % For each point in the burning area, we
 114 %     a) Build a line through this point and ignition point
 115 %     b) Find two(or one) points between which, line intersects boundary area,
 116 %        should be 2 variants
 117 %     This is done by checking the F(x,y)=(y-y1)*(x2-x1)-(x-x1)*(y2-y1)
 118 %     Needed points are that satisfy 2 conditions:
 119 %      1) they are sequential;
 120 %      2) F(x1,y1)*F(x2,y2)<0
 121 %     (In ase with one point, F(x,y)=0)
 122 %     c) Find out which side is right (from 2 possible choices of points)
 123 %        By similiar way
 124 %     d) Compute tign of the point by using Thales's theorem
 125 %     (it is not the same as on Wikipedia, see the link below, page 4)
 126 %        http://legacy.lclark.edu/~istavrov/geo-minithales-07.pdf
 127 %
 128
 129 aa=0;
 130 eps=0.1;
 131 for i=1:grid_1
 132     for j=1:grid_2
 133         if (ON(i,j)>0)          % The point we check is on the boundary
 134             B(i,j)=time_now;
 135             k=-1;
 136         elseif IN(i,j)>0        % The point is inside the burning region
 137             a_old=line_sign(ign_pnt(1),ign_pnt(2),long(i,j,1),lat(i,j,1),bound(1,1),bound(1,2));
 138             k=2;
 139             while (k>0)&&(k<=bnd_size(1))
 140                 if (a_old==0)
 141                     a_new=a_old;
 142                     k=1;
 143                 else
 144                     a_new=line_sign(ign_pnt(1),ign_pnt(2),long(i,j,1),lat(i,j,1),bound(k,1),bound(k,2));
 145                 end
 146                 if a_old*a_new<0
 147                     % Check if the point is on the line between
 148                     % the ignition point and 2 boundary points
 149                     a1=line_sign(long(i,j,1),lat(i,j,1),bound(k,1),bound(k,2),ign_pnt(1),ign_pnt(2));
 150                     a2=line_sign(long(i,j,1),lat(i,j,1),bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2));
 151                     if a1*a2<0
 152
 153                         dist1=line_dist(bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2),ign_pnt(1),ign_pnt(2));
 154                         dist2=line_dist(bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2),long(i,j,1),lat(i,j,1));
 155                         B(i,j)=time_now*(dist1-dist2)/dist1;
 156                         k=-1;
 157                     end
 158                 elseif a_new==0
 159                     % Case if the line goes exactly through the boundary point
 160                     % Check if the point lies between ignition point and boundary point
 161                     b1=sqrt((long(i,j,1)-ign_pnt(1))^2+(lat(i,j,1)-ign_pnt(2))^2);
 162                     b2=sqrt((long(i,j,1)-bound(k,1))^2+(lat(i,j,1)-bound(k,2))^2);
 163                     b3=sqrt((bound(k,1)-ign_pnt(1))^2+(bound(k,2)-ign_pnt(2))^2);
 164                     if (b1+b2<b3+eps)&&(b1+b2>b3-eps)
 165                         B(i,j)=time_now*b1/b3;
 166                         k=-1;
 167                     else
 168                         a_new=line_sign(ign_pnt(1),ign_pnt(2),long(i,j,1),lat(i,j,1),bound(k+1,1),bound(k+1,2));
 169                         k=k+1;
 170                     end
 171                 end
 172                 a_old=a_new;
 173                 k=k+1;
 174             end
 175         else
 176
 177            % point is outside the burning region
 178            % B(i,j)=time_now+1;
 179         a_old=line_sign(ign_pnt(1),ign_pnt(2),long(i,j,1),lat(i,j,1),bound(1,1),bound(1,2));
 180             k=2;
 181             while (k>0)&&(k<=bnd_size(1))
 182                 if (a_old==0)
 183                     a_new=a_old;
 184                     k=1;
 185                 else
 186                     a_new=line_sign(ign_pnt(1),ign_pnt(2),long(i,j,1),lat(i,j,1),bound(k,1),bound(k,2));
 187                 end
 188                 if a_old*a_new<0
 189                     % The point and second boundary point should be on the
 190                     % same side from the line going between ignition point
 191                     % and first boundary point
 192                     a1=line_sign(bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2),ign_pnt(1),ign_pnt(2));
 193                     a2=line_sign(bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2),long(i,j,1),lat(i,j,1));
 194                     if a1*a2<0
 195                         dist1=line_dist(bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2),ign_pnt(1),ign_pnt(2));
 196                         b1=sqrt((long(i,j,1)-ign_pnt(1))^2+(lat(i,j,1)-ign_pnt(2))^2);
 197                         %   dist2=line_dist(bound(k,1),bound(k,2),bound(k-1,1),bound(k-1,2),long(i,j,1),lat(i,j,1));
 198                         B(i,j)=time_now*b1/dist1;
 199
 200                         k=-1;
 201                     end
 202                 elseif a_new==0
 203                     % Case if the line goes exactly through the boundary point
 204                     % Check if the boundary point lies between ignition point and the point
 205                     b1=sqrt((long(i,j,1)-ign_pnt(1))^2+(lat(i,j,1)-ign_pnt(2))^2);
 206                     b2=sqrt((long(i,j,1)-bound(k,1))^2+(lat(i,j,1)-bound(k,2))^2);
 207                     b3=sqrt((bound(k,1)-ign_pnt(1))^2+(bound(k,2)-ign_pnt(2))^2);
 208                     if (b2+b3<b1+eps)&&(b2+b3>b1-eps)
 209                         B(i,j)=time_now*b1/b3;
 210                         k=-1;
 211                     else
 212                         a_new=line_sign(ign_pnt(1),ign_pnt(2),long(i,j,1),lat(i,j,1),bound(k+1,1),bound(k+1,2));
 213                         k=k+1;
 214                     end
 215                 end
 216                 a_old=a_new;
 217                 k=k+1;
 218           end
 219
 220         end
 221     end
 222 end
 223
 224 % check if ign_pnt(1), ign_pnt(2) - integers,
 225 % than set B(ignition point)=0
 226 if (rem(ign_pnt(1),1)==0) && (rem(ign_pnt(2),1)==0)
 227     B(ign_pnt(1),ign_pnt(1))=0;
 228 end
 229
 230 % Writing the data to the file data_out.txt
 231 fid = fopen('data_out.txt', 'w');
 232 dlmwrite('data_out.txt', B, 'delimiter', '\t','precision', '%.4f');
 233 fclose(fid);
 234
 235 write_array_2d('data_out1.txt',B)
 236
 237 A=B;
 238
 239 % Plot the results
 240 % questio: if ignition point coordinates are real and not on the mesh, than
 241 % plot will not print them
 242 %size(B)
 243 %grid_1
 244 %grid_2
 245 %x=1:1:grid_2;
 246 %y=1:1:grid_1;
 247 %figure(3)
 248 %surf(x,y,B)
 249 %B
 250