cycling/graph_dets.m

   1 function path_struct = graph_dets(w,cull)
   2 %w = read_wrfout_tign(wrfout)
   3 %cull - data reduction factor, reduces sizes of matrices
   4 %       example new_matrix = old_matrix(1:cull:end,1:cull:end)
   5
   6
   7 [fire_name,save_name,prefix] = fire_choice();
   8 red = subset_domain(w);
   9 time_bounds(2) = red.max_tign;
  10 time_bounds(1) = red.min_tign;
  11
  12 % figures
  13 fig.fig_map=0;
  14 fig.fig_3d=0;
  15 fig.fig_interp=0;
  16 p = sort_rsac_files(prefix);
  17
  18 g_str = 'g_graph.mat';
  19 if ~exist(g_str,'file')
  20     %loading L2 data
  21     g = subset_l2_detections(prefix,p,red,time_bounds,fig);
  22     save(g_str, 'g', '-v7.3');
  23 else
  24     g = [];
  25     reload_dets = input_num('Reload detections? 1 = yes',0);
  26     if reload_dets == 1
  27         g = subset_l2_detections(prefix,p,red,time_bounds,fig);
  28         save(g_str, 'g', '-v7.3');
  29     else
  30         load g_graph.mat;
  31     end
  32 end
  33
  34 close all
  35 pts = [];
  36 %minimum detection confidence level
  37 min_con = 70;
  38 %make unique ignition point
  39 for i = 1:length(g)% fprintf('Detections collected \n')
  40     % figure(1),scatter3(pts(:,2),pts(:,1),pts(:,3));
  41     % title('full scatter')
  42     if sum(g(i).det(3:5)) > 0
  43         fires = g(i).conf >= min_con;
  44         lons = mean(g(i).lon(fires));
  45         lats = mean(g(i).lat(fires));
  46         confs = mean(double(g(i).conf(fires)));
  47         times = g(i).time-0.25;
  48         frps = mean(g(i).power(fires));
  49         gran = i;
  50         pts = [lats,lons,times,confs,frps,gran];
  51
  52     end
  53     if norm(pts) > 0
  54         break
  55     end
  56 end
  57
  58 %can change end time for comparisons
  59 end_time = red.end_datenum;
  60 for i = 1:length(g)
  61     % don't use times after model end
  62     if (sum(g(i).det(3:5)) > 0) && (g(i).time < end_time)
  63         fires = g(i).conf >= min_con;
  64         lons = g(i).lon(fires);
  65         lats = g(i).lat(fires);
  66         times = g(i).time*ones(size(lons));
  67         confs = double(g(i).conf(fires));
  68         frps = g(i).power(fires);
  69         gran = i*ones(size(lons));
  70         pts = [pts;[lats',lons',times',confs',frps',gran']];
  71     end
  72 end
  73 % fprintf('Detections collected \n')
  74 % figure(1),scatter3(pts(:,2),pts(:,1),pts(:,3));
  75 % title('full scatter')
  76
  77 %prune the data
  78
  79 n_points = pts(1:cull:end,:,:,:,:,:);
  80 % figure(2),scatter3(n_points(:,2),n_points(:,1),n_points(:,3)-red.start_datenum);
  81 % title('Partial scatter')
  82
  83 %make edge weights
  84 n = length(n_points);
  85 %adjacency / distance matrix
  86 a = zeros(n,n);
  87 %velocity matrix
  88 v = a;
  89 %time matrix
  90 t = a;
  91 %cone volume matrix
  92 cv = a;
  93 %%% figure out way to get max_t automatically
  94 % maximum allowed time between nodes in the graph to allow them to be
  95 % connected
  96 max_t = 1.9*(24*3600);
  97
  98 %maybe change later
  99 pts = n_points;
 100
 101 %convert from points in the scattered data of L2 data to nearest
 102 %neighbors on the fire grid
 103 grid_pts = fixpoints2grid(red,n_points);
 104
 105
 106 %% computing distance between points using GPS coords
 107 E = wgs84Ellipsoid;distance(ig_point,i_point,E);
 108 %max distance from ignition
 109 max_d = 0;
 110 %error in time of fire from time of detection
 111 time_err = 0.2;
 112 distant_point = 1;
 113 ig_point = [pts(1,1),pts(1,2)];
 114 for i = 1:n
 115     time = pts(i,3);
 116     %%% local time for help in figuring out day/night
 117     locs(i) = local_time(time);
 118     i_point = [pts(i,1),pts(i,2)];
 119     %find furthest detection from ignition
 120     new_d = distance(ig_point,i_point,E);
 121     if new_d > max_d
 122         max_d = new_d;
 123         distant_point = i;
 124     end
 125     %distance from all points
 126     for j = 1:n
 127         time_diff = max(time_err,pts(j,3)-time)*(24*3600);
 128         %if (time_diff > 0  && time_diff < max_t)
 129
 130         %% maybe fix this, time_err is weird
 131         t(i,j) = time_diff;  %pts(j,3)-time;
 132         if time_diff < max_t
 133             j_point = [pts(j,1),pts(j,2)];
 134             a(i,j) = distance(i_point,j_point,E);
 135             v(i,j) = a(i,j)/time_diff;
 136         end
 137     end
 138 end
 139
 140 %fix up triangular matrices
 141 t = t-t';
 142 t_mask = t < -time_err*(24*3600);
 143 t(~t_mask) = abs(t(~t_mask));
 144 a = a+a';
 145 raw_dist = a;
 146 %a(t_mask)=0;
 147 % work on the v matrix?
 148 fprintf('Matrix ready \n');
 149 %scatter ignition and distant point
 150 % figure(2),hold on,zlabel('Days from start'),xlabel('Lon'),ylabel('Lat')
 151 % scatter3(pts(1,2),pts(1,1),pts(1,3)-red.start_datenum,'*r');
 152 % scatter3(pts(distant_point,2),pts(distant_point,1),pts(distant_point,3)-red.start_datenum,'*r');
 153 % hold off
 154
 155 %calculate max velocity of fire
 156 v_std = std(v(:)); %change to v(v>0)?
 157 v_mean = mean(v(:));
 158 %max_v = abs(max_d/(pts(distant_point,3)-pts(1,3)));
 159 max_v = v_mean;
 160 v_mask = v < max_v;
 161 %filter detections too far apart, fix this KLUGE
 162
 163 t1 = 0.9*(24*3600); %1.2 days for fast growth cone
 164 speed_penalty = 1.2;
 165 fast = 1.6;
 166 slow = 0.4;
 167 speeders = 0;
 168 for i = 1:n
 169     if (locs(i) > 8) && (locs(i) < 18)
 170         sp1 = fast;
 171     else
 172         sp1 = slow;
 173     end
 174     for j = 1:n
 175         if (locs(j) > 8) && (locs(j) < 18)
 176             sp2 = fast;
 177         else
 178             sp2 = slow;
 179         end
 180         speed_penalty = (sp1+sp2)/2;
 181         %fprintf('speed penalty = %0.2f \n',speed_penalty)
 182         % lft haf of speed cone
 183         if (a(i,j) > 0) && (v(i,j) > speed_penalty*max_v) && (t(i,j) < t1)
 184             a(i,j) = 0;
 185             speeders = speeders + 1;
 186         end
 187         %right half of speed cone
 188         max_r = max_v*t1;
 189         cone_slope = -max_r/(max_t-t1);
 190         if (a(i,j) > max_r-cone_slope*(t(i,j)-t1)) && (t(i,j) >= t1)
 191             a(i,j) = 0;
 192             speeders = speeders + 1;
 193             %fprintf('Second type speeder removed \n')
 194         end
 195     end
 196 end
 197 fprintf('%d speeders removed \n',speeders)
 198
 199 %% make directed graph from a
 200 fg = digraph(a);
 201 figure(3),plot(fg);
 202 path_count = 0;
 203
 204 % new_points = pts;
 205 start_pt = 1;
 206 for i = 1:start_pt
 207     for j = 1:cull:n
 208         distant_point = j;
 209         % finds shortest path between points i,j
 210         % p is the points in the path, d is the total distance
 211         [p,d] = shortestpath(fg,i,j);
 212         if ~isempty(p)
 213             path_count = path_count + 1;
 214             paths(path_count).p = p;
 215             paths(path_count).d = d;
 216             %path confidenc is geometric mean of detections in path
 217             paths(path_count).c = prod(pts(p,4))^(1/length(p));
 218             %fprintf('%d points in path \n',length(p))
 219         end
 220         figure(2),hold on
 221         %l2 points
 222         %plot3(pts(p,2),pts(p,1),pts(p,3)-red.start_datenum,':r');
 223         %grid points
 224         plot3(grid_pts(p,4),grid_pts(p,3),pts(p,3)-red.start_datenum,'g');
 225                 for k = 1:length(p)
 226                     scatter3(pts(p(k),2),pts(p(k),1),pts(p(k),3)-red.start_datenum,'*r');
 227                 end
 228         hold off
 229 %         % add a new point to the list by interpolation
 230 %         for k = 1:length(p)-1
 231 %             new_pt = ([pts(p(k),1),pts(p(k),2),pts(p(k),3)]+[pts(p(k+1),1),pts(p(k+1),2),pts(p(k+1),3)])/2;
 232 %             new_points = [new_points;new_pt];
 233 %         end
 234     end
 235     path_struct.raw_dist = raw_dist;
 236     path_struct.paths = paths;
 237     path_struct.graph = fg;
 238     path_struct.distances = a;
 239     path_struct.speeds = v;
 240     path_struct.points = double(pts);
 241     path_struct.red = red;
 242 %    path_struct.new_points = new_points;
 243     path_struct.grid_pts = grid_pts(:,3:4);
 244     path_struct.idx = grid_pts(:,1:2);
 245 end
 246 end
 247
 248 function lt = local_time(time)
 249     shift = 7;
 250     lt = 24*((time-shift)-floor(time-shift));
 251 end
 252
 253