cycling/evaluate_sim.m

   1 function score = evaluate_sim( wrfout_s, wrfout_c, wrf_time_step )
   2 %score = evaluate_sim( wrfout, perim )
   3 % inputs:
   4 %   wrfout_s - wrfout from a simulation w/o cycling or spinup or ROS adjustment, string with file path
   5 %   wrfout_c - wrfout from a simulation with cycling or spinup or ROS adjustment, string with file path
   6 %   wrf_time_step - optional string with time step to be read
   7 % output:
   8 %   score - number which gives how good the simulation matches the perimter
   9 %   data
  10 % note: needs the function kml2mat to run
  11
  12 %to do
  13 %
  14 % scatter detections first
  15 % use alpha on detections?
  16 % use wrfout with earliest time comparable to perimeter time by searching
  17 % through the time steps of the wrfouts and comparing with the perimeter
  18 % time. input will then be directory and domain level instead of wrfout
  19 % name
  20
  21 close all
  22
  23 min_fire_confidence = 7;
  24
  25 %exp_num = input_num('Compare No ROS adjust vs ROS adjust = [1] Cycle vs No cycle = [2] Two diff cycles = [3]',1,0);
  26 exp_num = input_num('Compare old method vs new method = [1] Cycle vs No cycle = [2] Two diff cycles = [3]',1,0);
  27 if exp_num == 3
  28    cycle_s = input_num('Cycle number for wrfout_s?',0);
  29    cycle_c = input_num('Cycle number for wrfout_c?',1);
  30 end
  31 %read wrfout section
  32 if nargin > 2
  33     w_c =read_wrfout_tign(wrfout_c,wrf_time_step);
  34     red_c = subset_domain(w_c,1);
  35     clear('w_c')
  36     w_s =read_wrfout_tign(wrfout_s,wrf_time_step);
  37     red_s = subset_domain(w_s,1);
  38     clear('w_s');
  39 else
  40     w_c =read_wrfout_tign(wrfout_c);
  41     red_c = subset_domain(w_c);
  42     clear('w_c');
  43     w_s =read_wrfout_tign(wrfout_s);
  44     red_s = subset_domain(w_s);
  45     clear('w_s');
  46 end
  47
  48
  49
  50 %collect information about simulations
  51 max_tign = max(red_s.max_tign,red_c.max_tign);
  52 min_tign = min(red_s.min_tign,red_c.min_tign);
  53
  54 % need to plot detections too, this will be moved elsewhere
  55 sim = input_num('Patch [1] or Camp [2] or Cougar [3] ? ',1,1);
  56 if sim == 1
  57     det_prefix = '../TIFs/';
  58     perim = '../PERIMs/patch_perims/';
  59 elseif sim == 2
  60     det_prefix = '../campTIFs/';
  61     perim = '../PERIMs/camp_perims/';
  62 else
  63     det_prefix = '../cougarTIFs/';
  64     perim = '../PERIMs/cougar_perims/';
  65 end
  66
  67 det_list=sort_rsac_files(det_prefix);
  68 fig.fig_map=0;
  69 fig.fig_3d=0;
  70 fig.fig_interp=0;
  71
  72
  73 if perim(end) == 'l'
  74     %read kml section, *.kml file
  75     a = kml2struct(perim);
  76 else
  77     % or
  78     fprintf('Reading directory of shape files \n')
  79     a = shape2struct(perim);
  80 end
  81
  82 %find perimeters from perim file
  83 p_count = 0;
  84
  85 %%%%%%% start 13 for special point %%%%%%%%
  86 %n gives  size of grid to use
  87 n = 100;
  88 for i = 1:length(a)
  89     %i, a(i)
  90     if strcmp(a(i).Geometry,'Polygon')
  91         p_count = p_count + 1;
  92         %i,a(i)
  93
  94         % get perimeter time
  95         a(i).p_string = a(i).Name(end-12:end);
  96
  97 %        only for kml file
  98 %         if a(i).p_string(1) ~= '1'
  99 %             a(i).p_string(1) = '0';
 100 %         end
 101 %         formatIn = 'mm-dd-yyyy HHMM';
 102         formatIn = 'yyyymmdd HHMM';
 103         %perim times are local, need to convert to UTC
 104         zone_shift = 6;
 105         if strcmp(a(i).Name(1:2),'ca') | strcmp(a(i).Name(1:2),'wa')
 106             %a(i).p_string = a(i).Name(end-12:end);
 107             zone_shift = 8;
 108             %formatIn = 'yyyymmdd HHMM';
 109             %a(i).p_string = a(i).p_string(end-12:end)
 110         end
 111         a(i).p_time = datenum(a(i).p_string,formatIn)+zone_shift/24;
 112
 113         %set decimate to an  postive integer to use just a subset of points
 114         %  in perimeter
 115         fprintf('Perimeter %s: %d points in the perimeter \n',a(i).p_string,length(a(i).Lon));
 116         if length(a(i).Lat) > 500
 117             decimate = 10;
 118             lats = a(i).Lat(1:decimate:end);
 119             lons = a(i).Lon(1:decimate:end);
 120         else
 121             lats = a(i).Lat;
 122             lons = a(i).Lon;
 123         end
 124
 125         %create regularly spaced data
 126         dx = (a(i).BoundingBox(2,1)-a(i).BoundingBox(1,1))/n;
 127         dy = (a(i).BoundingBox(1,2)-a(i).BoundingBox(2,2))/n;
 128
 129         xa = linspace(a(i).BoundingBox(1,1),a(i).BoundingBox(2,1),n);
 130         ya = linspace(a(i).BoundingBox(2,2),a(i).BoundingBox(1,2),n);
 131
 132         %find data inside of perimeter
 133         [x,y] = meshgrid(xa,ya);
 134         x = x(:);
 135         y = y(:);
 136         [in,on] = inpolygon(x,y,lons,lats);
 137         fires = logical(in+on);
 138         data = reshape(fires,n,n);
 139         %make all high confidence fires
 140         data = uint8(9.0*data);
 141         a(i).data = data;
 142         geotransform = [ a(i).BoundingBox(1,1) dx 0  a(i).BoundingBox(2,2) 0 dy];
 143         a(i).geotransform = geotransform;
 144         %save the file for use in data assimilation
 145         %save(a(i).TIF_name,'data','geotransform');
 146         %plot results
 147         plot_on = 0;
 148         if plot_on
 149             figure
 150             hold on
 151             plot(lons,lats);
 152             scatter(x,y);
 153             scatter(x(fires),y(fires),'*','r');
 154             title(a(i).Name);
 155             %plot(lons,lats)
 156             hold off
 157             %figure,mesh(data)
 158         end %if plot_on
 159
 160         %store perimter structure
 161         p_struct(p_count).time = a(i).p_time;
 162         p_struct(p_count).lats = lats;
 163         p_struct(p_count).lons = lons;
 164         p_struct(p_count).file = replace(a(i).Name,' ','_');
 165         p_struct(p_count).Name = a(i).Name;
 166
 167     end
 168 end %for
 169
 170 %save perim_struct.mat a
 171 fprintf('There were %i perimeters found in the data set\n',p_count)
 172
 173 %sort the struct by time first --> last
 174 T = struct2table(p_struct);
 175 sort_table = sortrows(T, 'time');
 176 p_struct = table2struct(sort_table);
 177
 178
 179 %compare perimiters with the wrfout
 180 % need to filter out perimeters which are before first ignition time, score
 181 % would be zero
 182
 183
 184 for i = 1:p_count
 185     if p_struct(i).time <= max_tign && p_struct(i).time >= min_tign
 186         fprintf('Comparing wrfout files with the infrared perimeter \n')
 187
 188         sim_fires_s = red_s.tign <= p_struct(i).time;
 189         sim_fires_c = red_c.tign <= p_struct(i).time;
 190
 191         max_lat_s = max(red_s.fxlat(sim_fires_s));
 192         min_lat_s = min(red_s.fxlat(sim_fires_s));
 193         max_lat_c = max(red_c.fxlat(sim_fires_c));
 194         min_lat_c = min(red_c.fxlat(sim_fires_c));
 195         max_lon_s = max(red_s.fxlong(sim_fires_s));
 196         min_lon_s = min(red_s.fxlong(sim_fires_s));
 197         max_lon_c = max(red_c.fxlong(sim_fires_c));
 198         min_lon_c = min(red_c.fxlong(sim_fires_c));
 199
 200         %put both tign and perim on the same grid
 201         max_lat = max([max_lat_s,max_lat_c,max(p_struct(i).lats(:))])+0.01;
 202         min_lat = min([min_lat_s,min_lat_c,min(p_struct(i).lats(:))])-0.01;
 203         max_lon = max([max_lon_s,max_lon_c,max(p_struct(i).lons(:))])+0.01;
 204         min_lon = min([min_lon_s,min_lon_c,min(p_struct(i).lons(:))])-0.01;
 205
 206         mesh_dim = 200;
 207         x_new = linspace(min_lon,max_lon,mesh_dim);
 208         y_new = linspace(min_lat,max_lat,mesh_dim);
 209         %scatter(x_new,y_new)
 210         [x_grid,y_grid] = meshgrid(x_new,y_new);
 211         x_vect = x_grid(:);
 212         y_vect = y_grid(:);
 213         tign_grid_s = griddata(red_s.fxlong,red_s.fxlat,red_s.tign,x_grid,y_grid);
 214         tign_grid_c = griddata(red_c.fxlong,red_c.fxlat,red_c.tign,x_grid,y_grid);
 215         %figure,mesh(x_grid,y_grid,tign_grid)
 216
 217         %mask of fire areas
 218         sim_fires_s = tign_grid_s <= p_struct(i).time;
 219         sim_fires_c = tign_grid_c <= p_struct(i).time;
 220
 221
 222         %compute areas of simulations
 223         area_sim_s = sum(sim_fires_s(:));
 224         x_fires_sim_s = x_grid(sim_fires_s);
 225         y_fires_sim_s = y_grid(sim_fires_s);
 226         area_sim_c = sum(sim_fires_c(:));
 227         x_fires_sim_c = x_grid(sim_fires_c);
 228         y_fires_sim_c = y_grid(sim_fires_c);
 229         shrink = 0.95;
 230         %find boundaries
 231         sim_boundary_s = boundary(x_fires_sim_s,y_fires_sim_s,shrink);
 232         sim_x_s = x_fires_sim_s(sim_boundary_s);
 233         sim_y_s = y_fires_sim_s(sim_boundary_s);
 234         sim_boundary_c = boundary(x_fires_sim_c,y_fires_sim_c,shrink);
 235         sim_x_c = x_fires_sim_c(sim_boundary_c);
 236         sim_y_c = y_fires_sim_c(sim_boundary_c);
 237
 238         %mesh(sim_fires)
 239
 240         % put perim on the grid in use...
 241         [in_p,on_p] = inpolygon(x_grid(:),y_grid(:),p_struct(i).lons(:),p_struct(i).lats(:));
 242         perim_fires = logical(in_p+on_p);
 243         area_perim = sum(perim_fires(:));
 244         x_fires_perim = x_grid(perim_fires);
 245         y_fires_perim = y_grid(perim_fires);
 246         %shrink = 0.75; already declared
 247         perim_boundary = boundary(x_fires_perim,y_fires_perim,shrink);
 248         perim_x = x_fires_perim(perim_boundary);
 249         perim_y = y_fires_perim(perim_boundary);
 250
 251         perim_fires_m = reshape(perim_fires,mesh_dim,mesh_dim);
 252
 253         %calculate score
 254         diff_s = abs(perim_fires_m - sim_fires_s);
 255         diff_c = abs(perim_fires_m - sim_fires_c);
 256         %figure, mesh(diff);
 257         area_diff_s = sum(diff_s(:));
 258         area_diff_c = sum(diff_c(:));
 259
 260         score_c(i) = (area_sim_c+area_perim-area_diff_c)/(area_sim_c+area_perim);
 261         score_s(i) = (area_sim_s+area_perim-area_diff_s)/(area_sim_s+area_perim);
 262
 263         %collect detection data here
 264         if i == 1
 265             time_bounds(1) = min_tign;
 266         else
 267             time_bounds(1) = p_struct(i-1).time;
 268         end
 269         time_bounds(2) = p_struct(i).time; % will use perimeter time
 270         g_c = load_subset_detections(det_prefix,det_list,red_c,time_bounds,fig);
 271         %g_s = load_subset_detections(det_prefix,det_list,red_s,time_bounds,fig);
 272         % need to loop throught the g_* files and collect all detections
 273
 274         for j = 1:length(g_c)
 275             if j == 1
 276                 det_g_c = g_c(j).data >= min_fire_confidence;
 277                 scatter_lon = g_c(j).xlon(det_g_c);
 278                 scatter_lat = g_c(j).xlat(det_g_c);
 279             end
 280             det_g_c = g_c(j).data >= min_fire_confidence;
 281             lon_g_c = g_c(j).xlon(det_g_c);
 282             lat_g_c = g_c(j).xlat(det_g_c);
 283             %scatter_det = [scatter_det det_g_c(:)];
 284             if size(lon_g_c) > 100
 285                 %decimate = uint8(log(size(lon_g_c)));
 286                 decimate = 2;
 287             end
 288             % if more than 100 detections, plot just a subset
 289             if sum(det_g_c(:)) > 0
 290                 scatter_lon = [scatter_lon(:);lon_g_c(1:decimate:end)];
 291                 scatter_lat = [scatter_lat(:);lat_g_c(1:decimate:end)];
 292             end
 293         end
 294         %plot perims w/o cycling
 295         figure(2*i);
 296         hold on
 297         %scatter detections here
 298         alpha(0.5);
 299         scatter(scatter_lon,scatter_lat,'m','*');        alpha(1.0);
 300         plot(sim_x_s,sim_y_s,'b');
 301         plot(perim_x ,perim_y,'r');
 302         %xticks(-112.72:0.02:-112.58)
 303         if exp_num == 2
 304             title_str = ('Perimeter Observation and Forecast without cycling');
 305             legend({'Satellite Fire Detections','Forecast without cycling','Infrared perimeter'});
 306             save_str = [p_struct(i).file '_s']; %% _s is for single run
 307         elseif exp_num == 1
 308               title_str = ('Perimeter Observation and Forecast using old method');
 309               legend({'Satellite Fire Detections','Forecast from old method','Infrared perimeter'});
 310               save_str = [p_struct(i).file '_old'];
 311 %             title_str = ('Perimeter Observation and Forecast without spin-up');
 312 %             legend({'Satellite Fire Detections','Forecast without spin-up','Infrared perimeter'});
 313 %             save_str = [p_struct(i).file '_No_spinup'];
 314 %             title_str = ('Perimeter Observation and Forecast without ROS  adjust');
 315 %             legend({'Satellite Fire Detections','Forecast without ROS adjust','Infrared perimeter'});
 316 %             save_str = [p_struct(i).file '_No_ROS_adj'];
 317         else %spinup_compare == 3
 318             title_str = sprintf('Perimeter Observation and Forecast for Cycle %d ',cycle_s);
 319             %title_str = ('Perimeter Observation and Forecast without spin-up');
 320             legend({'Satellite Fire Detections','Forecast','Infrared perimeter'});
 321             save_str = sprintf('%s_cycle_%d',p_struct(i).file,cycle_s);
 322         end
 323
 324         %legend({'Forecast without cycling','Infrared perimeter'});
 325         xlabel('Lon')
 326         ylabel('Lat')
 327         xlim([min_lon max_lon])
 328         ylim([min_lat max_lat])
 329
 330         score_str = sprintf('Score = %f',score_s(i));
 331         title({title_str,p_struct(i).Name,score_str});
 332 %         if spinup_compare ~= 1
 333 %             save_str = [p_struct(i).file '_s'];
 334 %         else
 335 %             save_str = [p_struct(i).file '_No_spinup'];
 336 %         end
 337
 338         savefig(save_str);
 339         saveas(gcf,[save_str '.png']);
 340         ax(i,1) = gca;
 341         hold off
 342         %plot perims w/ cycling
 343         figure(2*i+1);
 344         hold on
 345         %scatter detections here
 346         alpha(0.5);
 347         scatter(scatter_lon,scatter_lat,'m','*');
 348         alpha(1.0);
 349         plot(sim_x_c,sim_y_c,'b');
 350         plot(perim_x ,perim_y,'r');
 351         %xticks(-112.72:0.02:-112.58)
 352         if exp_num == 2
 353             legend({'Satellite Fire Detections','Forecast with cycling','Infrared perimeter'});
 354             title_str = ('Perimeter Observation and Forecast with cycling');
 355             save_str = [p_struct(i).file '_c'];
 356         elseif exp_num == 1
 357               title_str = ('Perimeter Observation and Forecast using new method');
 358               legend({'Satellite Fire Detections','Forecast from new method','Infrared perimeter'});
 359               save_str = [p_struct(i).file '_new'];
 360 %             legend({'Satellite Fire Detections','Forecast using spin-up','Infrared perimeter'});
 361 %             title_str = ('Perimeter Observation and Forecast using spin-up');
 362 %             save_str = [p_struct(i).file '_with_spinup'];
 363 %             legend({'Satellite Fire Detections','Forecast using ROS adjust','Infrared perimeter'});
 364 %             title_str = ('Perimeter Observation and Forecast using ROS adjust');
 365 %             save_str = [p_struct(i).file '_with_ROS_adj'];
 366         else %spinup_compare == 3
 367             title_str = sprintf('Perimeter Observation and Forecast for Cycle %d ',cycle_c);
 368             %title_str = ('Perimeter Observation and Forecast without spin-up');
 369             legend({'Satellite Fire Detections','Forecast','Infrared perimeter'});
 370             save_str = sprintf('%s_cycle_%d',p_struct(i).file,cycle_c);
 371         end
 372
 373         %legend({'Forecast with cycling','Infrared perimeter'});
 374         xlabel('Lon')
 375         ylabel('Lat')
 376         xlim([min_lon max_lon])
 377         ylim([min_lat max_lat])
 378 %         if spinup_compare ~= 1
 379 %             save_str = [p_struct(i).file '_c'];
 380 %         else
 381 %             save_str = [p_struct(i).file '_with_spinup'];
 382 %         end
 383
 384         score_str = sprintf('Score = %f',score_c(i));
 385         title({title_str,p_struct(i).Name,score_str});
 386         savefig(save_str);
 387         saveas(gcf,[save_str '.png']);
 388         ax(i,2) = gca;
 389         hold off
 390
 391     else
 392         if p_struct(i).time < min_tign
 393             fprintf('Minimum tign after perimeter time, score is zero. Not plotting. \n');
 394             score_s(i) = 1e-16;
 395             score_c(i) = 1e-16;
 396             %             sim_x_s = [];
 397             %             sim_y_s = [];
 398             %             sim_x_c = [];
 399             %             sim_y_c = [];
 400         else
 401             fprintf('Perimiter time after simulation end. No score assigned. \n')
 402         end
 403     end
 404
 405
 406 end
 407
 408
 409 %do one plot with everything
 410 % [m,n] =size(ax)
 411 % fnew = figure;
 412 % for j = 1:m
 413 %     ax_copy(j,1) = copyobj(ax(j,1),fnew);
 414 %     ax_copy(j,2) = copyobj(ax(j,2),fnew);
 415 %     subplot(m,n,2*j-1,ax_copy(j,1))
 416 %     subplot(m,n,2*j,ax_copy(j,2))
 417 %
 418 %
 419 % end %for j
 420
 421 whos
 422 %m_s = score_s > 0;
 423 m_c = score_c > 0.02;
 424
 425 mean_score_s = mean(score_s(m_c));
 426 mean_score_c = mean(score_c(m_c));
 427
 428
 429 score.score_s = score_s;
 430 score.score_c = score_c;
 431 score.mean_score_s = mean(score_s(m_c));
 432 score.mean_score_c = mean(score_c(m_c));
 433 score.file = p_struct.file;
 434 end
 435