cycling/cluster_paths.m

   1 function path_struct = cluster_paths(w,cull,grid_dist)
   2 % assign shortest paths, using clustering
   3 % inputs - w = read_wrfout_tign(f)
   4 %          cull - number for using smaller data sets cull = 1 --> use all
   5 %          data points.
   6 % output - cp , struct with path info
   7
   8 [fire_name,save_name,prefix,perim] = fire_choice();
   9 red = subset_domain(w);
  10 multi = input_num('Use multigrid? 1 = yes',0,1);
  11 if multi
  12     if exist('ps_multi.mat','file')
  13         load ps_multi.mat
  14         load an_multi.mat
  15         red2=ps_multi.red;
  16         new_an = interp2(red2.fxlat,red2.fxlong,an_multi,red.fxlat,red.fxlong);
  17         red.tign = new_an;
  18     end
  19 end
  20 %compute grid sizes
  21 E = wgs84Ellipsoid;
  22 dlon= distance(red.min_lat,red.min_lon,red.min_lat,red.max_lon,E);
  23 dlat= distance(red.min_lat,red.min_lon,red.max_lat,red.min_lon,E);
  24 if ~exist('grid_dist','var')
  25     grid_dist = 250;
  26 end
  27 new_m = round(dlon/grid_dist);
  28 new_n = round(dlat/grid_dist);
  29
  30 %shrink the size for large matrices
  31 target_size = max(new_m,new_n);
  32
  33 if max(size(red.tign)) > target_size
  34
  35     [m,n] = size(red.tign);
  36     %shrink_factor
  37     %sf = 4;
  38     max_dim = max(m,n);
  39     sf = target_size/max_dim;
  40     n =round(n*sf);
  41     m = round(m*sf);
  42     red_copy = red;
  43     red = subset_small(red,m,n);
  44 end
  45 time_bounds(2) = red.max_tign;
  46 time_bounds(1) = red.min_tign;
  47 new_end_time = input_num('Use alternate end time? Enter number of extra days, 0 if no.',0,1);
  48 if new_end_time ~=0
  49   time_bounds(2) = time_bounds(2)+new_end_time;
  50   red.max_tign = time_bounds(2);
  51   red.end_datenum = time_bounds(2);
  52   fprintf('New end time %s\n',datestr(time_bounds(2)));
  53 end
  54 % time_bounds(2) = 7.354591409722222e+05;
  55
  56 % figures
  57 fig.fig_map=0;
  58 fig.fig_3d=0;
  59 fig.fig_interp=0;
  60 p = sort_rsac_files(prefix);
  61 %time_bounds(2) = p.time(end);
  62 %time_bounds(1) = p.time(1);
  63
  64 %load satellite detection data
  65 g_str = 'g_cluster.mat';
  66 if ~exist(g_str,'file')
  67     %loading L2 data
  68     g = subset_l2_detections(prefix,p,red,time_bounds,fig);
  69     save(g_str, 'g', '-v7.3');
  70 else
  71     g = [];
  72     reload_dets = input_num('Reload detections? 1 = yes',1,1);
  73     if reload_dets == 1
  74         g = subset_l2_detections(prefix,p,red,time_bounds,fig);
  75         save(g_str, 'g', '-v7.3');
  76     else
  77         load g_cluster.mat;
  78     end
  79 end
  80
  81 %load satellite ground detection data
  82 % get fire mask, fxlong, fxlat for each granule
  83 %pos_detects = collect_pos(prefix,p,red,time_bounds,fig)
  84
  85 %add functionality to pull in perimeter data here
  86 use_perims = input_num('Use perimeter data ? 1 = yes',0);
  87 if use_perims == 1
  88     %use just 40 points per peimeter
  89     p_points = input_num('How many perimeter points to use?',20);
  90     %p_points = 20;
  91     p_gran = perim2gran(p_points,perim);
  92     interp_perim = input_num('Interpolate perimeters to grid? yes = 1',1)
  93     if interp_perim == 1
  94        for i =1 length(p_gran)
  95           pts = [p_gran(i).lat',p_gran(i).lon'];
  96           n_pts = fixpoints2grid(w,pts);
  97           n_pts = unique(n_pts,'rows');
  98           l = length(n_pts);
  99           p_gran(i).power = 50*ones(1,l);
 100           p_gran(i).data  =  9*ones(1,l);
 101           p_gran(i).conf  = 95*ones(1,l);
 102           p_gran(i).lat = n_pts(:,3)';
 103           p_gran(i).lon = n_pts(:,4)';
 104           p_gran(i).mask = [];
 105        end
 106     end
 107     gl = length(g);
 108     rm_idx = zeros(1,length(p_gran));
 109     for i = 1:length(p_gran)
 110         %only add perimeters up to final granules time
 111         if p_gran(i).time < time_bounds(2);% g(gl).time
 112             g(length(g)+1)=p_gran(i);
 113         else
 114             fprintf('Perimeter time after simulation end, removing from set of perimeters \n')
 115            rm_idx(i) = 1;
 116         end
 117     end
 118     rm_idx = logical(rm_idx);
 119     p_gran(rm_idx) = [];
 120     %sort the data by time
 121     T = struct2table(g);
 122     sortedT = sortrows(T,'time');
 123     g = table2struct(sortedT);
 124     %select only a specified perimeter, delete data after - use for
 125     %    initializing a fire from a specified perimeter
 126     spec_perim = input_num('Specify a perimeter? 1 = yes',0)
 127     if spec_perim == 1
 128        for i = 1:length(p_gran)
 129            fprintf('%d  %s  \n',i,p_gran(i).file)
 130        end
 131        perim_num = input_num('Which perimeter to use? ',1);
 132        %delete granules past perimeter
 133        for i = length(g):-1:1
 134            %fprintf('%d Time diff: %f \n',i, g(i).time - p_gran(perim_num).time)
 135            if g(i).time > p_gran(perim_num).time
 136                g(i) = [];
 137            end
 138        end
 139        %filter points outside of perimeter make low confidence so they are
 140        %ignore in the  graph
 141        gl = length(g);
 142        for i = 1:gl-1
 143           in = inpolygon(g(i).lon,g(i).lat,g(gl).lon,g(gl).lat);
 144           scatter(g(i).lon,g(i).lat)
 145           hold on, scatter(g(gl).lon,g(gl).lat)
 146           g(i).conf(~in) = 20;
 147        end
 148     end
 149 end
 150
 151
 152 %close all
 153 pts = [];
 154 %minimum detection confidence level
 155 min_con = 70;
 156 %make unique ignition point
 157 for i = 1:length(g)% fprintf('Detections collected \n')
 158     % figure(1),scatter3(pts(:,2),pts(:,1),pts(:,3));
 159     % title('full scatter')
 160     if sum(g(i).det(3:5)) > 0
 161         fires = g(i).conf >= min_con;
 162         lons = mean(g(i).lon(fires));
 163         lats = mean(g(i).lat(fires));
 164         confs = mean(double(g(i).conf(fires)));
 165         times = g(i).time-0.05;
 166         frps = mean(g(i).power(fires));
 167         gran = i;
 168         pts = [lats,lons,times,confs,frps,gran];
 169     end
 170     if norm(pts) > 0
 171         break
 172     end
 173 end
 174
 175 %can change end time for comparisons
 176 if new_end_time ~= 0
 177     end_time = time_bounds(2);
 178 else
 179     end_time = red.max_tign;
 180 end
 181 for i = 1:length(g)
 182     % don't use times after model end
 183     if (sum(g(i).det(3:5)) > 0) && (g(i).time < end_time)
 184         fires = g(i).conf >= min_con;
 185         lons = g(i).lon(fires);
 186         lats = g(i).lat(fires);
 187         times = g(i).time*ones(size(lons));
 188         confs = double(g(i).conf(fires));
 189         frps = g(i).power(fires);
 190         gran = i*ones(size(lons));
 191         pts = [pts;[lats',lons',times',confs',frps',gran']];
 192     end
 193 end
 194
 195 %prune the data
 196 n_points = pts(1:cull:end,:,:,:,:,:);
 197 %
 198 %filter Nan fom data
 199 %should be handled in the perim2gran.m function
 200 % for i = length(n_points):-1:1
 201 %     if sum(isnan(n_points(i,:)))~= 0
 202 %         n_points(i,:) = [];
 203 %     end
 204 % end
 205
 206
 207 %% for computing distance between points using GPS coords
 208 % also used for finding aspect of the slope, for clustering
 209 E = wgs84Ellipsoid;
 210 %[aspect,slope,dy,dx] = gradientm(red.fxlat,red.fxlong,red.fhgt,E);
 211 clst_pts =  fixpoints2grid(red,n_points);
 212 % just use the index numbers, maintain the l2 data coords
 213 clst_pts(:,3:4) = n_points(:,1:2);
 214 %ig_pt = [mean(clst_pts(:,3)),mean(clst_pts(:,4))];
 215    ig_pt = [clst_pts(1,3),clst_pts(1,4)];
 216 for i = 1:length(clst_pts)
 217    pt_1 = [ig_pt(1,1),clst_pts(i,4)];
 218    pt_2 = [clst_pts(i,3),clst_pts(i,4)];
 219    %distances in lon and lat directions, with sign
 220    d_lon = -sign(clst_pts(1,4)-clst_pts(i,4))*distance(ig_pt,pt_1,E);
 221    d_lat = -sign(clst_pts(1,3)-clst_pts(i,3))*distance(pt_2,pt_1,E);
 222    cp(i,:) = [clst_pts(i,:),d_lat,d_lon];
 223     %% work out x-y coordinate with pt 1 as origin
 224 end
 225
 226 %remove data points too far from the main set,
 227 %cluster pts into 2 clusters
 228 [s_idx2,s_c2] = kmeans(cp(:,5:6),2);
 229 %find cluster with smallest number of pts
 230 c1 = sum(s_idx2 == 1);
 231 c2 = sum(s_idx2 == 2);
 232 fprintf('Two clusters computed %d and %d points in them \n',c1,c2)
 233 small_clust = 1;
 234 if c2 < c1
 235     small_clust = 2;
 236 end
 237 if norm(s_c2(2,:)) > 1e4; %sum(s_idx2==small_clust)/(c1+c2) < 0.05
 238     cp(s_idx2==small_clust,:) = [];
 239     pts(s_idx2==small_clust,:) = [];
 240     n_points(s_idx2==small_clust,:) = [];
 241     clst_pts(s_idx2==small_clust,:) = [];
 242
 243 end
 244
 245
 246
 247
 248 %cluster the data
 249 dt = 3*ceil(g(end).time - g(1).time);
 250 space_clusters = 20; %days
 251 %more clusters for using perimeter data
 252 if use_perims == 1
 253     space_clusters = dt;
 254 end
 255
 256 %[s_idx,s_c] = kmeans(pts(:,1:2),space_clusters);
 257 %clustering using aspect, not good
 258 [s_idx,s_c] = kmeans(cp(:,5:6),space_clusters);
 259
 260 % find optimal cluster k
 261 % max_clusts = 20;
 262 % s_pts = pts(:,1:2);
 263 % klist=2:max_clusts;%the number of clusters you want to try
 264 % myfunc = @(X,K)(kmeans(X, K));
 265 % eva = evalclusters(s_pts,myfunc,'CalinskiHarabasz','klist',klist)
 266 % classes=kmeans(s_pts,eva.OptimalK);
 267 % dt = eva.OptimalK;
 268
 269 % %spatial clusters scatter plot
 270 % figure,scatter(pts(s_idx==1,1),pts(s_idx==1,2));
 271 % hold on% dt = round(g(end).time - g(1).time);% dt = round(g(end).time - g(1).time);
 272 % space_clusters = dt; %days
 273 % [s_idx,s_c] = kmeans(pts(:,1:2),space_clusters);
 274
 275
 276 % for i = 2:dt
 277 %   scatter(pts(s_idx==i,1),pts(s_idx==i,2));
 278 % end
 279 % hold off
 280
 281 %scatter 3d  lat/lon
 282 % figure,scatter3(pts(s_idx==1,2),pts(s_idx==1,1),pts(s_idx==1,3));
 283 % hold on
 284 % for i = 2:dt
 285 %   scatter3(pts(s_idx==i,2),pts(s_idx==i,1),pts(s_idx==i,3));
 286 % end
 287 % hold off
 288
 289 %scatter 3d  ldistances in lat/lon directions
 290 figure(177),scatter3(cp(s_idx==1,6),cp(s_idx==1,5),pts(s_idx==1,3)-red.start_datenum);
 291 title('Clustering of Data')
 292 xlabel('East-West Distance [m]')
 293 ylabel('North-South Distance [m]')
 294 zlabel('Time [days]')
 295 hold on
 296 for i = 2:space_clusters
 297   scatter3(cp(s_idx==i,6),cp(s_idx==i,5),pts(s_idx==i,3)-red.start_datenum);
 298 end
 299 hold off
 300
 301 %make edge weights
 302 n = length(n_points);
 303 %adjacency / distance matrix
 304 a = zeros(n,n);
 305 %velocity matrix
 306 v = a;
 307 %time matrix
 308 t = a;
 309 %cone volume matrix
 310 %cv = a;
 311 %%% figure out way to get max_t automatically
 312 % maximum allowed time between nodes in the graph to allow them to be
 313 % connected
 314 %max_t = 1.9*(24*3600);
 315
 316 %maybe change later
 317 pts = n_points;
 318
 319 %convert from points in the scattered data of L2 data to nearest
 320 %neighbors on the fire grid
 321 grid_pts = fixpoints2grid(red,n_points);
 322
 323 %% computing distance between points using GPS coords
 324
 325 %make cluster center distance matrix
 326 clust_dist=zeros(space_clusters);
 327 for i = 1:space_clusters
 328     i_clust = [s_c(i,1),s_c(i,2)];
 329     for j = 1:space_clusters
 330         j_clust = [s_c(j,1),s_c(j,2)];
 331         clust_dist(i,j) =  sqrt((i_clust(1)-j_clust(1))^2+((i_clust(2)-j_clust(2))^2));
 332     end
 333 end
 334
 335 %max distance from ignition
 336 max_d = 0;
 337 %error in time of fire from time of detection
 338 %time_err = 0.2;
 339 %distant_point = 1;
 340 ig_point = [pts(1,1),pts(1,2)];
 341 for i = 1:n
 342     time = pts(i,3);
 343     %%% local time for help in figuring out day/night
 344     locs(i) = local_time(time);
 345     i_point = [pts(i,1),pts(i,2)];
 346     %find furthest detection from ignition
 347     new_d = distance(ig_point,i_point,E);
 348     if new_d > max_d
 349         max_d = new_d;
 350         distant_point = i;
 351     end
 352     %distance from all points
 353     for j = 1:n
 354         time_diff = (pts(j,3)-time)*(24*3600);
 355
 356         t(i,j) = time_diff;  %pts(j,3)-time;
 357         if time_diff > 0
 358             j_point = [pts(j,1),pts(j,2)];
 359             a(i,j) = distance(i_point,j_point,E);
 360             v(i,j) = a(i,j)/max(time_diff,0.1);
 361             if a(i,j) > 2e4
 362                 %fprintf('Points far apart  %d and %d\n',i,j)
 363                 a(i,j) = 0;
 364                 v(i,j) = 0;
 365             end
 366 %             if v(i,j) > 2
 367 %                 fprintf('fast ROS beteween points %d and %d\n',i,j)
 368 %             end
 369         end
 370     end
 371 end
 372
 373 %fix up triangular matrices
 374 %t_mask = t >= 0;
 375 % a = a+a';
 376 raw_dist = a;
 377 fprintf('matrices done\n')
 378
 379 %start filtering distances and change those within same cluster
 380 cluster_mult = 0.25;
 381 for i = 1:n
 382     for j = 1:n
 383 %         % make points in same cluster close
 384         if (a(i,j) > 0) && (s_idx(i) == s_idx(j)) %in same cluster
 385             a(i,j) = cluster_mult*a(i,j);
 386         end
 387         % make points in different clusters further apart
 388         % if (a(i,j) > 0) && (s_idx(i)~=s_idx(j))
 389         %    a(i,j) = a(i,j) + clust_dist(s_idx(i),s_idx(j));
 390         % end
 391     end
 392 end
 393
 394
 395 %make paths
 396 fg = digraph(a);
 397 % figure(3),plot(fg);
 398 path_count = 0;
 399
 400 % new_points = pts;
 401 start_pt = 1;
 402 for i = 1:start_pt
 403     for j = 1:cull:n
 404         % finds shortest path between points i,j
 405         % p is the points in the path, d is the total distance
 406         [p,d] = shortestpath(fg,i,j);
 407         if ~isempty(p)
 408             path_count = path_count + 1;
 409             paths(path_count).p = p;
 410             paths(path_count).d = d;
 411             %path confidenc is geometric mean of detections in path
 412             paths(path_count).c = prod(pts(p,4))^(1/length(p));
 413             %fprintf('%d points in path \n',length(p))
 414         end
 415         %only plot if n is less than 400
 416         if n < 400
 417             figure(178),hold on
 418             title('Shortest Paths')
 419             xlabel('Lon')
 420             ylabel('Lat')
 421             zlabel('Time [days]')
 422             %l2 points
 423             %plot3(pts(p,2),pts(p,1),pts(p,3)-red.start_datenum,':r');
 424             %grid points
 425             plot3(grid_pts(p,4),grid_pts(p,3),pts(p,3)-red.start_datenum,'g');
 426             for k = 1:length(p)
 427                 scatter3(pts(p(k),2),pts(p(k),1),pts(p(k),3)-red.start_datenum,'*r');
 428             end
 429             hold off
 430         end
 431         %         % add a new point to the list by interpolation
 432         %         for k = 1:length(p)-1
 433         %             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;
 434         %             new_points = [new_points;new_pt];
 435         %         end
 436     end
 437     path_struct.raw_dist = raw_dist+raw_dist';
 438     path_struct.paths = paths;
 439     path_struct.graph = fg;
 440     path_struct.distances = a;
 441     path_struct.speeds = v;
 442     path_struct.points = double(pts);
 443     path_struct.red = red;
 444     %    path_struct.new_points = new_points;
 445     path_struct.grid_pts = grid_pts(:,3:4);
 446     path_struct.idx = grid_pts(:,1:2);
 447     %cluster information
 448     path_struct.s_idx = s_idx;
 449     path_struct.s_c = s_c;
 450 end
 451
 452
 453 end %function
 454
 455 function lt = local_time(time)
 456     shift = 7;
 457     lt = 24*((time-shift)-floor(time-shift));
 458 end