#bias[abs(bias)<5]=0  # set all bias<0.5 to zero to avoit minor adjustment and added striping in data

    #bias=bias*10

#    dc=d-bias  # Make 'corrected' version of data

#    NAvalue(bias)=0  #set bias==0 to NA to facilate plotting, areas that are NA were not adjusted

#    res=stack(list(dc=dc,d=d,bias=bias)) #,EstP=EstP,EstD1=EstD1,EstD2=EstD2

    res=stack(list(bias=bias)) #,EstP=EstP,EstD1=EstD1,EstD2=EstD2

#    rm(dt,dc)

    rm(dt)

exts=extf(xmin=-180,xmax=180,ymin=-30,ymax=30,size=10,overlap=0)

#exts=extf(xmin=-10,xmax=20,ymin=0,ymax=30,size=10,overlap=0)

#exts=extf(xmin=-60,xmax=60,ymin=-30,ymax=30,size=10,overlap=0)

## loop over sensor-months to create full grid of corrected values

for( i in 1:length(df2)){

    outfile=paste("tmp/",sub(".tif","",basename(file)),"_",sprintf("%03d",ti),".tif",sep="")

    writeLines(paste("Starting: ",outfile," tile:",ti," ( out of ",nrow(exts),")"))

    ## acount for scale of data is 10000*CF

    d=d*.01

    res=vsnr(d,gabor=psi,alpha=2,p=2,epsilon=1,prec=5e-6,maxit=50,C=1)

    if(!is.null(outfile)) writeRaster(res,file=outfile,overwrite=T,datatype='INT2S',options=c("COMPRESS=LZW", "PREDICTOR=2"))#,NAflag=-127)

    print(paste("Finished ",outfile))

}

system("rm tmp/test2.tif; gdal_merge.py -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes 'tmp/test_[0-9]*[.]tif' -o tmp/test2.tif")

names(res)=c("dc","d","bias")

NAvalue(res)=-32768

tplot=F

if(tplot){

    gcol=colorRampPalette(c("black","white"))

    levelplot(stack(list(Original=res[["d"]],Corrected=res[["dc"]])),col.regions=gcol(100),cuts=99,maxpixels=1e6)

    levelplot(stack(list(Original=res[["d"]],Corrected=res[["dc"]])),col.regions=gcol(100),cuts=99,maxpixels=1e6,ylim=c(10,13),xlim=c(-7,-1))

}

// MCD09CF

///////////////////////////////////////////////////////////

// Generate a cloud frequency climatology from M*D09 data

var driveFolder="ee_mcd09cf";

var run="g2"

// limit overall date range  (only dates in this range will be included)

var datestart=new Date("2000-01-01")  // default time zone is UTC

var datestop=new Date("2014-02-28")

// specify which months (within the date range above) to process

var verbose=false       // print info about collections along the way (slows things down)

var drawmap=false       // add image to map

// define regions

var globe = '[[-180, -89.95], [-180, 89.95], [180, 89.95], [180, -89.95]]';  

var sahara = '[[-18, 0], [-18, 30], [15, 30], [15, 0]]';  

// choose region to export (if exportDrive==true)

var region = globe

/////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////

// Generate some filtering variables based on input above and date

var currentdate = new Date();

var date= currentdate.getFullYear()+''+("0" + (currentdate.getMonth() + 1)).slice(-2)+''+currentdate.getDate();

print(date)

// identify start and stop years

var yearstart=datestart.getUTCFullYear();

var yearstop=datestop.getUTCFullYear();

var nYears=years.length;

//function to combine terra and aqua and limit by date

function getMCD09(modcol,datestart,datestop){

    var getdates=ee.Filter.date(datestart,datestop);

    return(ee.ImageCollection(modcol).filter(getdates));

}

//function to subset a collection by year and month

function fyearmonth(collection,year,month){

  // define filters

    var getyear=ee.Filter.calendarRange(year,year,"year");

    var getmonth=ee.Filter.calendarRange(month,month,"month");

  // extract values

    var output=collection.filter(getmonth).filter(getyear);

    return(output);

}

///////////////////////////////////////////////////////////////////////////

// Loop through months and get cloud frequency for each month in year range

var monthmean=function(collection,month){

  //For this month, build array to hold means for every year

    var mod = new Array (nYears);

  // loop over years and and put the mean monthly CF in the array

    for (var i=0; i<nYears; i ++) {

	mod[i]= fyearmonth(collection,years[i],month). // filter by year-month

	map(getcf).                              // extract cloud frequency

	mean();                                  // take the mean

    }

    if(verbose){print('Processing '+nYears+' years of data for Month '+month)}

  // build an image collection for all years for this month using the array

    return(ee.ImageCollection(mod));

}

    var i, pattern = 0;

    for (i = start; i <= end; i++) {

	pattern += Math.pow(2, i);

        .bitwise_and(pattern)

        .right_shift(start);

// function to extract MOD09 internal cloud flag

var getcf=function(img){ 

    var img2=getQABits(img.select(['state_1km']),10, 10, 'cloud'); // extract cloud flag from bit 10

    return(img2.

           mask(img2.gte(0)).           // mask out NAs

           gte(1).                      // Convert to logical

           multiply(ee.Image(10000)))};   // Multiply by 10000 to facilate taking mean as an integer 0-10000

    var img2=img.select(['num_observations_1km']);

  return img2.

        mask(img2.gte(0)).           // mask out NAs

        select([0],['nObs']).        // rename to nObs

        multiply(ee.Image(10000))};    // multiply to 100 to facilate taking mean as an integer

    var img2=img.select(['num_observations_1km']);

  return img2.

        mask(img2.gte(0)).           // mask out NAs

        gte(1).                      // Convert to logical to get percentage of days with n>0

        multiply(ee.Image(10000)).      // multiply to 100 to facilate taking mean as an integer

        select([0],['pObs'])};        // rename to nObs

