detection/load_subset_detections.m

   1 function g = load_subset_detections(prefix,p,red,time_bounds,fig)
   2 %
   3
   4 itime=find(p.time>=time_bounds(1) & p.time<=time_bounds(2));
   5 d=p.file(itime);      % files within the specified time
   6 t=p.time(itime);
   7 fprintf('Selected %i files in the given time bounds, from %i total.\n',...
   8     length(d),length(p.time))
   9
  10 k=0;
  11 for i=1:length(d),
  12
  13     file=d{i};
  14     fprintf('%s file %s \n',stime(t(i),red),file);
  15     %v=readmod14(prefix,file,'silent'); original line
  16     if file(end) == 't'
  17         v=readmod14(prefix,file);
  18         l2 = 0;
  19         %[v.lon,v.lat] = meshgrid(v.lon,v.lat);
  20
  21     % reading l2 data instead of tifs
  22     else
  23         % take only as matched pairs
  24         l2 = 1;
  25         if mod(i,2) == 1
  26             file2 = d{i+1};
  27             fprintf('send to readl2data , i = %d \n',i)
  28             v = readl2data(prefix,file,file2);
  29             % loop over unclassified data to filter
  30             % this will turn matrices to vectors
  31             unclass = [0,1,2,6];
  32             for uu = 1:length(unclass)
  33                mask = v.data(:) ~= unclass(uu);
  34                v.data = v.data(mask);
  35                v.lon = v.lon(mask);
  36                v.lat = v.lat(mask);
  37                %check for filter fail
  38                fail = sum(v.data(:)==unclass(uu));
  39                if fail
  40                    fprintf('Data corruption in %, class %d',file,unclass(uu));
  41                end
  42             end
  43         end
  44     end
  45     % select fire detection within the domain
  46     xj=find(v.lon > red.min_lon & v.lon < red.max_lon);
  47     xi=find(v.lat > red.min_lat & v.lat < red.max_lat);
  48     idx = intersect(xi,xj);
  49     ax=[red.min_lon red.max_lon red.min_lat red.max_lat];
  50     if isempty(xi) | isempty(xj)
  51         fprintf('outside of the domain\n');
  52     else
  53         x=[];
  54         if l2 == 0
  55             x.data=v.data(xi,xj);    % subset data
  56         else
  57             x.data=v.data(idx);
  58         end
  59         x.det(1)=sum(x.data(:)==3); % water
  60         x.det(2)=sum(x.data(:)==5);  % land
  61         x.det(3)=sum((x.data(:)==7)); % low confidence fire
  62         x.det(4)=sum((x.data(:)==8)); % medium confidence fire
  63         x.det(5)=sum((x.data(:)==9)); % high confidence fire
  64         if ~any(x.det)
  65             fprintf(' no data in the domain\n')
  66         else
  67             if l2 == 0
  68                 k=k+1;
  69             else
  70                 if mod(i,2) == 1
  71                     k = k+1;
  72                 end
  73             end
  74
  75             fprintf('water %i land %i fire low %i med %i high %i\n',x.det)
  76             x.axis=[red.min_lon,red.max_lon,red.min_lat,red.max_lat];
  77             x.file=v.file;
  78             x.time=v.time;
  79             if l2 == 0
  80                 x.lon=v.lon(xj);
  81                 x.lat=v.lat(xi);
  82             else
  83                 x.lon = v.lon(idx);
  84                 x.lat = v.lat(idx);
  85             end
  86             %interpolate to fire mesh, projection????
  87             if l2 == 0
  88                 [x.xlon,x.xlat]=meshgrid(x.lon,x.lat);
  89                 %use scatter interpolation instead, seems to make no
  90                 %difference but need to check
  91                 x.fxdata=griddata(double(x.xlon(:)),double(x.xlat(:)),double(x.data(:)),red.fxlong,red.fxlat,'nearest');
  92                 %x.fxdata=interp2(v.lon,v.lat,v.data,red.fxlong,red.fxlat,'nearest');
  93             else
  94                 %this next line breaks the plotting of detections
  95                 %downstream in fire_pixels3d.m
  96                 du = 0.005;
  97                 dv = 0.005;
  98                 u_space = red.min_lat:du:red.max_lat;
  99                 v_space = red.min_lon:dv:red.max_lon;
 100                 [x.xlon,x.xlat]=meshgrid(v_space,u_space);
 101                 %[x.xlon,x.xlat]=meshgrid(x.lon,x.lat);
 102                 data = griddata(double(x.lon(:)),double(x.lat(:)),double(x.data(:)),x.xlon,x.xlat,'nearest');
 103                 x.fxdata=griddata(double(x.lon(:)),double(x.lat(:)),double(x.data(:)),red.fxlong,red.fxlat,'nearest');
 104                 x.data = data;
 105                 x.lon = v_space;
 106                 x.lat = u_space;
 107             end
 108             if fig.fig_interp,  % map interpolated data to reduced domain
 109                 figure(fig.fig_interp)
 110                 cmap=cmapmod14;
 111                 c=reshape(cmap(x.fxdata+1,:),[size(x.fxdata),3]);
 112                 surf(red.fxlong,red.fxlat,zeros(size(red.fxlat)),c,'EdgeAlpha',0.2);
 113                 title(['Detection interpolated to fire mesh ',stime(x.time,red)])
 114             end
 115             g(k)=x;   % store the data structure
 116             if fig.fig_map,
 117                 figure(fig.fig_map);clf
 118                 showmod14(x)
 119                 hold on
 120                 contour(red.fxlong,red.fxlat,red.tign,[v.time v.time],'-k');
 121                 %fprintf('image time            %s\n',datestr(x.time));
 122                 if plot_also_wind,
 123                     if x.time >= ss.min_time && x.time <= ss.max_time,
 124                         step=interp1(ss.time,ss.num,x.time);
 125                         step0=floor(step);
 126                         if step0 < ss.steps,
 127                             step1 = step0+1;
 128                         else
 129                             step1=step0;
 130                             step0=step1-1;
 131                         end
 132                         w0=step1-step;
 133                         w1=step-step0;
 134                         uu=w0*ss.uh(:,:,step0)+w1*ss.uh(:,:,step1);
 135                         vv=w0*ss.vh(:,:,step0)+w1*ss.vh(:,:,step1);
 136                         fprintf('wind interpolated to %s from\n',datestr(x.time))
 137                         fprintf('step %i %s weight %8.3f\n',step0,datestr(ss.time(step0)),w0)
 138                         fprintf('step %i %s weight %8.3f\n',step1,datestr(ss.time(step1)),w1)
 139                         fprintf('wind interpolated to %s from\n',datestr(x.time))
 140                         sc=0.006;quiver(w.xlong,w.xlat,sc*uu,sc*vv,0);
 141                     end
 142                 end
 143                 hold off
 144                 axis(ax)
 145                 drawnow
 146                 % M(k)=getframe(gcf);
 147                 % print(fig.fig_map,'-dpng',['fig',v.timestr]);
 148             end
 149             if fig.fig_3d,
 150                 hold on; fire_pixels_3d(fig.fig_3d,x,red.base_time)
 151             end
 152         end
 153     end
 154 end
 155 if ~exist('g','var'),
 156     fprintf('No files found.')
 157     g=struct([]);
 158 end
 159 fprintf('%i detections selected\n',length(g))