cycling/ros_adjust.m

   1 function adjustment = ros_adjust(forecast, analysis, time, fuel_map, p, w)
   2 % adjustment = ros_adjust(forecast, analysis, time)
   3 % function returns a rate of spread adjustment factor
   4 % inputs:
   5 %       forecast - tign for forecast
   6 %       analysis - tign for analysis
   7 %       time - time step for which comparison between forecast
   8 %           and analysis is to be made, probably use observation end time
   9 %           for the cycle in question
  10 %       fuel_map - matrix with fuel types, use w.nfuel_cat
  11 %       p  - struct with all the data about the simulation, eventually
  12 %       (possibly) reove the first three paramters of this function since
  13 %       they are in p anyway...
  14 %       w - contains the computational grid...
  15 %
  16 % output:
  17 %       adjustment - factor with which to multiply ROS adjustments in
  18 %           namelist.fire
  19
  20 fuel_type = 14;
  21 adj = ones(fuel_type,1);
  22
  23 threshold = 10; %minimum number of grid cells with given fuel type to
  24
  25 % use detections to estimate burn area, nomimal and high confidence only,
  26 % g is the struct with detection data created in detect_fit_lev2 with the
  27 % following :
  28 % g = load_subset_detections(prefix,p,red,time_bounds,fig);
  29 % move this inside so need to save data as .mat file is removed
  30
  31 %% first create a list of all detection pixels
  32 load g.mat
  33 %vectors for lon, lat of fire pixels
  34 fire_lon = [];
  35 fire_lat =[];
  36 for k = 1:length(g)
  37     %fprintf('Reading tif %i \n',k);
  38     mask = g(k).data >= 7;
  39     x = g(k).xlon(mask);
  40     fire_lon = [fire_lon;x(:)];
  41     y = g(k).xlat(mask);
  42     fire_lat = [fire_lat;y];
  43 end
  44 %% create and draw polygon around the detection pixels
  45 shrink = 0.5;
  46 fire_perim = boundary(fire_lon,fire_lat,shrink);
  47 x_perim = fire_lon(fire_perim);
  48 y_perim = fire_lat(fire_perim);
  49
  50 figure
  51 hold on
  52 plot(x_perim,y_perim);
  53 scatter(fire_lon,fire_lat);
  54 title('Detections and estimated perimeter');
  55 hold off
  56 %% fit this fire perimiter onto the computational grid
  57
  58 %find data inside of perimeter
  59 x_grid = w.fxlong(:);
  60 y_grid = w.fxlat(:);
  61 [in,on] = inpolygon(x_grid,y_grid,x_perim,y_perim);
  62 fires = logical(in+on);
  63 g_fires = sum(fires(:));
  64 w_fires = w.tign_g <= time;
  65 w_count = sum(w_fires(:));
  66 new_adjust = sqrt(g_fires/w_count);
  67
  68 %% computing new ROS to use by exponetial model where
  69 %% log(burn area) = 0.51x+2.539
  70
  71 adjrw=input_num('current adjrw : ',1);
  72 exp_adj = 0.51*log(g_fires/w_count)+adjrw;
  73
  74
  75
  76 %plotting some results
  77 figure
  78 hold on
  79 contour(w.fxlong,w.fxlat,w.tign_g)
  80 scatter(w.fxlong(fires),w.fxlat(fires),'*')
  81 title('Contours of forecast and estimate of burn are from detections')
  82 hold off
  83
  84 figure
  85 hold on
  86 plot(x_perim,y_perim);
  87 scatter(fire_lon,fire_lat);
  88 contour(w.fxlong,w.fxlat,w.tign_g)
  89 title('Detections and estimated perimeter with contours from forecast');
  90 hold off
  91
  92
  93
  94
  95 %old method
  96 a=analysis<=time;
  97 a_count = sum(a(:));
  98 f=forecast<=time;
  99 f_count = sum(f(:));
 100 adjustment = sqrt(a_count/f_count);
 101
 102
 103
 104 fprintf('Recomended single ROS adjustment factor by old method: %f \n',adjustment);
 105 fprintf('Recomended single ROS adjustment factor by new method: %f \n',new_adjust);
 106 fprintf('Recomended new adjrw by exponential method: %f \n',exp_adj);
 107
 108 %trying adjustments by fuel types.....
 109 for i = 1:fuel_type
 110
 111     %fprintf('compute ROS adjustemnt for fuel type %i\n',i);
 112     mask = fuel_map == i;
 113     a_new = a.*mask;
 114     a_new_count = sum(a_new(:));
 115     f_new = f.*mask;
 116     f_new_count = sum(f_new(:));
 117     if f_new_count >=10
 118         adj(i) = sqrt(a_new_count/f_new_count);
 119     end
 120     %fprintf('ROS adjusment factor for fuel type %i : %f \n\n',i,adj(i));
 121
 122
 123 end
 124
 125 fprintf('adjr0 = %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n', ...
 126     adj(1),adj(2),adj(3),adj(4),adj(5),adj(6),adj(7));
 127 fprintf('        %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n', ...
 128     adj(8),adj(9),adj(10),adj(11),adj(12),adj(13),adj(14));
 129
 130 fprintf('adjrw = %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n', ...
 131     adj(1),adj(2),adj(3),adj(4),adj(5),adj(6),adj(7));
 132 fprintf('        %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n', ...
 133     adj(8),adj(9),adj(10),adj(11),adj(12),adj(13),adj(14));
 134
 135 fprintf('adjrs = %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n', ...
 136     adj(1),adj(2),adj(3),adj(4),adj(5),adj(6),adj(7));
 137 fprintf('        %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n', ...
 138     adj(8),adj(9),adj(10),adj(11),adj(12),adj(13),adj(14));
 139
 140
 141
 142
 143