cycling/compare_perim.m

   1 function compare_perim(w,tign_new)
   2
   3 %w = read_wrfout_tign(f), f - wrfout file
   4 %tign_new = squish(w) - analysis or inital estimate to compare against perimeter
   5
   6
   7 %subset the fire domain to make it easier to work with
   8 red = subset_domain(w);
   9
  10
  11
  12 %% find which fire and data locations to use
  13 % sim = input_num('Patch [1] or Camp [2] or Cougar [3] ? ',3);
  14 % if sim == 1
  15 %     prefix = '../TIFs/';
  16 %     perim = '../PERIMs/patch_perims/';
  17 % elseif sim == 2
  18 %     prefix = '../campTIFs/';
  19 %     perim = '../PERIMs/camp_perims/';
  20 % else
  21 %     prefix = '../cougarTIFs/';
  22 %     perim = '../PERIMs/cougar_perims/';
  23 % end
  24
  25 [fire_name,save_name,prefix,perim] = fire_choice();
  26
  27 %% get perim data, arrange and store
  28 fprintf('Reading directory of shape files \n')
  29 a = shape2struct(perim);
  30
  31 %count perims, gather information
  32 %find perimeters from perim file
  33 p_count = 0;
  34
  35 %%%%%%% start 13 for special point %%%%%%%%
  36 %n gives  size of grid to use
  37 n = 100;
  38 for i = 1:length(a)
  39     %i, a(i)
  40     if strcmp(a(i).Geometry,'Polygon')
  41         p_count = p_count + 1;
  42         %i,a(i)
  43
  44         % get perimeter time
  45         a(i).p_string = a(i).Name(end-12:end);
  46         formatIn = 'yyyymmdd HHMM';
  47         %perim times are local, need to convert to UTC
  48         %patch fire
  49         zone_shift = 6;
  50         % cougar creek , camp fire
  51         if strcmp(a(i).Name(1:2),'ca') | strcmp(a(i).Name(1:2),'wa')
  52             %a(i).p_string = a(i).Name(end-12:end);
  53             zone_shift = 8;
  54             %formatIn = 'yyyymmdd HHMM';
  55             %a(i).p_string = a(i).p_string(end-12:end)
  56         end
  57         %datenum format as used by TIGN
  58         a(i).p_time = datenum(a(i).p_string,formatIn)+zone_shift/24;
  59
  60         %set decimate to an  postive integer to use just a subset of points
  61         %  in perimeter
  62         fprintf('Perimeter %s: %d points in the perimeter \n',a(i).p_string,length(a(i).Lon));
  63         if length(a(i).Lat) > 20000
  64             decimate = 8;
  65             lats = a(i).Lat(1:decimate:end);
  66             lons = a(i).Lon(1:decimate:end);
  67         else
  68             lats = a(i).Lat;
  69             lons = a(i).Lon;
  70         end
  71
  72         %create regularly spaced data
  73 %         dx = (a(i).BoundingBox(2,1)-a(i).BoundingBox(1,1))/n;
  74 %         dy = (a(i).BoundingBox(1,2)-a(i).BoundingBox(2,2))/n;
  75 %
  76 %         xa = linspace(a(i).BoundingBox(1,1),a(i).BoundingBox(2,1),n);
  77 %         ya = linspace(a(i).BoundingBox(2,2),a(i).BoundingBox(1,2),n);
  78
  79         %find data inside of perimeter
  80 %         [x,y] = meshgrid(xa,ya);
  81 %         x = x(:);
  82 %         y = y(:);
  83 %         [in,on] = inpolygon(x,y,lons,lats);
  84 %         fires = logical(in+on);
  85 %         data = reshape(fires,n,n);
  86         %make all high confidence fires
  87 %         data = uint8(9.0*data);
  88 %         a(i).data = data;
  89 %         geotransform = [ a(i).BoundingBox(1,1) dx 0  a(i).BoundingBox(2,2) 0 dy];
  90 %         a(i).geotransform = geotransform;
  91 %         %save the file for use in data assimilation
  92 %         %save(a(i).TIF_name,'data','geotransform');
  93 %         %plot results
  94         plot_on = 1;
  95         if plot_on && a(i).p_time < max(tign_new(:));
  96             figure
  97             contour(red.fxlong,red.fxlat,tign_new,[a(i).p_time a(i).p_time],'b');
  98             fire_mask = tign_new <= a(i).p_time;
  99             x_0 = min(red.fxlong(fire_mask));x_1=max(red.fxlong(fire_mask));
 100             y_0 = min(red.fxlat(fire_mask));y_1=max(red.fxlat(fire_mask));
 101             hold on
 102             plot(lons,lats,'r');
 103             xl_0 = min(lons);xl_1 = max(lons);
 104             yl_0 = min(lats);yl_1 = max(lats);
 105             title(a(i).Name);xlabel('Lon'),ylabel('Lat')
 106             xlim([min(xl_0,x_0) max(xl_1,x_1)]),ylim([min(y_0,yl_0) max(yl_1,y_1)])
 107             if a(i).p_time > red.end_datenum
 108                 legend('Interpolation forecast','Infrared Perimeter')
 109             else
 110             legend('Interpolation','Infrared Perimeter')
 111             end
 112             hold off
 113             %figure,mesh(data)
 114
 115         end %if plot_on
 116
 117         %store perimter structure
 118         p_struct(p_count).time = a(i).p_time;
 119         p_struct(p_count).lats = lats;
 120         p_struct(p_count).lons = lons;
 121         p_struct(p_count).file = replace(a(i).Name,' ','_');
 122         p_struct(p_count).Name = a(i).Name;
 123
 124     end
 125 end %for
 126
 127 %save perim_struct.mat a
 128 fprintf('There were %i perimeters found in the data set\n',p_count)
 129
 130 %sort the struct by time first --> last
 131 T = struct2table(p_struct);
 132 sort_table = sortrows(T, 'time');
 133 p_struct = table2struct(sort_table);
 134
 135
 136 end