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