cycling/growth_score.m

   1 function [ growth_struct] = growth_score( wrfout )
   2 %function [ gs ] = growth_score( wrfout,prefix )
   3 % inputs:
   4 %   wrfout - string, parth to wrfout file for evaluation
   5 %   ts - string, time step in wrfout file
   6 % outputs:
   7 %   growth_struct- struct with the following:
   8 %      gs - mean error in area gowth rate of forecast vs. area
   9 %       from detections
  10 %      data_time - vector with times of sat observations, (days from sim
  11 %       start)
  12 %      sat_area, fore_area - vectors with estimated areas of fires (grid nodes
  13 %       as units)
  14
  15 % make reduced structure
  16 % if wrfout(end) ~= 't'
  17 %     w = read_wrfout_tign(wrfout);
  18 % else
  19 %     fprintf('wrfout is a mat file, loading \n');
  20 %     load(wrfout);
  21 % end
  22
  23 fire_choice = input_num('which fire? Patch: [0], Camp: [1], Cougar: [3]',0,1)
  24 %[fire_name,save_name,prefix,perim] = fire_choice();
  25 cycle = input_num('Which cycle? ',0)
  26 if fire_choice == 1
  27     prefix='../campTIFs/';
  28     title_str = sprintf('Camp Fire -- Cycle %d',cycle)
  29     save_str = sprintf('camp_growth_c%d',cycle);
  30 elseif fire_choice == 0
  31     prefix='../TIFs/';
  32     title_str = sprintf('Patch Springs Fire -- Cycle %d',cycle)
  33     save_str = sprintf('patch_growth_c%d',cycle);
  34 else
  35     prefix = '../cougarTIFs/';
  36     title_str = sprintf('Cougar Creek Fire -- Cycle %d',cycle)
  37     save_str = sprintf('cougar_growth_c%d',cycle);
  38 end
  39
  40 w = read_wrfout_tign(wrfout);
  41 red = subset_domain(w,1);
  42 time_bounds = [red.start_datenum red.end_datenum];
  43
  44 % figures
  45 fig.fig_map=0;
  46 fig.fig_3d=0;
  47 fig.fig_interp=0;
  48
  49 %find detections
  50
  51 if exist('g_full.mat','file')
  52     load g_full.mat
  53 else
  54     p = sort_rsac_files(prefix);
  55     g = load_subset_detections(prefix,p,red,time_bounds,fig);
  56     save('g_full.mat', 'g', '-v7.3')
  57     %save g_full.mat g;
  58 end
  59
  60 if exist('growth.mat','file')
  61     fprintf('Growth structure file exists \n')
  62     load growth.mat
  63 end
  64
  65
  66
  67 %granules with active fire detections
  68 data_count = 0;
  69
  70 min_confidence = 7;
  71 for i=1:length(g)
  72     fprintf('Loading granule %d \n',i)
  73     %
  74     % add to detections already on fire grid
  75
  76     if i == 1
  77         fire_mask = g(i).data >=min_confidence;
  78         sat_lons = g(i).xlon(fire_mask);
  79         sat_lats = g(i).xlat(fire_mask);
  80     end
  81     fire_mask = g(i).data >= min_confidence;
  82     lon_update = g(i).xlon(fire_mask);
  83     lat_update = g(i).xlat(fire_mask);
  84     if sum(fire_mask(:)) > 0
  85         fprintf('Active fire in domain\n')
  86         data_count = data_count+1;
  87         sat_lons = double([sat_lons(:);lon_update(:)]);
  88         sat_lats = double([sat_lats(:);lat_update(:)]);
  89
  90         %% posibly add new granules
  91         if exist('temp_struct','var')&length(temp_struct.data_time) < i
  92             temp_struct.data_time(i) = 0;
  93             temp_struct.sat_area(i) = 0;
  94             temp_struct.fore_area(i) = 0;
  95             temp_struct.sat_rate(i) = 0;
  96             temp_struct.fore_rate(i) = 0;
  97             temp_struct.data_file{i} = 'temp';
  98         end %% adding new granule
  99
 100
 101         %% get satellite fire area
 102         if exist('temp_struct','var') & strcmp(g(i).file,temp_struct.data_file(data_count))
 103             fprintf('Temp struct has area already \n')
 104             sat_area(data_count) = temp_struct.sat_area(data_count);
 105             %temp_struct has data_time as days since sim start
 106             data_time(data_count) = temp_struct.data_time(data_count)+red.start_datenum;
 107             data_file{data_count} = temp_struct.data_file{data_count};
 108         else
 109             fprintf('Computing Satellite area in %s \n',g(i).file)
 110             %draw polygon around the fire grid detections
 111             sat_boundary = boundary(sat_lons,sat_lats);
 112             [sat_in,sat_on] = inpolygon(red.fxlong(:),red.fxlat(:),sat_lons,sat_lats);
 113             %figure,scatter(
 114             sat_area(data_count) = sum(sat_in) + sum(sat_on);
 115             % recor time and file name
 116             data_time(data_count) = g(i).time;
 117             data_file{data_count} = g(i).file;
 118         end
 119
 120         fprintf('Computing forecast area \n')
 121         fore_mask = red.tign < g(i).time;
 122         fore_area(data_count) = sum(fore_mask(:));
 123         %diff_fore_area(i) = fore_area(i) - fore_area(i-1)
 124
 125
 126         fprintf('Computing Growth Rates \n')
 127         if data_count == 1
 128             diff_sat_area(data_count) = sat_area(data_count);
 129             diff_fore_area(data_count) = fore_area(data_count);
 130         else
 131             diff_sat_area(data_count) = sat_area(data_count)-sat_area(data_count-1);
 132             diff_fore_area(data_count) = fore_area(data_count)-fore_area(data_count-1);
 133         end
 134         %else diff_sat_area(1) = temp_sat_area
 135
 136         fprintf('Computing relative error \n')
 137         %also plot the sequences of growth rates in same axes
 138         end
 139 end
 140
 141 %make areas monotonic increasing
 142 for i = 2:data_count
 143     sat_area(i) = max(sat_area(i),sat_area(i-1));
 144     fore_area(i) = max(fore_area(i),fore_area(i-1));
 145 end
 146
 147 %convert data_time to days since simulation
 148 data_time = (data_time - red.start_datenum);
 149 %bad data in patch fire towards end...
 150 if fire_choice == 0
 151     term = min(36,length(data_time));
 152     %term = length(data_time);
 153 else
 154     term = length(data_time);
 155 end
 156
 157 fprintf('%d granules had active fire data\n',data_count);
 158 figure,plot(data_time(1:term),sat_area(1:term))
 159 hold on, plot(data_time(1:term),fore_area(1:term))
 160 legend('Area within detection perimeter','Area within forecast perimeter')
 161 xlabel('Simulation Time [days]')
 162 ylabel('Simulation Area [grid nodes]')
 163
 164 title(title_str);
 165 savefig(save_str);
 166 saveas(gcf,[save_str '.png']);
 167 hold off
 168
 169 % figure,plot(data_time(1:term),diff_sat_area(1:term))
 170 % hold on, plot(data_time(1:term),diff_fore_area(1:term))
 171 % legend('Change in IR Perimter Area','Change in Forecast Area')
 172 % xlabel('Simulation Time [days]')
 173 % ylabel('Change in Area [grid nodes]')
 174 % hold off
 175 %gs = mean(abs(sat_area-fore_area));
 176 %relative error
 177 if fire_choice == 0
 178     gs = mean(abs(sat_area(1:term)-fore_area(1:term))./sat_area(1:term))
 179 else
 180     gs = mean(abs(sat_area-fore_area)./sat_area)
 181 end
 182 %compute rates
 183 sat_rate = sat_area;
 184 fore_rate = fore_area;
 185 for i = 2:length(sat_area)
 186     sat_rate(i) = sat_area(i) - sat_area(i-1);
 187     fore_rate(i) = fore_area(i) - fore_area(i-1);
 188 end
 189
 190
 191 growth_struct.gs = gs;
 192 growth_struct.cycle = cycle;
 193 growth_struct.data_time = data_time;
 194 growth_struct.sat_area = sat_area;
 195 growth_struct.fore_area = fore_area;
 196 growth_struct.sat_rate = sat_rate;
 197 growth_struct.fore_rate = fore_rate;
 198 growth_struct.data_file = data_file;
 199
 200 if ~exist('growth.mat','file')
 201     temp_struct = growth_struct;
 202     save growth.mat temp_struct;
 203 end
 204
 205
 206 end
 207