%  script to read in SSEBOP geotiff files and write out netcdf, one year at
%  a time.

for yydir=2000:2017;   %loop over years
    yystr=int2str(yydir);
    disp(yystr);
    dirname=['/Projects/CMIP5C/OBS_DATA/SSEBOP/MonthlyET/' yystr '/'];
    %dirname='./';
    a=dir([dirname 'm*.tif']);   % list all the months in that year-directory.  This can handle partial years by only processing existing files.
    %
    %  ****Must set the output type here**** GeoTIFFs were originally
    %  written in Uint16. (unsigned integer).   I have decided not to use
    %  uint16 for output because it can cause problems with other software.
    %  
    %  it turns out that compression of netCDF4 files creates adequately
    %  small files. 
    
    outtype='single';
    %outtype='int16';
    
    % get lats and lons from a sample file;  
    refgridfname='/Projects/CMIP5C/OBS_DATA/SSEBOP/MonthlyET/2000/m200001.modisSSEBopETactual.tif';  % since the grid changes on 2016-08, we establish a reference grid.
    [Zref,R0]=geotiffread(refgridfname);Zref=single(Zref);
    [mm,nn]=size(Zref);
    [ii,jj]=meshgrid([1:nn],[1:mm]);[lat,lon]=intrinsicToGeographic(R0,ii,jj);  %this creates a lat-lon grid from the geographic reference data. (I use grid centers as the lat-lon coordinates)
    latout=lat(:,1); lonout=lon(1,:)';
    clear ii jj lat lon yr mon dd
    Zout=zeros(nn,mm,length(a),outtype);
    
    for k = 1:length(a);
        fname=a(k).name;
        [Z,R]=geotiffread([dirname fname]);
        [mm,nn]=size(Z);
        if strcmp(outtype,'single')==1,
            Z =single(Z);
            mv=single(9999);  % missing value
        else
            mv=int16(9999);  % missing value
        end
        
        if mm==2834,   % note that files after 2016-08 have different dimension and no missing value set. 
            Z2=zeros(size(Zref),outtype);
            Z2(6:end-5 ,2:end-4)=Z;
            Z2(Zref==9999)=9999;   % set ocean pixels to 9999;
            Z=Z2;
        end
        Zout(:,:,k)=Z';    %output these in lon,lat,time so that it will be read (in C-notation by ncdump ) as time,lat,lon....
        
        yr(k) = str2num(fname(2:5));
        mon(k) = str2num(fname(6:7));
        dd(k)=1;   % set monthly value to 1st  of the month.
    end
    
    
    d0=datenum(2000,1,1);   % reference date
    d1=datenum(yr,mon,dd);
    dout=d1-d0;
    
    time=dout;
    fnameout=[fname(9:end-4) '.' fname(2:5) '.nc'];
    
    dimname = {'lon','lat','time'};
    dim = {lonout,latout,time};
    varname={'ETa'};
    var={Zout};
    vartypes={outtype};
    fillvalue={mv};
    vardesc={'Actual Evapotranspiration'};
    descglobal='SSEBOP estimate of actual evapotranspiration from MODIS satellite instrument';
    write_ndvar_netcdf_typed(fnameout,dimname,dim,var,varname,vardesc,vartypes,descglobal,fillvalue);
    ncwriteatt(fnameout,'ETa','units','mm');
    %ncwriteatt(fnameout,'ETa','_FillValue',mv);
    ncwriteatt(fnameout,'ETa','missing_Value',mv);
    ncwriteatt(fnameout,'time','units','days since 2000-1-1 0:0:0');    
end


