
% read MODIS (5-km) swath outputs
% The file name is designated by the start time of 5-minutes swath. 
% The file is in binary format. 
% The geo-location information and flux for every pixel within a 5-minutes swath 
% are stored in the order of latitude (270*406), longitude (270*406) 
% and radiative flux (270*406). 
% Usually, the output grid size is 270*406 pixels. 
% However, there are swaths that are of 270*408 pixels or 270*410 pixels.
% Those cases are denoted by the additional suffix "_408" or "_410" to
% the file names. 
import matlab.io.hdf4.*  
   org = 'nasa'
 
   modis = 'terra';% 'aqua';%
   year  = 2015;
   
   %domain = [-125,-100,30,55];
   domain = [-150,-60,20,60]; %[-60,60] may not be exactly. The first point center around [-150,20]
   gridStep = 0.05; gridStepStr='0.05deg';
   
   D2R = pi/180;
   R2D = 180/pi;
   earthR = 6371.2;
   gridR = 25;%10;%15;%5;%50;%(km)%footprint.
   MaxNumPixPerBox = fix(pi*(gridR/5+1)^2);

	 %-- station locations
   %
   stnCell=snowStns('st');

   nStns  = size(stnCell,1);
   stnID  = strcat(stnCell(:,3),'_',stnCell(:,5)); 
   stnLat = cell2mat(stnCell(:,2));
   stnLon = cell2mat(stnCell(:,1));
   %stnLon(stnLon<0) = stnLon(stnLon<0)+360;%0~360.
   stnLon(stnLon>180) = stnLon(stnLon>180)-360;%-180..180.
   stnFlux = nan(MaxNumPixPerBox,4);
   stFlx   = nan(1440/5,365,nStns);
   
 
   % dataPath=['/data/srb9/5kmMODIS/Results/'];
   dataPath=['/data/srb1/wchen/WORK/MODIS/SRB_SW/case_5km_cp_run/results/'];
   yyyy = num2str(year, '%04d');
   filepath =[dataPath,'/',yyyy,'/',modis,'/'];
   
  
   for doi=182:212
       doy   = num2str(doi, '%03d');
       fptn  = [filepath,'modis.',modis,'.',yyyy,doy,'.*.dat'];
       fList = dir(fptn);
       nFl   = length(fList);
   
      if nFl==0 
         disp('---read_MODIS_swath_new(): Empty directory!')       
      else
     
   for i=1:nFl
    [ds2] = read_MODIS_swath_v40(...
             modis, filepath, fList(i).name);
    imnt = (ds2.hor*12+ds2.mnt/5)+1;
    
    latR = ds2.latM.*D2R; %in radiance
    lonR = ds2.lonM.*D2R; %in radiance
      for istn = 1:nStns
   
         latS = stnLat(istn)*D2R; %; in radians
         lonS = stnLon(istn)*D2R; %; in radians

         dlat = abs(latS-latR);
         dlon = abs(lonS-lonR);

         %-- plane triangular distance
         %y = dlat.*earthR;
         %x = dlon.*earthR.*cos(latS);
         %dist = sqrt(x.^2+y.^2);
         %-- big circle distance.
         y = sqrt((cos(latR).*sin(dlon)).^2 + ...
                  (cos(latS).*sin(latR) - ...
                   sin(latS).*cos(latR).*cos(dlon)).^2);
         x = sin(latS).*sin(latR) + cos(latS).*cos(latR).*cos(dlon);
         dist = earthR*atan2(y,x);
         ind = find(dist<gridR);
         if ~isempty(ind)
             disp(['not empty',doy,num2str(ds2.hor),num2str(ds2.mnt)])
            dstIn = dist(ind);
            latIn = ds2.latM(ind);
            lonIn = ds2.lonM(ind);
            flxIn = ds2.varM(ind);

            [dstIn,ix] = sort(dstIn);
            latIn = latIn(ix);
            lonIn = lonIn(ix);
            flxIn = flxIn(ix);
         
            stnFlux(1:length(flxIn),:) = [dstIn,latIn,lonIn,flxIn];
            stFlx(imnt,doi,istn)=nanmean(flxIn);
                          

         end %if ~isempty(ind)
         
      end %for istn = 1:nStns
     end%nFL
      end %if nFL==0
	end % doi
    %keyboard
    save(['modis.',modis,'.4.st.',num2str(nStns),'sts.R',num2str(gridR,'%02d'),'.',num2str(year),'07.mat'],'stFlx')
function [outStns] = snowStns(tier)



if strcmp(tier, 'ARM')
   outStns = {-97.485 36.605 'ARM' 318 'SGP_C1'};


elseif strcmp(tier, 'SURFRAD')
   outStns = {-116.02 36.63 'SFR' 0    'dra';...
               -88.37 40.06 'SFR' 213  'bon';...
              -105.10 48.31 'SFR' 634  'fpk';...
               -89.52 34.37 'SFR' 98   'gwn';...
              -105.24 40.13 'SFR' 1689 'tbl';...
               -77.93 40.72 'SFR' 376  'psu';...
               -96.62 43.73 'SFR' 376  'sxf'};

elseif strcmp(tier, 'st')
   outStns = {-116.02 36.63 'SFR' 0    'dra';...
               -88.37 40.06 'SFR' 213  'bon';...
              -105.10 48.31 'SFR' 634  'fpk';...
               -89.52 34.37 'SFR' 98   'gwn';...
              -105.24 40.13 'SFR' 1689 'tbl';...
               -77.93 40.72 'SFR' 376  'psu';...
               -96.62 43.73 'SFR' 473 'sxf';...
               -97.49 36.61 'ARM' 214  'SGP'};
end


end

   %gridR = earthR*dLat*D2R/2;
