cycling/time_score.m

   1 function [ score,diff_var ] = time_score(perim_path,wrfout_path,wrfout_time_step)
   2 %function [ score ] = time_score(perim_kml,wrfout)
   3 % Function computes a score for a fire simulation based on fire arrival
   4 % times of the model compared with satellite observations
   5 % needs the kml2struct function to work
   6 % Inputs:
   7 %       perim - string, path to kml file with perim(s) or to a directory
   8 %          with kmls files or shape file, shape files should be downloaded
   9 %          from geomac as zip, then put into a common directory
  10 %       wrfout - string, path to a wrfout file for a simulation
  11 %       wrfout_time_step - optional, string with the timestep to be read
  12 % Outputs:
  13 %       score - average of the differences in fire arrival times for each
  14 %       pixel in the perimeter
  15
  16 %%%%% read the wrfout file and create interpolant
  17 if nargin > 2
  18     ts = wrfout_time_step;
  19     fprintf('Will read wrfout at %s \n',ts)
  20     w = read_wrfout_tign(wrfout_path,ts);
  21 else
  22     w = read_wrfout_tign(wrfout_path);
  23 end
  24
  25 %%% compute score using data likelihood
  26 like_score = input_num('Use likelihood approach? No = [0] Yes = [1]',0,1);
  27 close all
  28
  29 red = subset_domain(w,1);
  30 %% set cone_mask < to use only the cone, no flat part at top
  31 cone_mask = red.tign_g <= max(red.tign_g(:));
  32 flat_top = red.tign_g == max(red.tign_g(:));
  33 %Fw = scatteredInterpolant(w.fxlong(:),w.fxlat(:),w.tign_g(:));
  34 Fr = scatteredInterpolant(red.fxlong(cone_mask),red.fxlat(cone_mask),red.tign(cone_mask));
  35
  36 %%%% testing the interpolation object
  37 %zq = Fr(red.fxlong,red.fxlat);
  38 %figure,mesh(red.fxlong,red.fxlat,red.tign_g),title('org cone')
  39 %figure,mesh(red.fxlong,red.fxlat,zq),title('interpolated new cone')
  40
  41 %%% reading the perimeters
  42 if perim_path(end) == 'l'
  43     fprintf('Using a kml file for perimeters \n')
  44     %%%%%% convert kml file to a struct
  45     temp_struct = kml2struct(perim_path);
  46
  47
  48 elseif perim_path(end) == '/'
  49
  50     perim_dat = ['kml';'shp'];
  51     p_type = input_num('Type of perimeter file to use? (1) kml (2) shp', 2,1)
  52     fprintf('Reading %s files in  directory %s \n',perim_dat(p_type,:),perim_path)
  53     if p_type == 1
  54         d = dir([perim_path,'*.kml'])
  55     else
  56         p_struct = shape2struct(perim_path)
  57         temp_struct = p_struct;
  58     end
  59
  60
  61 end % if perim_path...
  62
  63
  64
  65 %%%%%% find the perimeters in the struct, make new struct, find times of
  66 %%%%%% perim in UTC
  67 perim_idx = [];
  68 %perim_time = [];
  69 perim_count = 0;
  70
  71 for i = 1:length(temp_struct)
  72     % find the perimters, which have the tag 'Polygon'
  73     if  strcmp(temp_struct(i).Geometry,'Polygon')
  74         perim_count = perim_count + 1;
  75
  76         perim_idx = [perim_idx,i];
  77
  78         if perim_path(end) == 'l'
  79             time_str = temp_struct(i).Name(end-13:end);
  80             perim_date(perim_count) = datetime(time_str,'InputFormat','MM-dd-yyyy HHmm');
  81         else
  82             time_str = temp_struct(i).Name(end-12:end);
  83             perim_date(perim_count) = datetime(time_str,'InputFormat','yyyyMMdd HHmm');
  84         end
  85
  86         %create new stuct with only perims
  87         perim_struct(perim_count).Lon = temp_struct(i).Lon;
  88         perim_struct(perim_count).Lat = temp_struct(i).Lat;
  89         perim_struct(perim_count).Name = temp_struct(i).Name;
  90         perim_struct(perim_count).time = datenum(perim_date(perim_count));
  91         %perim_time = datenum(perim_date);
  92         %perim_struct(perim_count) = temp_struct(i);
  93         %perim_struct(perim_count).time = 0;
  94         % look for nan/inf in data last entry seems to always be a nan....
  95         %inf_lat = ~isfinite(temp_struct(i).Lat);
  96         %sum(inf_lat)
  97     end %if
  98 end
  99
 100 %sort the perimeters
 101
 102 perim_time = datenum(perim_date);
 103 fprintf('There were %d perimetrs in the kml file \n',perim_count)
 104 [sorted_times,sort_idx] = sort(perim_date);
 105
 106 %% compute scores and plot the data
 107
 108 perim_scores = zeros(perim_count,1);
 109 perim_vars = zeros(perim_count,1);
 110 figure(perim_count+1)
 111 if size(red.fxlong) < 500
 112     mesh(red.fxlong,red.fxlat,(red.tign-red.start_datenum))
 113 else
 114     mesh(red.fxlong(1:10:end,1:10:end),red.fxlat(1:10:end,1:10:end),(red.tign(1:10:end,1:10:end)-red.start_datenum))
 115 end
 116 %legend('Forecast')
 117 title('Forecast and Perimeters');
 118 xlabel('Lon'),ylabel('Lat'),zlabel('Simulation Time [days]')
 119 for j = 1:perim_count
 120     i = sort_idx(j);
 121     %p_lon = temp_struct(perim_idx(i)).Lon(1:end-1);
 122     p_lon = perim_struct(i).Lon(1:end-1);
 123     %p_lat = temp_struct(perim_idx(i)).Lat(1:end-1);
 124     p_lat = perim_struct(i).Lat(1:end-1);
 125
 126     % only evaluate perims before the end of the simulation
 127     % scatter plot the perim data
 128     if perim_struct(i).time <= red.end_datenum - 0.2
 129         %% ??? perim_time(i)-max(red.max_tign)
 130         z = perim_struct(i).time*ones(size(p_lon));
 131         z_interp = Fr(p_lon,p_lat);
 132         %%% time difference is perimter time minus forecast time
 133         %%%   thus early forecast time gives a positive number and
 134         %%%   data likelihood would be higher
 135         diff = (z-z_interp)*24;
 136         diff = diff(~isnan(diff));
 137         % test only fire arrival times before simulation end
 138         cut_top = 0;
 139         if cut_top
 140             cone = diff < max(diff)-0.1;
 141             diff = diff(cone);
 142         end
 143
 144         %like_score = 0;
 145         if like_score > 0
 146             %make or load spline for data likelihood
 147             if  ~exist('p_spline','var')
 148                 if exist('p_spline.mat','file')
 149                     fprintf('Loading spline from mat file \n')
 150                     load p_spline.mat
 151                 else
 152                     fprintf('Making and saving spline \n')
 153                     [p_spline,~,~]= make_spline(72,400);
 154                     save p_spline.mat p_spline
 155                 end
 156             else
 157                 fprintf('Using spline in workspace \n')
 158             end %make spline
 159         end %like_score adjustments
 160
 161
 162         %%%% output results %%%%
 163         if like_score == 0
 164             rel_error = abs(diff(~isnan(diff)))./z(~isnan(diff));
 165             perim_scores(i) = mean(abs(diff));
 166             %perim_scores(i) = mean(diff);
 167             perim_vars(i) = var(diff);
 168             title_spec = sprintf('Histogram of errors %s',perim_struct(i).Name);
 169             figure(i),histogram(diff),title(title_spec)
 170             xlabel('Forecast difference from IR perimeter [hours]')
 171             ylabel('Number of perimeter points')
 172             figure(perim_count+1),hold on
 173             %scatter3(p_lon,p_lat,z_interp,'*')
 174             scatter3(p_lon,p_lat,(z-red.start_datenum),'.')
 175             hold off;
 176             fprintf('%s : Score %f \n', perim_struct(i).Name, perim_scores(i) );
 177             fprintf('   mean %f var %f \n',mean(diff),var(diff));
 178
 179         else  %%% using likelihood approach
 180             %using old like
 181             %stretch=[0.5,10,5,10];
 182             %dw = ones(size(diff));2
 183             %likes = like2(dw,diff,stretch);
 184             likes = p_spline(diff);
 185             %%exp(likes) give "probabilty"
 186             %likes = exp(likes);
 187             perim_scores(i) = mean(likes);
 188             perim_vars(i) = var(likes);
 189             title_spec = sprintf('Histogram of likelihoods %s',perim_struct(i).Name);
 190             figure(i),histogram(likes),title(title_spec)
 191             xlabel('Perimter point likelihood')
 192             ylabel('Number of perimeter points')
 193             figure(perim_count+1),hold on
 194             %scatter3(p_lon,p_lat,z_interp,'*')
 195             scatter3(p_lon,p_lat,(z-red.start_datenum),'.')
 196             hold off;
 197             fprintf('%s : Score %f \n', perim_struct(i).Name, perim_scores(i) );
 198             fprintf('   mean %f var %f \n',mean(likes),var(likes));
 199
 200         end
 201
 202     else
 203         fprintf('%s perimter after simulation end \n',perim_struct(i).Name);
 204     end
 205
 206
 207 end
 208 hold off
 209
 210 %% use all perims foir the final score ???
 211 % perim_struct.Name
 212 % score_mask = zeros(length(perim_struct),1);
 213 % for k = 1:length(perim_struct)
 214 %     query_string = sprintf('Use %s in score?',perim_struct(k).Name);
 215 %     score_mask(k) = input_num(query_string,1);
 216 % end
 217 % score_mask = logical(score_mask);
 218
 219 score_mask = logical([1 1 1 0 0]);
 220 score = perim_scores(score_mask);
 221 diff_var = perim_vars(score_mask);
 222 score = mean(perim_scores(perim_scores ~= 0));
 223
 224
 225
 226 end % function