perimeter_new/perimeter_in_jm.m

   1 function tign=perimeter_in_jm(long,lat,fire_area,wrfout,time,interval,time_step,num_wrf, input_type)
   2
   3 % Volodymyr Kondratenko           July 19 2013
   4 % Rebuilt Jan Mandel May 2016
   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 tign=get_tign_from_dif_eq(wrfout,fire_area,distance,time,interval,time_step,num_wrf,long,lat);
  49
  50 fid = fopen('output_tign_test.txt', 'w');
  51     dlmwrite('output_tign_test.txt', tign, 'delimiter', '\t','precision', '%.4f');
  52     fclose(fid);
  53
  54 end
  55
  56
  57 function [tign]=get_tign_from_dif_eq(wrfout,fire_area,distance,time,interval,time_step,num_wrf,long,lat)
  58
  59 time_now=interval*(time_step*(num_wrf-1)+time) % because time starts with 00:00
  60 time_perimeter=time_now;
  61 time_end=time_now;  % how long will propagate
  62
  63 [tign,fire_mask_out,fire_mask_in]=initial_tign(fire_area,time_now,time_end);
  64 perimeter_mask=fire_mask_in & fire_mask_out;
  65 ros_old=read_ros_from_wrfout(wrfout{num_wrf},time);
  66 figure(2);plot_ros(long,lat,ros_old);
  67 time_max=max(tign(:));
  68 fprintf('propagating out from perimeter time %g to %g\n',time_now,time_max)
  69 [t,d]=propagate_init(tign,distance);
  70 [t,d]=propagate(t,d,1,~fire_area,fire_mask_out,distance,ros_old,time_max,1);
  71 tign=t(:,:,2,2);
  72 err=big(tign(perimeter_mask)-time_now); fprintf('tign change on the perimeter %g\n',err)
  73 figure(1);mesh(long,lat,tign);drawnow
  74 [t,d]=propagate_init(tign,distance);
  75
  76 for ts=(time_step*(num_wrf-1)+time):-1:2 % ts -time step
  77     cur_num_wrf=ceil((ts-1)/time_step);
  78     cur_ts=mod((ts-1),time_step);
  79     if (cur_ts==0)
  80         cur_ts=time_step;
  81     end
  82     ros_new=read_ros_from_wrfout(wrfout{cur_num_wrf},cur_ts);
  83     cur_time_beg=interval*(time_step*(cur_num_wrf-1)+cur_ts);
  84     fprintf('propagating back in time to %g\n',cur_time_beg)
  85     [t,d]=propagate(t,d,-1,fire_area,fire_mask_in,distance,ros_new,cur_time_beg,1);
  86     tign=t(:,:,2,2);
  87     err=big(tign(perimeter_mask)-time_now); fprintf('tign change on the perimeter %g\n',err)
  88     tt=tign;
  89     tt( tt > time_perimeter | tt==0)=NaN;
  90     figure(3);mesh(long,lat,tt);title(num2str(cur_time_beg));
  91     hold on; contour3(long,lat,tt,[time_perimeter,time_perimeter]),hold off
  92     figure(1);plot_ros(long,lat,ros_new);
  93     drawnow
  94 end
  95 end
  96
  97 function distance=get_distances(long,lat)
  98 % computing 4d array of distances between a point and its 8 neighbors
  99 %
 100 % input:
 101 %   long(i,j), lat(i,j) geographical coordinates of node [i,j], i=1:m, j=1:n, [m,n]=size(long)=size(lat)
 102 %
 103 % output
 104 %   distance(i+1,j+1,a+2,b+2) = geographical distance between node [i,j] and [i+a,j+b] , a,b=-1:1
 105
 106 distance=zeros(size(long,1),size(long,2),3,3);
 107 i=2:size(long,1)-1;
 108 j=2:size(long,2)-1;
 109 for a=-1:1
 110    for b=-1:1
 111       %for j=2:size(long,2)-1
 112       %   for i=2:size(long,1)-1
 113       %      % distance between node [i,j] and [i+a,j+b]
 114             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);
 115       %   end
 116       %end
 117    end
 118 end
 119 end
 120