cycling/squish4.m

   1 function tign_new = squish4(ps,a,b,gq)
   2 %% ps is struct with paths, red, graph,distances, etc.
   3 % a, b - parameters of the function smooth up function
   4 %% ps = graph_dets(w)
   5
   6 if ~exist('a','var')
   7     a = 100;
   8     b = 1/3;
   9 end
  10
  11 %combine l2 time with detection points fixed to grid
  12 pts = ps.grid_pts;
  13 pts(:,3) = ps.points(:,3);
  14 pts(1,3) = pts(1,3)-1/2;
  15 %forecast
  16 tign = ps.red.tign;
  17 % t_max = max(tign(:))-.1;
  18 % tign(tign>=t_max) =t_max+1;
  19 t0 = min(tign(:));
  20 [m,n]=size(ps.red.tign);
  21
  22 % flat_tign = input_num('Use flat start? 1 = yes',0)
  23 % if flat_tign
  24 %     tign = ps.red.end_datenum*ones(m,n);
  25 % end
  26
  27 %%% data likelihood spline %%%
  28 %set up data likelihood spline
  29 [p_like_spline,~,~] = make_spline(100,1000);
  30 %test
  31 % t = 4*(-10:20)*3600;
  32 % ts = p_like_spline(t);
  33 % figure,plot(t,exp(ts))
  34 % hold on,plot(t,1-exp(ts))
  35 %make vector of data likelihoods
  36 if ~exist('gq','var')
  37     gq = input_num('use ground detections? yes = [1]',1,1);
  38 end
  39 if gq
  40     smooth_ground = 1;%max(m,n)/50;
  41     fprintf('Collecting ground detection data\n')
  42     %grounds = ground_detects(ps.red);
  43     %[jgrid,igrid]=meshgrid([1:length(ps.red.jspan)]',[1:length(ps.red.ispan)]');
  44     [jgrid,igrid]=meshgrid([1:n]',[1:m]');
  45     %make polygon around detections
  46     %infire = inpolygon(grounds.land(:,4),grounds.land(:,3),pts(:,2),pts(:,1));
  47     %infire = inpolygon(grounds.land(:,4),grounds.land(:,3),grounds.land(infire,4),grounds.land(infire,3));
  48     infire = inpolygon(igrid(:),jgrid(:),ps.idx(:,1),ps.idx(:,2));
  49     %remove holes in mask
  50     infire = inpolygon(igrid(:),jgrid(:),igrid(infire),jgrid(infire));
  51 end
  52 %try
  53
  54 tmax=max(tign(:));
  55 tign_ground = tign;
  56 tign_flat = ps.red.end_datenum*ones(size(ps.red.tign));
  57 %figure(159),hold on,scatter(pts(:,2),pts(:,1),'*r')
  58 time_err = 0.0;
  59 g_diff = (tign_flat-tign+time_err)*24*3600;
  60 g_likes = p_like_spline(g_diff);
  61 beta = 1/10;
  62 beta_vect = exp(g_likes);
  63 use_beta_likes = 1;
  64 %%% ground detection likelihood
  65 % g_times = zeros(pts_length,1);
  66 % for i = 1:pts_length
  67 %     g_times(i) = ps.red.tign(ps.idx(i,1),ps.idx(i,2));
  68 % end
  69 ground_steps = 1;
  70 if gq
  71     data_area = sum(infire);
  72     for i = 1:ground_steps
  73         fire_mask = tign_ground < tmax-0.1;
  74         fire_area = sum(fire_mask(:));
  75         out_side = fire_mask(:)-infire;
  76         out_side(out_side<0)=0;
  77         out_sum(i) = sum(out_side(:));
  78         out_fraction = out_sum(i)/data_area;
  79         if i > 2 && out_sum(i) == out_sum(i-1)
  80             smooth_ground = 0.9*smooth_ground;
  81             fprintf('decreasing ground smoothing \n')
  82         end
  83         if out_fraction < 0.1
  84             break
  85         end
  86         fprintf('Fire area: %f pixel_out/data area: %f  pixels_out %f\n',fire_area,out_sum(i)/data_area,out_sum(i));
  87         if use_beta_likes == 1
  88             tign_ground(~infire) = beta_vect(~infire).*tign_ground(~infire)+(1-beta_vect(~infire)).*tign_flat(~infire);
  89         else
  90             tign_ground(~infire) = beta*tign_ground(~infire)+(1-beta)*tign_flat(~infire);
  91         end
  92         %flatten the top
  93         t_mask = tign_ground > tmax;
  94         tign_ground(t_mask) = tmax;
  95         %tign_temp = imgaussfilt(tign_ground,smooth_ground );
  96         tign_temp = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_ground,a,b);
  97         tign_ground(~infire) = tign_temp(~infire);
  98 %         a = 1/10;%1/2-1/(2*i);
  99 %         tign_ground(infire) = a*tign_ground(infire)+(1-a)*tign_temp(infire);
 100 %         figure(73),scatter3(ps.red.fxlong(~infire),ps.red.fxlat(~infire),tign_ground(~infire))
 101 %         pause(1/2)
 102         figure(159)
 103         contourf(ps.red.fxlong,ps.red.fxlat,tign_ground,20,'k'),hold on
 104         scatter(pts(1:10:end,2),pts(1:10:end,1),'*r'),hold off
 105         t_str = sprintf('Perimeter Shrinking \n Iteration %d',i);
 106         save_str = sprintf('perim_shrink_%d',i);
 107         xlabel('Lon'),ylabel('Lat'),title(t_str)
 108 %         xl=[-121.5094 -121.2496];xlim(xl)
 109 %         yl = [46.0367   46.2035];ylim(yl)
 110 %         savefig(save_str);
 111 %         saveas(gcf,[save_str '.png']);
 112         %figure(160),mesh(ps.red.fxlong,ps.red.fxlat,tign_ground);
 113         %pause(.5)
 114         %t_min(i) = min(tign_ground(:));
 115     end
 116 end
 117 %plot(t_min)
 118 % figure,scatter(igrid(infire),jgrid(infire))
 119 % time_mask = grounds.land(:,5) > max(ps.points(:,3));
 120 % ground_mask = logical(time_mask.*~infire);
 121 tign_new = tign_ground;
 122 %%%plot and save original perimeter
 123 % figure(159),hold on
 124 % scatter(pts(:,2),pts(:,1),'*r')
 125 % contourf(ps.red.fxlong,ps.red.fxlat,ps.red.tign,20,'k'),hold off
 126 % t_str = sprintf('Perimeter Shrinking \n Iteration %d',0);
 127 % save_str = sprintf('perim_shrink_%d',0);
 128 % xlabel('Lon'),ylabel('Lat'),title(t_str)
 129 % xlim(xl);
 130 % ylim(yl);
 131 % savefig(save_str);
 132 % saveas(gcf,[save_str '.png']);
 133
 134 %%%%%%% end ground detection work %%%%%%
 135
 136
 137 idx = ps.idx;
 138 fig_num = 33;
 139 pts_length = length(ps.grid_pts);
 140 %max time to look at detection data
 141 max_l2_time = max(max(pts(:,3)));
 142
 143
 144
 145
 146 %tign_flat is constant fire arrival time
 147 %tign_flat=ps.red.end_datenum*ones(size(tign));
 148 title_str = 'Analysis';
 149
 150 %make matrix of forecast tign time for the points in the graph
 151 t_times = zeros(pts_length,1);
 152 for i = 1:pts_length
 153     t_times(i) = tign(idx(i,1),idx(i,2));
 154 end
 155
 156
 157
 158 %data likelikehoods
 159 time_diff = (pts(:,3)-t_times)*24*3600;
 160 likes = p_like_spline(time_diff);
 161 alpha_vect = exp(likes);
 162
 163 %mask = ones(size(tign));
 164 %change back to z = tign
 165 % if exist('ps.tn','var')
 166 %     figure(68),contour3(ps.red.fxlong,ps.red.fxlat,ps.tn-ps.red.start_datenum,20)
 167 % else
 168 %     figure(68),contour3(ps.red.fxlong,ps.red.fxlat,tign-ps.red.start_datenum,20)
 169 % end
 170 % hold on
 171 %lons = ps.red.fxlong;
 172 %lats = ps.red.fxlat;
 173
 174 %just squish the ignition point
 175 %bigger k ==> more smoothing
 176
 177 norms=[];
 178
 179 %random multiplier, increase for larger grids
 180 %perturbs points on path in x-y plane
 181 % try computing this as a fraction of grid size
 182 rm = round(m/100)+1;
 183
 184 % random multiplier, keep the same
 185 % perturbs points downward in time to
 186 rt = 1;
 187 % weight for tign_new
 188
 189 %alhpa blends  estimate of tign at a point with old estimate
 190 % new_estimate = alph*current_setimate + (1-alpha)*old_estimate
 191 % for data assimilation, alpha will be computed from exp(likelihood)
 192 % alpha = 0.5;
 193 %constant for smooth in rlx_shp
 194 alpha_2 = 0.7; %smaller alph_2 ==> smoother
 195 %number of loops to run
 196 smoothings = 20;
 197 for k = 1:smoothings
 198 %     figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new)
 199 %     title(title_str)
 200     %hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3),'*r'),hold off
 201     %pause(3/k)
 202     %plotting a slice
 203     %figure(fig_num+17),plot(tign_new(:,202));
 204     for i = 1:length(ps.paths)
 205         p = ps.paths(i).p;
 206         %         figure(73),hold on
 207         %         plot3(pts(p,2),pts(p,1),pts(p,3)-ps.red.start_datenum,'r')
 208         %         hold off
 209         %plot3(pts(p,2),pts(p,1),tign(idx(p,1),idx(p,2))-ps.red.start_datenum,'g')
 210         for j = 1:length(p) %p(j) is the current detection vertex
 211             tign_old = tign_new;
 212             %mesh indices for path points, perturbed
 213             p_i = idx(p(j),1);%+rm*round(randn);
 214             p_j = idx(p(j),2);%+rm*round(randn);
 215             %%% make mean of old and new, in small block around path point
 216             %alpha is now the data likelikehood
 217             %alpha = alpha_vect(p(j));
 218             alpha = 0.0;
 219             tign_new(p_i-round(rm*rand):p_i+round(rm*rand),p_j-round(rm*rand):p_j+round(rm*rand)) = alpha*tign_new(p_i,p_j) + (1-alpha)*pts(p(j),3)-rt;
 220             %%%% alternate strategy.
 221             %             if k == 1 && j > 1
 222             %             t1 = min(tign_new(p_i,p_j),tign_new(p_i-1,p_j-1));
 223             %             t2 = max(tign_new(p_i,p_j),tign_new(p_i-1,p_j-1));
 224             %             dt12 = t2-t1;
 225             %             dg = pts(p(j),3)-pts(p(j-1),3);
 226             %             time_shift = 0.5*(dg-dt12);
 227             %             t1 = t1-time_shift;
 228             %             t2 = t2+time_shift;
 229             %             tign_new(p_i,p_j) = max(t1,t2);
 230             %             tign_new(p_i-1,p_j-1) = min(t1,t2);
 231             %             end
 232
 233             %tign_flat(p_i-round(rm*rand):p_i+round(rm*rand),p_j-round(rm*rand):p_j+round(rm*rand)) = 0.5*(tign_new(p_i,p_j) + pts(p(j),3)-rt*rand);
 234
 235             %no need to interpolate - already done
 236             % interpolate new point between adjacent path detections
 237             if j > 1
 238                 %weighted average will move new point close to that with
 239                 %higher FRP
 240                 %                 frp1 = ps.points(p(j-1),5);
 241                 %                 frp2 = ps.points(p(j),5);
 242                 %                 w1 = frp1/(frp1+frp2);
 243                 %                 w2 = frp2/(frp1+frp2);
 244                 w1 = 0.5;
 245                 w2 = 0.5;
 246                 %assign tign for all in small, block around midpoint
 247                 %             frp1 = ps.points(p(j-1),5);
 248                 %             frp2 = ps.points(p(j),5);
 249                 %             w1 = frp1/(frp1+frp2);
 250                 %             w2 = frp2/(frp1+frp2);
 251                 new_lat = w1*pts(p(j-1),1)+w2*pts(p(j),1);
 252                 new_lon = w1*pts(p(j-1),2)+w2*pts(p(j),2);
 253                 new_t = w1*pts(p(j-1),3)+w2*pts(p(j),3);
 254                 [new_i,new_j,new_lat,new_lon]= fixpt(ps.red,[new_lat,new_lon]);
 255                 tign_new(new_i-round(rm*rand):new_i+round(rm*rand),new_j-round(rm*rand):new_j+round(rm*rand)) = new_t-rt*rand;
 256
 257                 %tign_new(new_i-round(rm*rand):new_i+round(rm*rand),new_j-round(rm*rand):new_j+round(rm*rand)) = new_t-rt*rand;
 258
 259                 %mask(new_i,new_j)=0;
 260             end
 261         end
 262     end
 263     %size of local averaging to apply aoutomate by grid size?
 264     %patch = max(1,round(sqrt(smoothings-k)));
 265     patch = 4*ceil(max(m,n)/100);
 266     fprintf('Using patch size %d in rlx_shp function \n',patch);
 267     %patch = 2;
 268
 269     %%%%% make use of ground detections %%%%
 270     %tign_new(~infire) = beta*tign_new(~infire)+(1-beta)*tign_flat(~infire);
 271
 272     %smooth the tign
 273     %tign_new(tign_new < t0) = t0;
 274     tign_new = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_new,a,b);
 275     if k < smoothings/5
 276         tign_new = imgaussfilt(tign_new,1);
 277         tign_new = rlx_shp(tign_new,alpha_2,patch);
 278     end
 279 %   tign_flat = rlx_shp(tign_flat,alpha_2,patch);
 280
 281     %collect information about tign at the detection points
 282     for i = 1:pts_length
 283         t_times(i) = tign_new(ps.idx(i,1),ps.idx(i,2));
 284         %flat_times(i) = tign_flat(ps.idx(i,1),ps.idx(i,2));
 285     end
 286     norms(k,1) = norm(t_times-ps.points(:,3),2);
 287     norm_tign_new = norms(k,1);
 288
 289
 290     %only do norm for times before final detection time
 291     time_mask = tign_new < pts(end,3);  %max(max(pts(:,3)));
 292     norms(k,2) = norm(tign_new(time_mask)-tign_old(time_mask));%/norm(tign_old(time_mask));
 293     figure(fig_num+3);plot(norms(1:k,1));
 294     tstr = sprintf('Norm of difference between successive \n TIGN after each interpolation');
 295     title(tstr)
 296     figure(fig_num+4);plot(1:k,norms(1:k,2))
 297     tstr = sprintf('Norm of difference between times of \n detections and TIGN at detectionon locations');
 298     title(tstr);
 299     xlabel('Iterations of Interpolation')
 300     %temp_var(k) = min(tign_new(:))-ps.red.start_datenum;
 301     %figure(fig_num+5);plot(1:k,temp_var(1:k)*24),title('Change in ignition time'),xlabel('iteration'),ylabel('hours')
 302     fprintf('Loop %d complete norm of diff = %f \n', k,norms(k))
 303     norm_break = 1;
 304     if k > 2 && norms(k,norm_break) >= norm(k-1,norm_break) && norms(k,norm_break) >= norms(k-2,norm_break)
 305         fprintf('graph norm increase \n')
 306         rm = 0;
 307         rt = 0;
 308         %break
 309     end
 310 end
 311 %tign_new = [];
 312 %smooth the result
 313
 314
 315 % %%% try using ground detections after paths....
 316 ground_steps = 4;
 317 tign_ground = tign_new;
 318 if gq
 319     data_area = sum(infire);
 320     for i = 1:ground_steps
 321         fire_mask = tign_ground < tmax-0.1;
 322         fire_area = sum(fire_mask(:));
 323         out_side = fire_mask(:)-infire;
 324         out_side(out_side<0)=0;
 325         out_sum(i) = sum(out_side(:));
 326         out_fraction = out_sum(i)/data_area;
 327         if i > 2 && out_sum(i) == out_sum(i-1)
 328             smooth_ground = 0.9*smooth_ground;
 329             fprintf('decreasing ground smoothing \n')
 330         end
 331         if out_fraction < 0.1
 332             %break
 333         end
 334         fprintf('Fire area: %f pixel_out/data area: %f  pixels_out %f\n',fire_area,out_sum(i)/data_area,out_sum(i));
 335         tign_ground(~infire) = beta*tign_ground(~infire)+(1-beta)*tign_flat(~infire);
 336         %tign_ground(~infire) = beta_vect(~infire).*tign_ground(~infire)+(1-beta_vect(~infire)).*tign_flat(~infire);
 337         t_mask = tign_ground > tmax;
 338         tign_ground(t_mask) = tmax;
 339         tign_temp = imgaussfilt(tign_ground,1/6);
 340         %tign_temp = imgaussfilt(tign_ground,smooth_ground );
 341         %tign_temp = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_ground,a,b);
 342         tign_ground(~infire) = tign_temp(~infire);
 343 %         a = 1/10;%1/2-1/(2*i);
 344 %         tign_ground(infire) = a*tign_ground(infire)+(1-a)*tign_temp(infire);
 345 %         figure(73),scatter3(ps.red.fxlong(~infire),ps.red.fxlat(~infire),tign_ground(~infire))
 346 %         pause(1/2)
 347         figure(159)
 348         contourf(ps.red.fxlong,ps.red.fxlat,tign_ground,20,'k'),hold on
 349         scatter(pts(1:5:end,2),pts(1:5:end,1),'*r'),hold off
 350         t_str = sprintf('Perimeter Shrinking \n Iteration %d',i);
 351         save_str = sprintf('perim_shrink_%d',i);
 352         xlabel('Lon'),ylabel('Lat'),title(t_str)
 353 %         xl=[-121.5094 -121.2496];xlim(xl)
 354 %         yl = [46.0367   46.2035];ylim(yl)
 355 %         savefig(save_str);
 356 %         saveas(gcf,[save_str '.png']);
 357         %figure(160),mesh(ps.red.fxlong,ps.red.fxlat,tign_ground);
 358         %pause(.5)
 359         %t_min(i) = min(tign_ground(:));
 360     end
 361 end
 362 tign_new = tign_ground;
 363 %% end of post smoothing
 364
 365 tign_new = smooth_up(ps.red.fxlong,ps.red.fxlat,tign_new,a,b);
 366
 367
 368
 369 % %fix the ripple on the top
 370 % r_mask = tign_new >=max(tign_new(:))-0.05;
 371 % tign_new(r_mask) = max(tign_new(:));
 372
 373
 374 figure(fig_num),mesh(ps.red.fxlong,ps.red.fxlat,tign_new-t0)
 375 title(title_str)
 376 hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3)-t0,'*r'),hold off
 377
 378 %plot the forecast for comparison
 379 figure,mesh(ps.red.fxlong,ps.red.fxlat,ps.red.tign-t0)
 380 title('Forecast')
 381 hold on,scatter3(ps.grid_pts(:,2),ps.grid_pts(:,1),ps.points(:,3)-t0,'*r')
 382
 383 filled_contours(ps,tign_new)
 384
 385 end
 386
 387