perimeter_new/perimeter_in.m

   1 function tign=perimeter_in(long,lat,fire_area,wrfout,time,interval,time_step,num_wrf, input_type)
   2
   3 % Volodymyr Kondratenko           July 19 2013
   4 % Ideas line 261
   5
   6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   7
   8 % Input:
   9 %    long,lat      longtitude and latitude converted to meters
  10 %    fire_area     fire map,[0,1] array, where 0- not
  11 %                  burning area, >0 burning area, 1 - area that was burnt
  12 %    wrfout        name of the wrfout file, that is being used to read ros
  13 %    time_now      the perimeter time (s)
  14 %    time          update the wind every count time steps
  15 %    interval      time step in wrfout in seconds
  16 %
  17 % Output:
  18 %    Final Matrix of times of ignition will be printed to 'output_tign.txt' % JM use save tign instead, create tign.mat
  19
  20 % Algorithm state is stored in arrays A D C
  21 %
  22 % A contains rows [i,j] of indices of nodes not burning that have at least one burning neighbor
  23 %   and time of ignition > time_now
  24 % Computing 4d array of distances between a point and its 8 neighbors
  25 %ideas
  26 % Change to this later (if it speeds up the algorithm)
  27                % [row,col]=find(C(A(jj,1)-1:A(jj,1)+1,A(jj,2)-1:A(jj,2)+1)==0);
  28                % for index=1:size(row,1)
  29                %   dx=row(index)-2;
  30                %   dy=col(index)-2;
  31 % Iine 180, maybe I hould do it each step toavoid extra calculations
  32 data_steps='started perimeter_in';
  33
  34 distance=get_distances(long,lat);
  35
  36 %case when the fire_area is given in a file
  37 if (input_type==2)
  38  xv=fire_area(:,1);
  39  yv=fire_area(:,2);
  40  xv=xv*100000;
  41  yv=yv*100000;
  42  lat1=lat*100000;
  43  long1=long*100000;
  44  [IN,ON] = inpolygon(long1,lat1,xv,yv);
  45  fire_area = IN(:,:,1);
  46 end
  47
  48
  49 tign=get_tign_from_dif_eq(wrfout,fire_area,distance,time,interval,time_step,num_wrf,data_steps);
  50
  51 fid = fopen('output_tign_test.txt', 'w');
  52     dlmwrite('output_tign_test.txt', tign, 'delimiter', '\t','precision', '%.4f');
  53     fclose(fid);
  54
  55 end
  56
  57
  58 function [tign]=get_tign_from_dif_eq(wrfout,fire_area,distance,time,interval,time_step,num_wrf,data_steps)
  59
  60 % Getting matrix A from the initial tign_g, where
  61 [A,C,D,tign]=get_perim_from_initial_tign(fire_area,time,interval,time_step,num_wrf);
  62 %contour(C);title(sprintf('Original perimeter')); drawnow
  63 clear fire_area
  64
  65 I=zeros(size(distance));            % Matrix of distances
  66 P=A; %Points on the perimeter, who has at least one neighbor, that was not updated yet.
  67 C_old=C;
  68 ros_old=read_ros_from_wrfout(wrfout{num_wrf},time);
  69 for ts=(time_step*(num_wrf-1)+time):-1:2 % ts -time step
  70    if (ts<time) % At the first step we initialize A and D from get_perim_from_initial_tign
  71       A=[];
  72       [A(:,1),A(:,2)]=find(C(2:end-1,2:end-1)==1);
  73       A(:,1)=A(:,1)+1;
  74       A(:,2)=A(:,2)+1;
  75       for k=1:size(A,1)
  76          D(A(k,1),A(k,2))=interval;
  77       end
  78       if ~isempty(P)
  79           P_new=[];
  80         for l=1:size(P,1)
  81            if (any(any(C(P(l,1)-1:P(l,1)+1,P(l,2)-1:P(l,2)+1)==0))==1)
  82              P_new=[P_new; P(l,:)];
  83              D(P(l,1),P(l,2))=interval;
  84            end
  85         end
  86         P=P_new;
  87         A=[A;P];
  88       end
  89
  90    end
  91  %  if size(A,1)==0
  92  %     data_steps=sprintf('%s\n finished at the time %i',data_steps,ts);
  93  %     fid = fopen(myfile, 'w');
  94  %     fprintf(fid,'%s',data_steps);
  95  %     fclose(fid);
  96  %     break;
  97  %  end
  98   % data_steps=sprintf('%s \n Step %i',data_steps,ts);
  99   % data_steps=sprintf('%s \n first element of A',data_steps);
 100   % data_steps=sprintf('%s \n %s',data_steps,mat2str(A(1,:)));
 101   % data_steps=sprintf('%s \n C around [%i %i]',data_steps,pnt_a,pnt_b);
 102   % data_steps=print_big_mat(data_steps,pnt_a,pnt_b,C);
 103   % data_steps=sprintf('%s \n tign around [%i %i]',data_steps,pnt_a,pnt_b);
 104   % data_steps=print_big_mat(data_steps,pnt_a,pnt_b,tign);
 105
 106    cur_num_wrf=ceil((ts-1)/time_step);
 107    cur_ts=mod((ts-1),time_step);
 108    if (cur_ts==0)
 109        cur_ts=time_step;
 110    end
 111    ros_new=read_ros_from_wrfout(wrfout{cur_num_wrf},cur_ts);
 112    tign_old=tign;
 113    [tign,C,D,I,data_steps]=get_tign_one_timestep(tign,ros_old,ros_new,A,C,D,I,distance,interval,ts,data_steps);
 114    t=tign;t(t==0)=NaN;figure(1); mesh(t);title(sprintf('step %i, tign',ts)); colorbar; drawnow
 115    figure(2);tt=tign;tt(tign==tign_old)=NaN;mesh(tt);title(sprintf('step %i, tign diff',ts)); colorbar; drawnow
 116
 117    data_steps=sprintf('%s \n Main cycle for step %i is over',data_steps,ts);
 118    % figure(2); contour(C);title(sprintf('step %i, Matrix C, before subfunction',ts)); drawnow
 119    if any(any(D~=0))
 120       data_steps=sprintf('%s \n Error: D needs to be=0',data_steps);
 121    end
 122    ros_old=ros_new;
 123    %figure(1); contour(tign);title(sprintf('step %i, tign',ts)); drawnow
 124    %figure(2); contour(C);title(sprintf('step %i, Matrix C',ts)); drawnow
 125    changed=sum(C(:)~=C_old(:));
 126    data_steps=sprintf('%s\n After a cycle tign changed in %i points',data_steps,changed);
 127 %   fid = fopen(myfile, 'w');
 128 %   fprintf(fid,'%s',data_steps);
 129 %   fclose(fid);
 130    if changed>0,
 131         C_old=C;
 132    end
 133 end
 134 index=find(C==0);
 135 data_steps=sprintf('%s\n All the steps are done, but there are still %i points, that were not updated yet',data_steps,size(index));
 136 %fid = fopen(myfile, 'w');
 137 %fprintf(fid,'%s',data_steps);
 138 %fclose(fid);
 139 end
 140
 141
 142 function [tign,C,D,I,data_steps]=get_tign_one_timestep(tign,ros_old,ros_new,B,C,D,I,distance,interval,ts,data_steps)
 143
 144 step=1;
 145 data_steps=sprintf('%s\n subcycle in step %i',data_steps,ts);
 146 while any(any(D>0))
 147    tign_old=tign;
 148    data_steps=sprintf('%s\n substep %i',data_steps,step);
 149    data_steps=sprintf('%s\n points whose neighbors are being updated (size(B,1)=) %i',data_steps,size(B,1));
 150    D_old=D;
 151    for j=1:size(B,1)
 152       for dx=-1:1
 153          for dy=-1:1
 154             if (C(B(j,1)+dx,B(j,2)+dy)==0)
 155                F=0.25*D(B(j,1),B(j,2))*(ros_old(B(j,1),B(j,2),2-dx,2-dy)+ros_new(B(j,1),B(j,2),2-dx,2-dy) + ...
 156                  ros_old(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)+ros_new(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy));
 157              I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)=(step==1)*I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy);
 158              if (I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)>0)&&(step>1)
 159                   data_steps=sprintf('%s\n Error(5): I happens to be more than 0',data_steps);
 160              end
 161                if (I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)+F>=distance(B(j,1),B(j,2),2-dx,2-dy))
 162                   tign(B(j,1)+dx,B(j,2)+dy)=min(tign(B(j,1),B(j,2)),ts*interval)- ...
 163                     ((distance(B(j,1),B(j,2),2-dx,2-dy)-I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy))/F)*(D(B(j,1),B(j,2)));
 164                   C(B(j,1)+dx,B(j,2)+dy)=1;
 165                   D(B(j,1)+dx,B(j,2)+dy)=tign(B(j,1)+dx,B(j,2)+dy)-(ts-1)*interval;
 166                   if (D(B(j,1)+dx,B(j,2)+dy)<0)
 167                      data_steps=sprintf('%s\n Error(3) D happens to be less than 0',data_steps);
 168                      data_steps=sprintf('%s\n At point (i,j)%d%d',data_steps,i,j);
 169                      data_steps=sprintf('%s\n B(j,1),B(j,2), B(j,1)+dx,B(j,2)+dy %d%d%d%d',data_steps,B(j,1),B(j,2), B(j,1)+dx,B(j,2)+dy);
 170                      data_steps=sprintf('%s\n D(B(j,1)+dx,B(j,2)+dy %d)',data_steps,D(B(j,1)+dx,B(j,2)+dy));
 171                   end
 172                else
 173                   I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)=I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)+F;
 174                end
 175             elseif (C(B(j,1)+dx,B(j,2)+dy)==1)
 176               F=0.25*D(B(j,1),B(j,2))*(ros_old(B(j,1),B(j,2),2-dx,2-dy)+ros_new(B(j,1),B(j,2),2-dx,2-dy) + ...
 177                  ros_old(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)+ros_new(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy));
 178                if (I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)+F>=distance(B(j,1),B(j,2),2-dx,2-dy))
 179                   tign_new=min(tign(B(j,1),B(j,2)),ts*interval)- ...
 180                     ((distance(B(j,1),B(j,2),2-dx,2-dy)-I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy))/F)*(D(B(j,1),B(j,2)));
 181                   I(B(j,1)+dx,B(j,2)+dy,2-dx,2-dy)=0;
 182                   if tign_new>tign(B(j,1)+dx,B(j,2)+dy)
 183                      tign(B(j,1)+dx,B(j,2)+dy)=tign_new;
 184                      D(B(j,1)+dx,B(j,2)+dy)=tign(B(j,1)+dx,B(j,2)+dy)-(ts-1)*interval;
 185                      if (D(B(j,1)+dx,B(j,2)+dy)<0)
 186                         data_steps=sprintf('%s\n Error(4): D happens to be less than 0',data_steps);
 187                      end
 188                   end
 189                end
 190             end
 191          end
 192       end
 193       D(B(j,1),B(j,2))=0;
 194    end
 195    % figure(3); contour(C);title(sprintf('step %i, Matrix C in subfunction substep %i',ts2,step)); drawnow
 196    step=step+1;
 197    [row,col]=find((D_old==D)&(D_old~=0));
 198    if (size(row,1)~=0)
 199       data_steps=sprintf('%s\n Error find((D_old==D_new)&(D_old~=0)) is not equal, D(row(1),col(1))= %f',data_steps,D(row(1),col(1)));
 200    end
 201    % I do this part in the end, since I initialize A for the first step
 202    B=[];
 203    [B(:,1),B(:,2)]=find(D(2:end-1,2:end-1)>0);
 204    B(:,1)=B(:,1)+1;
 205    B(:,2)=B(:,2)+1;
 206    if any(tign(:)~=tign_old(:)),
 207         figure(3);tt=tign;tt(tign==tign_old)=NaN;mesh(tt);title(sprintf('step %i, tign new only')); colorbar; drawnow
 208    end
 209 end
 210 [row,col]=find(C==1);
 211 for b=1:size(row,1)
 212    if ~any(any(C(row(b)-1:row(b)+1,col(b)-1:col(b)+1)==0))
 213       C(row(b),col(b))=2;
 214    end
 215 end
 216 end
 217
 218 function distance=get_distances(long,lat)
 219 % computing 4d array of distances between a point and its 8 neighbors
 220 %
 221 % input:
 222 %   long(i,j), lat(i,j) geographical coordinates of node [i,j], i=1:m, j=1:n, [m,n]=size(long)=size(lat)
 223 %
 224 % output
 225 %   distance(i+1,j+1,a+2,b+2) = geographical distance between node [i,j] and [i+a,j+b] , a,b=-1:1
 226
 227 distance=zeros(size(long,1),size(long,2),3,3);
 228 i=2:size(long,1)-1;
 229 j=2:size(long,2)-1;
 230 for a=-1:1
 231    for b=-1:1
 232       %for j=2:size(long,2)-1
 233       %   for i=2:size(long,1)-1
 234       %      % distance between node [i,j] and [i+a,j+b]
 235             distance(i,j,a+2,b+2)=sqrt((long(i+a,j+b,1)-long(i,j,1)).^2+(lat(i+a,j+b,1)-lat(i,j,1)).^2);
 236       %   end
 237       %end
 238    end
 239 end
 240 %IDEAS
 241 % some ideas to look at
 242 % maybe use mirroring?
 243 % maybe keep numbering [i,j] always same as original outside of routines?
 244 % now need to remember which array is shifted and bordered by zeros and which one is not
 245 % if I was doing it: keep distance size (m,n,3,3) but not compute distance(1,1,...) etc
 246 % and keep it zero, same with delta_tign, etc.
 247 end
 248
 249 function [A,C,D,tign]=get_perim_from_initial_tign(fire_area,time,interval,time_step,num_wrf)
 250
 251 % in:
 252 % tign    ignition time
 253 % out:
 254 % C = 3 - area outside of the fire perimeter;
 255 %     2 - area whose tign was already computed and it was used to
 256 %         compute the tign of its neighbors + not burning points that lie
 257 %         on the perimeter;
 258 %     1 - area whose tign was computed, but may still updated from
 259 %         neighbors. This area will be used to compute the neighbors either
 260 %         next time step or during the subcycle;
 261 %     0 - area inside the fire perimeter, that needs to be computed;
 262 % A       rows [i,j] of indices of nodes not burning that have at least one
 263 %         burning neighbor. (equal to area where to C=2, but only in the
 264 %         first time step)
 265 C=4*(fire_area==0);
 266 D=zeros(size(fire_area));
 267 tign=zeros(size(fire_area));
 268 format long
 269 for i=2:size(fire_area,1)-1
 270    for j=2:size(fire_area,2)-1
 271       % if (i,j) is not burning
 272       if (fire_area(i,j)==0)
 273          % if any neighbor is burning
 274          if (any(any(fire_area(i-1:i+1,j-1:j+1)>0))==1)
 275             % add [i,j] to A
 276             C(i,j)=3;
 277             D(i,j)=interval;
 278           % why?  tign(i,j)=interval*(time_step*(num_wrf-1)+time-1); % I do that because time starts with 00:00
 279           % Why did I have above line before
 280             tign(i,j)=interval*(time_step*(num_wrf-1)+time); % I do that because time starts with 00:00
 281          end
 282       end
 283    end
 284 end
 285 [row,col]=find(C==3);
 286 A=[row,col];
 287 end
 288