cycling/moe_2d.m

   1 function [mx,my,s] = moe_2d(w1,w2)
   2 %w1 is wrfout from forecast
   3 %w2 is wrfout with "ground truth" fire arrival time - optional if perimeter
   4 %              is used
   5 perim_use = 0;
   6 if ~exist('w2','var')
   7     [fire_name,save_name,prefix,perim] = fire_choice();
   8     perim_use = 1;
   9 end
  10 [m,n] = size(w1.fxlong);
  11 %perim_use = input_num('Compare to Perimeter? Yes = 1',1,1)
  12
  13 %new grid size
  14 %compute grid sizes
  15 if isfield(w1,'analysis')
  16     w1.tign_g = w1.analysis;
  17     fprintf('Using w.analysis for comparison\n')
  18 end
  19 resize_grid = input_num('Resize grids for comparison? 0 = no',1,1);
  20 if resize_grid
  21     grid_space = input_num('Grid spacing?',250,1);
  22     red = subset_domain(w1,1);
  23     E = wgs84Ellipsoid;
  24     dlon= distance(red.min_lat,red.min_lon,red.min_lat,red.max_lon,E);
  25     dlat= distance(red.min_lat,red.min_lon,red.max_lat,red.min_lon,E);
  26     new_n = round(dlon/grid_space);
  27     new_m = round(dlat/grid_space);
  28     new_red = subset_small(red,new_m,new_n);
  29     w1.fxlong = new_red.fxlong;
  30     w1.fxlat = new_red.fxlat;
  31     w1.tign_g = new_red.tign_g;
  32     [m,n] = size(w1.fxlong);
  33 end
  34 if perim_use == 1
  35
  36     % load data if it has been processes already for another comparison
  37     load_perim = input_num('Load perimeter data? Yes = 1',0,0);
  38     if load_perim == 1 && exist('perim_moe.mat','file')
  39         load perim_moe.mat
  40     else
  41         perim_points = input_num('Number of perim points to use?',200)
  42         perim_struct = perim2gran(perim_points,perim);
  43         for i = 1:length(perim_struct)
  44             fprintf('%d : %s %s\n',i,perim_struct(i).file,datestr(perim_struct(i).time))
  45         end
  46         fprintf('End time of estimate %s\n',datestr(new_red.end_datenum));
  47         perim_choice = input_num('Which perimeter to use ?',3);
  48         %interpolate the perimeter to the grid
  49         perim_points = fixpoints2grid(w1,[perim_struct(perim_choice).lat',perim_struct(perim_choice).lon']);
  50         perim_points = unique(perim_points,'rows');
  51         %make polygon comparison
  52         %[x,y] = meshgrid(1:n,1:m);
  53         %in_perim = inpolygon(y,x,perim_points(:,1),perim_points(:,2));
  54         %in_perim = inpolygon(y,x,x(in_perim),y(in_perim));
  55         %in_perim = inpolygon(y,x,x(in_perim),y(in_perim));
  56
  57         in_perim = inpolygon(w1.fxlong,w1.fxlat,perim_points(:,4),perim_points(:,3));
  58         in_perim = inpolygon(w1.fxlong,w1.fxlat,w1.fxlong(in_perim),w1.fxlat(in_perim));
  59         figure,scatter(w1.fxlong(in_perim),w1.fxlat(in_perim),'k');
  60         hold on,scatter(perim_points(:,4),perim_points(:,3)),hold off
  61         title('old')
  62         %new method
  63         p = make_poly(perim_points(:,4),perim_points(:,3),2);
  64         in_perim2 = inpolygon(w1.fxlong,w1.fxlat,p(:,1),p(:,2));
  65         in_perim2 = inpolygon(w1.fxlong,w1.fxlat,w1.fxlong(in_perim2),w1.fxlat(in_perim2));
  66         figure,scatter(w1.fxlong(in_perim2),w1.fxlat(in_perim2),'k');
  67         hold on,scatter(perim_points(:,4),perim_points(:,3)),hold off
  68         title('New')
  69         p_use = input_num('Use which perim? New = 1',1);
  70         if p_use == 1
  71             in_perim = in_perim2;
  72         end
  73
  74         fprintf('Choose default bounds 2\n')
  75         r = subset_domain(w1,1);
  76         tign_g = r.tign_g; %datenum2time(w1.tign_g,r);
  77         perim_time = perim_struct(perim_choice).time;
  78         t_g = datenum2time(perim_time,r);
  79         t_g = r.end_time;
  80         x = []; y = [];
  81         save perim_moe.mat perim_points in_perim x y m n perim_struct perim_choice r perim_time t_g tign_g
  82     end
  83
  84
  85     perim_flat = zeros(m,n);
  86     perim_flat(in_perim) = 1;
  87     perim_area = sum(perim_flat(:));
  88
  89 else
  90     t_g = min(max(w1.tign_g(:)),max(w2.tign_g(:)))-1*3600;
  91     tign_g = w1.tign_g;
  92     %area of ground truth
  93     in_perim = w2.tign_g<=t_g;
  94
  95     perim_flat = zeros(m,n);
  96     perim_flat(in_perim) = 1;
  97     perim_area = sum(perim_flat(:));
  98
  99
 100 end %if perim_use...
 101
 102 %start comparison with w.tign_g
 103
 104 w1_mask = tign_g < t_g-100;
 105 w1_flat = zeros(m,n);;
 106 w1_flat(w1_mask) = 1;
 107 w1_area = sum(w1_flat(:));
 108
 109 %find overlap
 110 overlap = perim_flat+w1_flat;
 111 overlap(overlap<1.5) = 0;
 112 overlap(overlap>0) = 1;
 113 overlap = logical(overlap);
 114
 115 false_neg = perim_flat-w1_flat;
 116 false_neg(false_neg < 0) = 0;
 117 false_neg = logical(false_neg);
 118 false_pos = w1_flat - perim_flat;
 119 false_pos(false_pos < 0) = 0;
 120 false_pos = logical(false_pos);
 121
 122 %MOE
 123 mx = 1 - sum(false_neg(:))/perim_area;
 124 my = 1 - sum(false_pos(:))/w1_area;
 125 %s --- sorenson index
 126 s = 2*sum(overlap(:))/(w1_area+perim_area);
 127
 128 plot_on = 1;
 129 if plot_on == 1
 130     figure(129)
 131     close(gcf)
 132     figure(129),scatter(w1.fxlong(false_neg),w1.fxlat(false_neg),'r*')
 133     hold on,scatter(w1.fxlong(false_pos),w1.fxlat(false_pos),'k*'),
 134     scatter(w1.fxlong(overlap),w1.fxlat(overlap),'g*')
 135     legend('False Negative','False Positive','Overlap')
 136     tstr = sprintf('MOE = (%f,%f) \n Sorenson = %f',mx,my,s);
 137     title(tstr);
 138 end
 139
 140 end