perimeter/perimeter_in.m

   1 function result=perimeter_in(long,lat,fire_area,wrfout,time_now,time,interval,count);
   2
   3 % Volodymyr Kondratenko           July 19 2013
   4
   5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   6
   7 % Input:
   8 %    long,lat      longtitude and latitude converted to meters
   9 %    fire_area     fire map,[0,1] array, where 0- not
  10 %                  burning area, >0 burning area, 1 - area that was burnt
  11 %    wrfout        name of the wrfout file, that is being used to read ros
  12 %    time_now      the perimeter time (s)
  13 %    time          update the wind every count time steps
  14 %    interval      time step in wrfout in seconds
  15 %    count         update the wind every count time steps
  16 %  JM NEVER CHANGE INPUT VARIABLES
  17 %
  18 % Output:
  19 %    Final Matrix of times of ignition will be printed to 'output_tign.txt' % JM use save tign instead, create tign.mat
  20 %
  21 %
  22 [n,m]=size(long);
  23
  24
  25 'perimeter_in'
  26
  27 % JM algorithm state is stored in arrays A D C
  28 %
  29 % A contains rows [i,j] of indices of nodes not burning that have at least one burning neighbor
  30 %   and time of ignition > time_now
  31
  32 % Reading Rate of spread
  33 ros=read_data_from_wrfout(wrfout,time); % JM should be read_ros_from_wrfout
  34
  35 % Getting matrix A from the initial tign_g, where
  36 A=get_perim_from_initial_tign(fire_area);
  37 % Computing 4d array of distances between a point and its 8 neighbors
  38 distance=get_distances(long,lat);
  39 % JM do not change this later, use adjustment array instead if you must
  40
  41 % Computing 4d array of differences of tign (delta_tign) between a point
  42 % and its 8 neighbors
  43 [delta_tign]=get_delta_tign(distance,ros);
  44
  45 % Initializing tign
  46 % Everything outside of the burning area is set to time_now,
  47 % inside is set to 0
  48 tign_in=zeros(n,m);
  49 tign_in=(fire_area(:,:)==0).*time_now;  %Should I use |a-b|<eps ? %JM never test reals on equality
  50
  51 % D - if D[i,j]=1, then the neighbors of the point [i,j]
  52 % will be updated only in the next timestep (only after we update ros)
  53 D=zeros(n,m);
  54
  55 changed=1;
  56 tic
  57     sprintf('size of A- %i',size(A,1))
  58 for istep=1:max(size(tign_in)),
  59     if size(A,1)==0,
  60                 'printed'
  61                 break
  62     end
  63
  64     time_toc=toc;
  65     str= sprintf('%f -- How long does it take to run step %i',time_toc,istep-1);
  66     tign_last=tign_in;
  67     contour(tign_in(2:end-1,2:end-1));drawnow
  68     [tign_in,distance,A,D]=tign_update(tign_in,A,D,delta_tign,time_now,distance,interval,ros);
  69 %%%    [tign_in,distance,A,D]=tign_update(tign_in,A,D,delta_tign,time_now,distance,interval,ros);
  70     [row,col]=find(D>0)
  71     contour(tign_in(2:end-1,2:end-1));title(sprintf('step %i',istep)),drawnow
  72
  73     changed=sum(tign_in(:)~=tign_last(:));
  74
  75     sprintf('%s \n step %i inside tign changed at %i points \n %f -- norm of the difference',str,istep,changed,norm(tign_in-tign_last))
  76     sprintf('size of A- %i',size(A,1))
  77 % % %     if (size(A,1)==0)
  78 % % %         if (time-count)>0
  79 % % %             'getting new ros'
  80 % % %             A=[];
  81 % % %             [A(:,1),A(:,2)]=find(D>0);
  82 % % %             sprintf('size of A- %i',size(A,1))
  83 % % %             D=zeros(n,m);
  84 % % %             time_now=time_now-count*interval
  85 % % %             time=time-1;
  86 % % %             ros=read_data_from_wrfout(wrfout,time);
  87 % % %             delta_tign=get_delta_tign(distance,ros);
  88 % % %
  89 % % %         elseif (time-count)<0
  90 % % %
  91 % % %         'time-count<0'
  92 % % %         end
  93 % % %     end
  94 % % %
  95 end
  96
  97 result=tign_in;
  98
  99 fid = fopen('output_tign.txt', 'w');
 100     dlmwrite('output_tign.txt', result, 'delimiter', '\t','precision', '%.4f');
 101     fclose(fid);
 102
 103 if changed~=0,
 104     'did not find fixed point inside'
 105    end
 106 end
 107
 108 %%% function [tign,distance,A,D]=tign_update(tign,A,D,delta_tign,time_now,distance,interval,ros)
 109 function [tign,distance,A,D]=tign_update(tign,A,D,delta_tign,time_now,distance,interval,ros)
 110
 111 % A - set of points whose neighbors will be updated
 112 % C - if C[i,j]=1, then the neighbors of the point [i,j]
 113 % will be updated in the next iteration within the same timestep
 114 % D - if D[i,j]=1, then the neighbors of the point [i,j]
 115 % will be updated only in the next timestep (only after we update ros)
 116
 117 C=zeros(size(tign));
 118
 119 for i=1:size(A,1)
 120     for dx=-1:1
 121         for dy=-1:1
 122                 if (tign(A(i,1)+dx,A(i,2)+dy)<time_now)==1
 123                     tign_new=tign(A(i,1),A(i,2))-0.5*(delta_tign(A(i,1)+dx,A(i,2)+dy,-dx+2,-dy+2)+delta_tign(A(i,1),A(i,2),2-dx,2-dy));
 124                     if (tign(A(i,1)+dx,A(i,2)+dy)<tign_new)&&(tign_new<=time_now)
 125                            tign(A(i,1)+dx,A(i,2)+dy)=tign_new;
 126                             C(A(i,1)+dx,A(i,2)+dy)=1;
 127
 128                             if (time_now-tign_new>interval)
 129                                 D(A(i,1),A(i,2))=1;
 130                             end
 131 %
 132 %                         if (time_now-tign_new<interval)
 133 %
 134 %                         else
 135 %                             D(A(i,1),A(i,2))=1;
 136 %                         end
 137
 138 % % %                         if (time_now-tign_new<interval)
 139 % % %                             tign(A(i,1)+dx,A(i,2)+dy)=tign_new;
 140 % % %                             C(A(i,1)+dx,A(i,2)+dy)=1;
 141 % % %                         else
 142 % % %                             distance(A(i,1),A(i,2),2-dx,2-dy)=distance(A(i,1),A(i,2),2-dx,2-dy)-interval*ros(A(i,1),A(i,2),2-dx,2-dy);
 143 % % %                             distance(A(i,1)+dx,A(i,2)+dy,2-dx,2-dy)=distance(A(i,1)+dx,A(i,2)+dy,2-dx,2-dy)-interval*ros(A(i,1)+dx,A(i,2)+dy,2-dx,2-dy);
 144 % % %                             D(A(i,1)+dx,A(i,2)+dy)=1;
 145 % % %                         end
 146
 147                     end
 148
 149             end
 150         end
 151     end
 152 end
 153 A=[];
 154 [A(:,1),A(:,2)]=find(C>0);
 155 end
 156
 157 function distance=get_distances(long,lat)
 158 % computing 4d array of distances between a point and its 8 neighbors
 159 %
 160 % input:
 161 %   long(i,j), lat(i,j) geographical coordinates of node [i,j], i=1:m, j=1:n, [m,n]=size(long)=size(lat)
 162 %
 163 % output
 164 %   distance(i+1,j+1,a+2,b+2) = geographical distance between node [i,j] and [i+a,j+b] , a,b=-1:1
 165
 166     distance=zeros(size(long,1),size(long,2),3,3);
 167
 168     for i=2:size(long,1)-1
 169       for j=2:size(long,2)-1
 170         for a=-1:1
 171             for b=-1:1
 172                 % distance between node [i,j] and [i+a,j+b]
 173                 distance(i,j,a+2,b+2)=sqrt((long(i+a,j+b,1)-long(i,j,1))^2+    ...
 174                           (lat(i+a,j+b,1)-lat(i,j,1))^2);
 175             end
 176         end
 177       end
 178     end
 179     % note that distance(i,i,...) makes no sense if i=2,m+1 or j=2,n+1 in the direction to the outside
 180
 181     % some ideas to look at
 182     % maybe use mirroring?
 183     % maybe keep numbering [i,j] always same as original outside of routines?
 184     % now need to remember which array is shifted and bordered by zeros and which one is not
 185     % if I was doing it: keep distance size (m,n,3,3) but not compute distance(1,1,...) etc
 186     % and keep it zero, same with delta_tign, etc.
 187
 188
 189 end
 190
 191 % maybe get_delta_tign ?
 192
 193 function delta_tign=get_delta_tign(distance,ros)
 194 % computing 4d array of differences of tign (delta_tign) between a point and its 8 neighbors
 195 % input:
 196 %    distance(i+1,j+1,a+2,b+2)  geographical distance between nodes [i,j] and [i+a,j+b]
 197 %                           from get_distances
 198 %    ros(i,j,a+2,b+2)   rate of spread at node [i,j] in the direction towards [i+a,j+b]
 199 %
 200 % output:
 201 %    delta_tign(i,j,a+2,b+2) time the fire takes to propagate from [i,j] to [i+a,j+b]
 202
 203 delta_tign=zeros(size(distance));
 204
 205 for i=2:size(delta_tign,1)-1
 206     for j=2:size(delta_tign,2)-1
 207         for a=-1:1
 208             for b=-1:1
 209                 delta_tign(i,j,a+2,b+2)=distance(i,j,a+2,b+2)/ros(i-1,j-1,a+2,b+2);
 210             end
 211         end
 212     end
 213 end
 214
 215 end
 216
 217 function A=get_perim_from_initial_tign(fire_area);
 218 % in:
 219 %    tign         ignition time
 220 % out:
 221 %    A            rows [i,j] of indices of nodes not burning that have at least one
 222 %                 burning neighbor
 223 A=[];
 224 format long
 225 for i=2:size(fire_area,1)-1
 226     for j=2:size(fire_area,2)-1
 227         % if (i,j) is not burning
 228         if (fire_area(i,j)==0)
 229             % if any neighbor is not burning
 230             if (any(any(fire_area(i-1:i+1,j-1:j+1)>0))==1)
 231             % add [i,j] to A
 232             A=[A;[i,j]];
 233             end
 234        end
 235     end
 236 end
 237 end
 238
 239
 240
 241 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 242 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 247 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 249 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 250 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 253 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 255 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 256
 257 % dead code
 258 function result=perimeter_in_2(long,lat,ros,time_now,bound,wrfout,interval,count)
 259
 260 % Volodymyr Kondratenko           December 8 2012
 261
 262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 263
 264 % The function creates the initial matrix of times of ignitions
 265 %                       given the perimeter and time of ignition at the points of the perimeter and
 266 %                       is using wind and terrain gradient in its calculations
 267 % Input: We get it after reading the data using function read_file_perimeter.m
 268 %
 269 %        long                   FXLONG*UNIT_FXLONG, longtitude coordinates of the mesh converted into meters
 270 %        lat                    FXLAT*UNIT_FXLAT, latitude coordinates of the mesh converted into meters
 271 %        uf,vf                  horizontal wind velocity vectors of the points of the mesh
 272 %        dzdxf,dzdyf    terrain gradient of the points of the mesh
 273 %        time_now               time of ignition on the fireline (fire perimeter)
 274 %        bound                  set of ordered points of the fire perimeter 1st=last
 275 %                                               bound(i,1)-horisontal; bound(i,1)-vertical coordinate
 276 %
 277 % Output:   Matrix of time of ignition
 278 %
 279 % Code:
 280
 281 %addpath('../../other/Matlab/util1_jan');
 282 %addpath('../../other/Matlab/netcdf');
 283 tic
 284 fuels % This function is needed to create fuel variable, that contains all the characteristics
 285       % types of fuel, this function lies in the same folder where you run the main_function.m
 286           % (where is it originally located/can be created?)
 287
 288 format long
 289
 290 bnd_size=size(bound);
 291 n=size(long,1);
 292 m=size(long,2);
 293
 294 %tign=zeros(n,m);      % "time of ignition matrix" of the nodes
 295 %A=zeros(n,m);         % flag matrix of the nodes, A(i,j)=1 if the time of ignition of the
 296                       % point (i,j) was updated at least once
 297
 298 'started'
 299 %  IN - matrix, that shows, whether the point is inside (IN(x,y)>0) the burning region
 300 %  or outside (IN(x,y)<0)
 301 %  ON - matrix that, shows whether the point is on the boundary or not
 302 %  Both matrices evaluated using "polygon", coefficients are multiplied by
 303 %  10^6, because the function looses acuracy when it deals with decimals
 304
 305 xv=bound(:,1);
 306 yv=bound(:,2);
 307 xv=xv*100000;
 308 yv=yv*100000;
 309 lat1=lat*100000;
 310 long1=long*100000;
 311 [IN,ON] = inpolygon(long1,lat1,xv,yv);
 312
 313 % Code
 314
 315 [delta_tign]=delta_tign_calculation(long,lat,ros);
 316
 317 % Calculates needed variables for rate of fire spread calculation
 318
 319 %%%%%%% First part %%%%%%%
 320 % Set everything inside to time_now and update the tign of the points outside
 321
 322 % Initializing flag matrix A and time of ignition (tign)
 323 % Extending the boundaries, in order to speed up the algorythm
 324 %A=[];
 325 %C=zeros(n+2,m+2);
 326 % A contains coordinates of the points that were updated during the last
 327 % step
 328
 329 IN_ext=(2)*ones(n+2,m+2);
 330 IN_ext(2:n+1,2:m+1)=IN(:,:,1);
 331
 332
 333 % Set all the points outside to time_now and update the points inside
 334
 335 % Initializing flag matrix A and time of ignition (tign)
 336 % Extending the boundaries, in order to speed up the algorythm
 337 A=[];
 338 C=zeros(n+2,m+2);
 339 for i=2:n+1
 340     for j=2:m+1
 341         if IN_ext(i,j)==0
 342             if sum(sum(IN_ext(i-1:i+1,j-1:j+1)))>0
 343             A=[A;[i,j]];
 344             end
 345        end
 346     end
 347 end
 348
 349 tign_in=zeros(n+2,m+2);
 350 tign_in(2:n+1,2:m+1)=(1-IN(:,:,1)).*time_now;
 351 changed=1;
 352
 353 time_old=time_now;
 354 % The algorithm stops when the matrix converges (tign_old-tign==0) or if the amount of iterations
 355 % % reaches the max(size()) of the mesh
 356 count
 357 interval
 358 count*interval
 359 for istep=1:max(size(tign_in)),
 360     if changed==0,
 361                 % The matrix of tign converged
 362                 'printed'
 363                 break
 364     end
 365
 366     tign_last=tign_in;
 367     time_toc=toc;
 368     str= sprintf('%f -- How long does it take to run step %i',time_toc,istep-1);
 369
 370
 371     if ((time_old-max(max(tign_in(A(:,1),A(:,2)))))>=(count*interval))&&((time-count)>0)
 372     'getting new ros'
 373       sprintf('time_old= %f',time_old)
 374       sprintf('tign_in(A(1,1),A(1,2))= %f',tign_in(A(1,1),A(1,2)))
 375       sprintf('min(min(tign_in(A(:,1),A(:,2))))= %f',min(min(tign_in(A(:,1),A(:,2)))))
 376       sprintf('max(max(tign_in(A(:,1),A(:,2))))= %f',max(max(tign_in(A(:,1),A(:,2)))))
 377       sprintf('time_old-max(max(tign_in(A(:,1),A(:,2)))= %f',time_old-max(max(tign_in(A(:,1),A(:,2)))))
 378
 379         count1=mod(time_old-max(max(tign_in(A(:,1),A(:,2)))),count*interval)
 380         time_old=time_old-count1*interval
 381         time=time-count1
 382        ros=read_data_from_wrfout(wrfout,size(long,1),size(long,2),time);
 383        delta_tign=delta_tign_calculation(long,lat,ros);
 384     end
 385
 386
 387     % tign_update - updates the tign of the points
 388     [tign_in,A,C]=tign_update(tign_in,A,IN_ext,delta_tign,time_now);
 389   % when it is outside the last parameter is 0, inside 1
 390     changed=sum(tign_in(:)~=tign_last(:));
 391 %    if (changed<=5)
 392 %       for i=1:size(A,1)
 393 %           A(i,:)
 394 %        end
 395 %    end
 396
 397     sprintf('%s \n step %i inside tign changed at %i points \n %f -- norm of the difference',str,istep,changed,norm(tign_in-tign_last))
 398    sprintf('size of A- %i',size(A,1))
 399 end
 400 final_tign=tign_in;
 401 %final_tign(2:n+1,2:m+1)=(IN(:,:,1)>0).*tign_in(2:n+1,2:m+1)+(IN(:,:,1)==0).*tign(2:n+1,2:m+1);
 402 result=final_tign(2:n+1,2:m+1);
 403
 404 fid = fopen('output_tign.txt', 'w');
 405     dlmwrite('output_tign.txt', result, 'delimiter', '\t','precision', '%.4f');
 406     fclose(fid);
 407
 408 if changed~=0,
 409     'did not find fixed point inside'
 410    end
 411 end
 412
 413
 414
 415
 416