cycling/cone_compare.m

   1 function [r1,r2,adjr0,outer] = cone_compare(ps,tign2)
   2 %compares feature of fire arrival cones
   3 %ps = tign_try(w), tign2 = squish(ps)
   4
   5
   6 if isfield(ps.red,'red')
   7     upsize = input_num('Interpolate to original "red" grid? 1 = yes',1,1);
   8     if upsize
   9         red_orig = ps.red.red;
  10         %interpolate back to original size grid
  11         F = scatteredInterpolant(double(ps.red.fxlat(:)),double(ps.red.fxlong(:)),tign2(:));
  12         tign2 = F(red_orig.fxlat,red_orig.fxlong);
  13         %figure,contour(red_orig.fxlong,red_orig.fxlat,tign2,20,'k')
  14         %title('resized to original grid spacing')
  15         ps.red = red_orig;
  16     end
  17 end
  18
  19 lon = ps.red.fxlong;
  20 lat = ps.red.fxlat;
  21 %forecast
  22 tign = ps.red.tign;
  23 %blur the data for smoother gradients
  24 st = 1;
  25 % tign = imgaussfilt(tign,st);
  26 % tign2 = imgaussfilt(tign2,st);
  27 % use snmooth_up function to do smoothing
  28 % tign = smooth_up(lon,lat,tign);
  29 % tign2 = smooth_up(lon,lat,tign2);
  30
  31 %compare areas of the two and plot over time
  32 [ac,gs] = area_compare(ps,tign2);
  33 close all
  34
  35 %compute area of fires
  36 %t_end = min(max(tign(:)),max(tign2(:)))-0.1;
  37 a1 = ac.area_a(end);
  38 a2 = ac.area_b(end);
  39 outer.a1 = a1;
  40 outer.a2 = a2;
  41 adjr0 = sign(a2-a1)*1/10*sqrt(abs(a1-a2)/a2);
  42 %adjr0 = 1/10*sqrt(abs(a1-a2)/a2);
  43 %make the top flate for each
  44 % tign(tign>=t_end)=t_end;
  45 % tign2(tign2>=t_end)=t_end;
  46 fprintf('Forecast Area: %d Data area: %d adjr0: %f\n',a1,a2,adjr0);
  47 % figure,contour(lon,lat,tign,[t_end t_end],'k')
  48 % hold on,contour(lon,lat,tign2,[t_end t_end],'b')
  49 % legend('forecast','estimate')
  50 % t_str = sprintf('Perimeters \n Forecast Area = %d    Data Area = %d',a1,a2);
  51 % title(t_str)
  52 cull = input_num('Thin data set? [1]',1,1);
  53 if cull ~= 1
  54     lon = lon(1:cull:end,1:cull:end);
  55     lat = lat(1:cull:end,1:cull:end);
  56     tign = tign(1:cull:end,1:cull:end);
  57     tign2 = tign2(1:cull:end,1:cull:end);
  58 end
  59
  60 %compute gradient step sizes
  61 E = wgs84Ellipsoid;
  62 [n,m] = size(lon);
  63 hx = distance(lat(1,round(m/2)),lon(1,1),lat(1,round(m/2)),lon(end,end),E)/m;
  64 hy = distance(lat(1,1),lon(round(n/2),1),lat(end,end),lon(round(n/2),1),E)/n;
  65 aspect_ratio = [1 hy/hx 1];
  66
  67
  68 %gradient in fire arrival time, unit vector
  69 [dx1,dy1] = fire_gradients(lon,lat,tign,1);
  70 % [dx1,dy1] = gradient(tign);
  71 % dx = dx1./sqrt(dx1.^2+dy1.^2);
  72 % dy = dy1./sqrt(dx1.^2+dy1.^2);
  73 % dx1=dx;
  74 % dy1=dy;
  75 % clear dx dy
  76
  77
  78 % figure,contour(lon,lat,tign,20,'k');
  79 % daspect(aspect_ratio)
  80 % hold on
  81 % quiver(lon,lat,dx1,dy1)
  82 % hold off
  83 % title('Contour and gradient of estimate')
  84 % dx1=dx1/hx;dy1=dy1/hy;
  85
  86 %unit vector
  87 [dx2,dy2] = fire_gradients(lon,lat,tign2,1);
  88 % figure,contour(lon,lat,tign2,20,'k');
  89 % hold on
  90 % quiver(lon,lat,dx2,dy2)
  91 % hold off
  92 %[dx2,dy2] = gradient(tign2);
  93 % dx2=dx2/hx;dy2=dy2/hy;
  94 %compute the gradient of the terrain
  95 %elev = ps.red.fhgt;
  96 %[aspect,slope,ey,ex] = gradientm(lat,lon,elev,E);
  97 % figure,contour(lon,lat,elev,'k')
  98 % hold on
  99 % quiver(lon(1:5:end,1:5:end),lat(1:5:end,1:5:end),ex(1:5:end,1:5:end),ey(1:5:end,1:5:end))
 100
 101 %slopes of fire directions by directional derivatives
 102 % sl1 = ex.*dx1+ey.*dy1;
 103 % sl2 = ex.*dx2+ey.*dy2;
 104 % sl_diff = sl1-sl2;
 105 % figure,histogram(sl_diff)
 106 % t_str=sprintf('Difference in slopes normal to fire front \n Mean = %f', mean(sl_diff(~isnan(sl_diff))));
 107 % xlabel('Slope differences')
 108 % ylabel('Number')
 109 % title(t_str)
 110
 111 %mask for only the fire cone
 112 t_msk1 = tign<max(tign(:))-0.1;
 113 t_msk2 = tign2<max(tign2(:))-0.1;
 114 t_msk = logical(t_msk1.*t_msk2);
 115
 116 %measure angles
 117 % theta1 = atan2(dy1(t_msk),dx1(t_msk));
 118 % theta2 = atan2(dy2(t_msk),dx2(t_msk));
 119 %no mask
 120 theta1 = atan2(dy1,dx1);
 121 theta2 = atan2(dy2,dx2);
 122
 123 td = theta1-theta2;
 124 td(td>pi) = td(td>pi)-2*pi;
 125 td(td<-pi) = td(td<-pi)+2*pi;
 126 td_msk = ~isnan(td);
 127 outer.avg_t_diff = mean(td(td_msk));
 128 outer.std_t_diff = std(td(td_msk));
 129 b_msk = abs(td)<pi/6;
 130 figure,histogram(td)
 131 format short
 132 str = ['05','10','20','40'];
 133 tstr= sprintf('Angle difference in gradients \n Mean : %f Std deviation: %f \n Polygon Interpolation',mean(td(td_msk)),std(td(td_msk)));%,str(str_num:str_num+1));
 134 title(tstr)
 135 xlabel('Difference of angles (radians)')
 136 ylabel('Number')
 137 % save_str = sprintf('patch_angle_diffs_%s_p',str(str_num:str_num+1))
 138 % savefig(save_str);
 139 % saveas(gcf,[save_str '.png']);
 140
 141 % figure
 142 % quiver(lon(t_msk),lat(t_msk),dx1(t_msk),dy1(t_msk))
 143 % % quiver(lon,lat,dy1,dx1)
 144 % hold on
 145 % quiver(lon(t_msk),lat(t_msk),dx2(t_msk),dy2(t_msk))
 146 % % quiver(lon,lat,dx2,dy2)
 147 % title('Gradients in fire surfaces, unit vectors')
 148 % legend('Forecast','Interpolated')
 149
 150 %get vectors which are not unit vectors
 151 [du1,dv1] = fire_gradients(lon,lat,tign,0);
 152 [du2,dv2] = fire_gradients(lon,lat,tign2,0);
 153
 154 % figure
 155 % quiver(lon(t_msk),lat(t_msk),du1(t_msk),du1(t_msk))
 156 % hold on
 157 % title('Gradients in fire surfaces')
 158 % legend('Forecast','Interpolated')
 159
 160 %get magnitudes of these vectors for comparison
 161 % m1 = sqrt(du1(t_msk).^2+dv1(t_msk).^2);
 162 % m2 = sqrt(du2(t_msk).^2+dv2(t_msk).^2);
 163 m1 = sqrt(du1.^2+dv1.^2);
 164 m2 = sqrt(du2.^2+dv2.^2);
 165 mdiff = m1-m2;
 166 g_msk = abs(mdiff)>=0.0;
 167 % figure,histogram(mdiff(g_msk));
 168 % avg_mdiff = mean(mdiff(g_msk));
 169 % std_mdiff = std(mdiff(g_msk));
 170 % tstr = sprintf('Histogram of difference in vector magnitudes \n mean = %f std = %f',avg_mdiff,std_mdiff);
 171 % title(tstr);
 172
 173 %use the slope from gradientm function instead
 174 %rate of spread
 175 rx1 = 1./du1/(24*3600);
 176 ry1 = 1./dv1/(24*3600);
 177 rx2 = 1./du2/(24*3600);
 178 ry2 = 1./dv2/(24*3600);
 179 r1 = sqrt(rx1.^2+ry1.^2);
 180 r2 = sqrt(rx2.^2+ry2.^2);
 181 % cut_off = est_max(ps,r2)
 182 % cut_off = min(0.1,cut_off);
 183 cut_off = 2;
 184 b1 = r1<cut_off;b2 = r2 < cut_off;b_msk = logical(b1.*b2);
 185 % figure
 186 % quiver(lon(b_msk),lat(b_msk),rx1(b_msk),ry1(b_msk))
 187 % hold on
 188 % quiver(lon(b_msk),lat(b_msk),rx2(b_msk),ry2(b_msk))
 189 % t_str =sprintf('ROS vectors for ROS < %f',cut_off);
 190 % title(t_str);
 191
 192 r_diff = r1-r2;
 193 %r_diff = imgaussfilt(r_diff,1/2);
 194 r_msk = abs(r_diff)<5;
 195 avg_r_diff = mean(r_diff(b_msk));
 196 std_r_diff = std(r_diff(b_msk));
 197 outer.avg_r_diff = avg_r_diff;
 198 outer.std_r_diff = std_r_diff;
 199 figure,histogram(r_diff(b_msk));
 200 tstr= sprintf('Differences in ROS, forecast-estimate \n Mean = %f  Std Dev. = %f \n Polygon Interpolation  %s % of Data',avg_r_diff,std_r_diff);%,str(str_num:str_num+1));
 201 title(tstr)
 202 xlabel('ROS (m/s)')
 203 ylabel('Number')
 204 %save_str = sprintf('patch_ros_diffs_%s_p',str(str_num:str_num+1));
 205 %savefig(save_str);
 206 %saveas(gcf,[save_str '.png']);
 207 r_fast = r_diff>0;%(avg_r_diff+1*std_r_diff);
 208 r_slow = r_diff<0;%(avg_r_diff-1*std_r_diff);
 209 % figure,scatter(lon(r_fast),lat(r_fast),'*r');
 210 % hold on,scatter(lon(r_slow),lat(r_slow),'b');
 211 % title('Locations for fuel adjustment')
 212 % legend('Forecast too fast','Forecast too slow')
 213
 214 %regression on slope and ros differences
 215 % r_diff(abs(r_diff)>2) = NaN;
 216 % sl_diff(abs(r_diff)>2) = NaN;
 217 % figure,scatter(sl_diff(~isnan(r_diff)),r_diff(~isnan(r_diff)));
 218 % mdl = fitlm(sl_diff(:),r_diff(:));
 219
 220 % figure,histogram(ps.red.nfuel_cat(r_fast)),xticks(1:14)
 221 % title('Fuel Types Where Fire is Burning too Fast')
 222 % xlabel('Fuel Type'),ylabel('Number')
 223 % figure,histogram(ps.red.nfuel_cat(r_slow)),xticks(1:14)
 224 % title('Fuel Types Where Fire is Burning too Slow')
 225 % xlabel('Fuel Type'),ylabel('Number')
 226
 227 % figure
 228 % quiver(lon(t_msk),lat(t_msk),rx1(t_msk),ry1(t_msk))
 229 % hold on
 230 % quiver(lon(t_msk),lat(t_msk),rx2(t_msk),ry2(t_msk))
 231 % title('ROS vectors')
 232 % legend('Ground Truth','Interpolated')
 233 % figure,histogram(ps.red.nfuel_cat(t_msk)),xticks(1:14)
 234 % title('Fuel types')
 235
 236
 237 %find where the fire is burning too fast
 238
 239 % fast = mdiff < avg_mdiff-1/2*std_mdiff;
 240 % slow = mdiff > avg_mdiff+1/2*std_mdiff;
 241 %figure,scatter(lon(slow),lat(slow)),
 242 % figure,scatter(lon(fast),lat(fast),'*r')
 243 %legend('Forecast is Slow','Forecat is fast')
 244
 245 %cluster the domain by fuel type and gradient difference
 246 % [s_idx,s_c] = kmeans([lon(:),lat(:),mdiff(:)],2);
 247 % figure,scatter(lon(s_idx==1),lat(s_idx==1),'r')
 248 % hold on,scatter(lon(s_idx==2),lat(s_idx==2),'b')
 249 % legend('Forecast too fast','Forecast too slow')
 250
 251 %compute slope of terrain
 252
 253 %compare by fuel type
 254 % for i = 1:13
 255 %     fuel_mask = ps.red.nfuel_cat == i;
 256 %     msk = logical(fuel_mask.*r_msk);
 257 %     fuel_rate_diff =  r1(msk)-r2(msk);
 258 %     mean_frd = mean(fuel_rate_diff);
 259 %     std_frd = std(fuel_rate_diff);
 260 %     format short
 261 %     fprintf('Fuel type: %d ROS difference: %f Std. Deviation: %f \n',i,mean_frd,std_frd);
 262 % end
 263
 264 %close all
 265
 266 r1 = mean(r1(b_msk));
 267 r2 = mean(r2(b_msk));
 268
 269
 270 end