library(foreach)

library(RcppOctave)

library(rgdal)

# final output will be written to data directory here:

    writeRaster(res2,file=toutfile,overwrite=T,datatype='INT2S',options=c("COMPRESS=LZW", "PREDICTOR=2"),NAvalue=-32768)

    ## remove temporary files

    rmr(d);rmr(d2);rmr(psi);rmr(res);rmr(res2)

    print(paste("Finished Temporary File: ",toutfile))

}

foreach( i=1:nrow(df2)) %dopar% {

    ifile=df2$path[i]

    outfile=paste(datadir,"/mcd09ctif/",df2$sensor[i],"_",df2$month[i],".tif",sep="")

    ## create VRT of first band of the full image 

    fvrt=sub("[.]tif",".vrt",ifile)

    system(paste("gdalbuildvrt -b 1 ",fvrt," ",ifile))

    system(paste("gdal_merge.py --config GDAL_CACHEMAX 10000 -init -32768 -n -32768 -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes  -o ",outfile," ",tfiles,sep="")) 

    writeLines(paste("Finished ",outfile))

}

datadir="/mnt/data2/projects/cloud/"

# Create combined (MOD+MYD) corrected mean CF

        f=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*[O|Y].*_",sprintf("%02d",i),"[.]tif$",sep=""),full=T)

        ## Define output and check if it already exists

        tmcd=paste(datadir,"/mcd09ctif/MCD09_",sprintf("%02d", i),"_uncompressed.tif",sep="")

        tmcd2=paste(datadir,"/mcd09ctif/MCD09_",sprintf("%02d", i),".tif",sep="")

        ## check if output already exists

        system(paste("gdalwarp -overwrite -r average -co COMPRESS=LZW -co ZLEVEL=9  ",ops," ",paste(f,collapse=" ")," ",tmcd))

        ## update metadata

        tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS MOD09GA and MYD09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",

            "TIFFTAG_DOCUMENTNAME='Collection 5 MCD09 Cloud Frequency'",

            paste("TIFFTAG_DATETIME='2013",sprintf("%02d", i),"15'",sep=""),

              "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")

        system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py ",tmcd," ",paste("-mo ",tags,sep="",collapse=" "),sep=""))

        # create final fixed image

        writeLines(paste("Finished month",i))

    }

    hd=c("<ColorInterp>Palette</ColorInterp>","<ColorTable>")

    ft="</ColorTable>"

    colR=colorRampPalette(c("#08306b","#0d57a1","#2878b8","#4997c9","#72b2d7","#a2cbe2","#c7dcef","#deebf7","#f7fbff"))

    ct=paste("<Entry c1=\"",cols$red,"\" c2=\"",cols$green,"\" c3=\"",cols$blue,"\" c4=\"255\"/>")

    cti=grep("ColorInterp",vrt)  # get index of current color table

    vrt2=c(vrt[1:(cti-1)],hd,ct,ft,vrt[(cti+1):length(vrt)])

    ## update missing data flag following http://lists.osgeo.org/pipermail/gdal-dev/2010-February/023541.html

    csi=grep("<ComplexSource>",vrt2)  # get index of current color table

    write.table(vrt2,file=outfilevrt,col.names=F,row.names=F,quote=F)

    #system(paste("gdal_translate  -of VRT -scale 0 10000 0 100  ",outfilevrt," ",outfilevrt2))

                                        #    system(paste("pkreclass  -i ",outfilevrt," ",paste("-mo ",tags,sep="",collapse=" ")," ",outfilevrt," ",outfile2))

    tags=c(paste("TIFFTAG_IMAGEDESCRIPTION='Monthly Cloud Frequency for 2000-2013 extracted from C5 MODIS M*D09GA PGE11 internal cloud mask algorithm (embedded in state_1km bit 10).",

        "TIFFTAG_DOCUMENTNAME='Collection 5 Cloud Frequency'",

        paste("TIFFTAG_DATETIME='2014'",sep=""),

        "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")

}

### Download files from google drive

## This only works if google-cli is installed and has already been authenticated 

download=T

##  Get list of available files

rerun=T  # set to true to recalculate all dates even if file already exists

## define month-sensors to process

i=1

#jobs=jobs[jobs$sensor=="MYD09",]

        s2=sub("GA","",s)

        ## Define output and check if it already exists

        ttif2=paste(datadir,"/mcd09tif/",s2,"_",sprintf("%02d", m),".tif",sep="")

        ## check if output already exists

        ## Merge to geotif in temporary directory

        ## 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",

        system(paste("gdalwarp -overwrite ",ops," ",tvrt," ",ttif1))

        ## Compress file and add metadata tags