cycling/tign_try.m

   1 function path_struct = tign_try(w,cpd,grid_dist)
   2 %scatter fake data on TIGN cone, try to recover the shape
   3 %cpd clusters per day of data
   4 %[fire_name,save_name,prefix] = fire_choice();
   5 clear_q = input_num('Delete mat files? 1 = yes',1,1);
   6 if clear_q
   7     !rm red.mat pts.mat
   8     %!rm pts.mat
   9 end
  10
  11 if ~exist('red.mat','file')
  12     red = subset_domain(w,1);
  13     save red.mat red
  14 else
  15     load red.mat
  16 end
  17
  18 %compute grid sizes
  19 E = wgs84Ellipsoid;
  20 dlon= distance(red.min_lat,red.min_lon,red.min_lat,red.max_lon,E);
  21 dlat= distance(red.min_lat,red.min_lon,red.max_lat,red.min_lon,E);
  22 if ~exist('grid_dist','var')
  23     grid_dist = 250;
  24 end
  25 new_m = round(dlon/grid_dist);
  26 new_n = round(dlat/grid_dist);
  27
  28 time_bounds(2) = red.max_tign;
  29 time_bounds(1) = red.min_tign;
  30
  31 %shrink the size for large matrices
  32 target_size = max(new_m,new_n);
  33
  34 if max(size(red.tign)) > target_size
  35
  36     [m,n] = size(red.tign);
  37     %shrink_factor
  38     %sf = 4;
  39     max_dim = max(m,n);
  40     sf = target_size/max_dim;
  41     n =round(n*sf);
  42     m = round(m*sf);
  43     red_copy = red;
  44     red = subset_small(red,m,n);
  45 end
  46
  47 % figures
  48 fig.fig_map=0;
  49 fig.fig_3d=0;
  50 fig.fig_interp=0;
  51
  52 [m,n] = size(red.tign);
  53 rm = 0.1;
  54 [ig_x,ig_y] = find(red.tign == min(red.tign(:)));
  55 ig_x = ig_x(1);
  56 ig_y = ig_y(1);
  57
  58 pts = [red.fxlat(ig_x,ig_y),red.fxlong(ig_x,ig_y),min(red.tign(:)),100,100,1];
  59 idx = [ig_x,ig_y];
  60
  61 % k can be expanded to loop through discrete lime levels like MODIS, VIIRS
  62 for k = 1 : 1
  63     % i,j are indices of the mesh
  64     for i = 1:m
  65         for j = 1:n
  66             if red.tign(i,j) < red.end_datenum - rm
  67                 if rand < cpd/100
  68                     %pts = [pts;[lats',lons',times',confs',frps',gran']];
  69                     pts = [pts;[red.fxlat(i,j),red.fxlong(i,j),red.tign(i,j),100,100,2*round(red.tign(i,j)-red.start_datenum)]];
  70                     idx = [idx;[i,j]];
  71                 end
  72             end
  73         end
  74     end
  75 end
  76 if length(pts) > 500
  77    pt_length = length(pts);
  78    msk = rand(pt_length,1);
  79    msk = msk > 500/pt_length;
  80    pts(msk,:) = [];
  81    idx(msk,:) = [];
  82 end
  83
  84 %add perimeter, optional
  85 [new_pts,new_idx] = make_perims(red.tign,red);
  86 pts = [pts;new_pts];
  87 idx = [idx;new_idx];
  88
  89
  90 %discretize the time
  91 dis_time = 0;
  92 if dis_time ==1
  93     sim_days = time_bounds(2)-time_bounds(1);
  94     over_passes = round(sim_days*3);
  95     time_bin = linspace(time_bounds(1),time_bounds(2),over_passes)+1/10*randn(1,over_passes);
  96     for i = 2:over_passes
  97         for j = 1:length(pts)
  98             if pts(j,3) > time_bin(i-1) && pts(j,3) <= time_bin(i)
  99                 pts(j,3) = time_bin(i);
 100             end
 101         end
 102     end
 103 end
 104
 105 %use same set of points from previous example / multi step
 106 % if exist('pts.mat','file')
 107 %     clear pts idx
 108 %     load pts.mat
 109 % else
 110 %     save pts.mat pts idx
 111 % end
 112
 113 figure,mesh(red.fxlong,red.fxlat,red.tign-red.start_datenum),hold on
 114 scatter3(pts(:,2),pts(:,1),pts(:,3)-red.start_datenum,'r*')
 115 [st1,st2]=sort(pts(:,3));
 116 pts(:,1) = pts(st2,1);
 117 pts(:,2) = pts(st2,2);
 118 pts(:,3) = pts(st2,3);
 119 idx(:,1) = idx(st2,1);
 120 idx(:,2) = idx(st2,2);
 121
 122 cull = 1;
 123 %cluster the data
 124 t1 = red.min_tign;
 125 t2 = red.max_tign;
 126 %clusters = round((t2-t1)*cpd);
 127 %clusters = round(length(pts)/20);
 128 %clusters = min(cpd,length(pts(:,3)));
 129 clusters = 20;
 130 [s_idx,s_c] = kmeans(pts(:,1:2),clusters);
 131 %scatter the clusters with coloring
 132 figure,scatter3(pts(s_idx==1,2),pts(s_idx==1,1),pts(s_idx==1,3)-red.start_datenum);
 133 hold on
 134 for i = 2:clusters
 135   scatter3(pts(s_idx==i,2),pts(s_idx==i,1),pts(s_idx==i,3)-red.start_datenum);
 136 end
 137 hold off
 138
 139 %make edge weights
 140 n = length(pts);
 141 %adjacency / distance matrix
 142 a = zeros(n,n);
 143 %velocity matrix
 144 v = a;
 145 %time matrix
 146 t = a;
 147 %cone volume matrix
 148 %cv = a;
 149 grid_pts = pts(:,1:2);
 150
 151 %distance matrices..
 152 %% computing distance between points using GPS coords
 153 E = wgs84Ellipsoid;
 154
 155 %make cluster center distance matrix
 156 clust_dist=zeros(clusters);
 157 %needed ????
 158 for i = 1:clusters
 159     i_clust = [s_c(i,1),s_c(i,2)];
 160     for j = 1:clusters
 161         j_clust = [s_c(j,1),s_c(j,2)];
 162         clust_dist(i,j) =  distance(i_clust,j_clust,E);
 163     end
 164 end
 165
 166 %max distance from ignition
 167 max_d = 0;
 168 %error in time of fire from time of detection
 169 time_err = 0.2;
 170 distant_point = 1;
 171 ig_point = [pts(1,1),pts(1,2)];
 172 for i = 1:n
 173     time = pts(i,3);
 174     %%% local time for help in figuring out day/night
 175     locs(i) = local_time(time);
 176     i_point = [pts(i,1),pts(i,2)];
 177     %find furthest detection from ignition
 178     new_d = distance(ig_point,i_point,E);
 179     if new_d > max_d
 180         max_d = new_d;
 181         distant_point = i;
 182     end
 183     %distance from all points
 184     for j = 1:n
 185         time_diff = (pts(j,3)-time)*(24*3600);
 186
 187         t(i,j) = time_diff;  %pts(j,3)-time;
 188         if time_diff > 0
 189             j_point = [pts(j,1),pts(j,2)];
 190             a(i,j) = distance(i_point,j_point,E);
 191             v(i,j) = a(i,j)/max(time_diff,0.1);
 192         end
 193     end
 194 end
 195 raw_dist = a;
 196
 197 %start filtering distance
 198 cluster_mult = 0.25;
 199 for i = 1:n
 200     for j = 1:n
 201         %         % make points in same cluster close
 202         if (a(i,j) > 0) && (s_idx(i) == s_idx(j)) %in same cluster
 203             a(i,j) = cluster_mult*a(i,j);
 204         end
 205         % make points in different clusters further apart
 206 %         if (a(i,j) > 0) && (s_idx(i)~=s_idx(j))
 207 %             a(i,j) = a(i,j) + clust_dist(s_idx(i),s_idx(j));
 208 %         end
 209     end
 210 end
 211
 212
 213 %make paths
 214 fg = digraph(a);
 215 %figure(3),plot(fg);
 216 path_count = 0;
 217
 218 % new_points = pts;
 219 start_pt = 1;
 220 figure(222),close(gcf)
 221 for i = 1:start_pt
 222     for j = 1:cull:n
 223         % finds shortest path between points i,j
 224         % p is the points in the path, d is the total distance
 225         [p,d] = shortestpath(fg,i,j);
 226         if ~isempty(p)
 227             path_count = path_count + 1;
 228             paths(path_count).p = p;
 229             paths(path_count).d = d;
 230             %path confidenc is geometric mean of detections in path
 231             paths(path_count).c = prod(pts(p,4))^(1/length(p));
 232             %fprintf('%d points in path \n',length(p))
 233         end
 234         figure(222),hold on
 235         %l2 points
 236         %plot3(pts(p,2),pts(p,1),pts(p,3)-red.start_datenum,':r');
 237         %grid points
 238         plot3(pts(p,2),pts(p,1),pts(p,3)-red.start_datenum,'g');
 239         for k = 1:length(p)
 240             scatter3(pts(p(k),2),pts(p(k),1),pts(p(k),3)-red.start_datenum,'*r');
 241         end
 242         hold off
 243         %         % add a new point to the list by interpolation
 244         %         for k = 1:length(p)-1
 245         %             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;
 246         %             new_points = [new_points;new_pt];
 247         %         end
 248     end
 249
 250 end
 251
 252 path_struct.raw_dist = raw_dist+raw_dist';
 253 path_struct.paths = paths;
 254 path_struct.graph = fg;
 255 path_struct.distances = a;
 256 path_struct.speeds = v;
 257 path_struct.points = double(pts);
 258 path_struct.red = red;
 259 %    path_struct.new_points = new_points;
 260 path_struct.grid_pts = grid_pts;
 261 path_struct.idx = idx;
 262
 263
 264
 265
 266 end
 267
 268
 269 function lt = local_time(time)
 270     shift = 7;
 271     lt = 24*((time-shift)-floor(time-shift));
 272 end