//////////////////////////////////////////////////////////////////////////////////////////////////////

//  Clip high-latitude nonsense

//  For an unknown reason there are pixels at high latitudes with nobs>0 in winter (dark) months.

//  These values were identified by viewing the monthly CF and nObs datasets and identifying the latitude

//  at which there was a stripe of fully missing data (nObs=0) followed by an erroneous region at higher latitudes

//  This returns a polygon for any month which can be used to clip the image

var monthbox=function(month) {

  // limits for each month (12 numbers correspond to jan-december limits)

    var ymin=[-90,-90,-90,-90,-70,-70,-70,-77,-77,-90,-90,-90];

    var ymax=[ 77, 90, 90, 90, 90, 90, 90, 90, 90, 90, 77, 69];

  // draw the polygon, use month-1 because month is 1:12 and array is indexed 0:11

    var ind=month-1;

    var box = ee.Geometry.Rectangle(-180,ymin[ind],180,ymax[ind]);

    return (box)};

////////////////////////////////////////////////////

////////////////////////////////////////////////////

////////////////////////////////////////////////////

////////////////////////////////////////////////////

//  Start processing the data

//var mcol='MOD09GA'

//var  tmonth=1

// loop over collections and months

var mcols = ['MOD09GA','MYD09GA'];

for (var c=0;c<mcols.length;c++) {

    var mcol=mcols[c];

    for (var tmonth = monthstart; tmonth <= monthstop; tmonth ++ ) {

	if(verbose){

	    print('Starting processing for:'+ mcol+' month '+tmonth);

	}

// build a combined M*D09GA collection limited by start and stop dates

	var MCD09all=getMCD09(mcol,datestart,datestop);

	var MCD09m=monthmean(MCD09all,tmonth);

	if(verbose){

	    print('MCD09m info:');

	    print(MCD09m.getInfo());

	}

// take overall mean and SD across all years for this month

	var MCD09_mean=MCD09m.mean();

	var MCD09_sd=MCD09m.reduce(ee.call("Reducer.sampleStdDev"));

///  Process MODIS observation frequency/proportion

// get number of obs and proportion of days with at least one obs

	var MCD09_pObs = MCD09all.filter(ee.Filter.calendarRange(tmonth,tmonth,"month")).map(getprop).mean();

	var MCD09_nObs = MCD09all.filter(ee.Filter.calendarRange(tmonth,tmonth,"month")).map(getcount).mean();

//////////////////////////////////////////////////////////////////////////////////////////////////////

//  Process monthly cloud trends

//var trend= MCD09m.formaTrend(null,2)

//if(verbose) print(trend.getInfo())

/////////////////////////////////////////////////////////////////////////////////////////////////////

// Build a single 8-bit image with all bands

// mask by nobs (only keep pixels where nobs > 0.75)

var MCD09=MCD09_mean.

            addBands(MCD09_sd).

            addBands(MCD09_nObs).

            addBands(MCD09_pObs).

            mask(MCD09_nObs.gte(25)).

            clip(monthbox(tmonth)).

            int16();

	if(exportDrive){

  //define export filename

	    var filename=date+'_'+run+'_'+yearstart+yearstop+'_'+mcol+'_'+tmonth;

	    if(verbose){  print('Exporting to: '+filename)}

	    exportImage(MCD09,filename,

			{'maxPixels':1000000000,

			 'driveFolder':driveFolder,

			 'crs': 'SR-ORG:6974', //4326

			 'scale': '926.625433055833',

			 'region': region

			});

	}

    } // close loop through months

} // close loop through collections

//  var palette="000000,00FF00,FF0000";

    var palette="000000,FFFFFF";

    addToMap(MCD09.select([0]),{min:0,max:10000,palette:palette}, "mean");

    addToMap(MCD09.select([1]),{min:0,max:5000,palette:palette}, "sd");

    addToMap(MCD09.select([2]),{min:0,max:15000,palette:palette}, "nObs");

    addToMap(MCD09.select([3]),{min:0,max:10000,palette:palette}, "pObs");

//  addToMap(trend.select('long-trend'), {min:-10, max:10, palette:palette}, 'trend');

}

### Download files from google drive

download=T

if(download) system(paste("google docs get 2014032* ",datadir,"/mcd09ee",sep=""))

        system(paste("gdalbuildvrt -b 1 -b 2 -srcnodata -32768 ",tvrt," ",paste(df$path[df$month==m&df$sensor==s],collapse=" ")))

        ## specify sourc projection because it gts it slightly wrong by default #-ot Int16 -dstnodata -32768

        ops=paste("-s_srs ",proj,"  -t_srs 'EPSG:4326' -multi -srcnodata -32768  -ot Int16 -dstnodata -32768 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",

            "-co BIGTIFF=YES --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -co COMPRESS=LZW -co PREDICTOR=2")

        ops2=paste("-ot Int16 -co COMPRESS=LZW -co PREDICTOR=2 -stats")

            "Band Descriptions: 1) Mean Monthly Cloud Frequency x 10000 2) Standard Deviation of Mean Monthly Cloud x 10000'"),