library(mgcv)

registerDoMC(2)

## start raster cluster

#beginCluster(10)

tmpfs="tmp/"#tempdir()

ramdisk=F

rasterOptions(tmpdir=tmpfs,overwrite=T, format="GTiff",maxmemory=1e9)

jobs=expand.grid(month=1:12,sensor=c("MOD09","MYD09"))

## Loop over existing months to build composite netcdf files

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

        ## get month

        m=jobs$month[i]

        date=df$date[df$month==m][1]

        ## get sensor

        s=jobs$sensor[i]

        tvrt=paste(tmpfs,"/",s,"_",m,".vrt",sep="")

        ttif1=paste(tmpfs,"/",s,"_",m,".tif",sep="")

        ttif2=paste(tmpfs,"/",s,"_",m,"_wgs84.tif",sep="")

        ncfile=paste("data/mcd09nc/",s,"_",m,".nc",sep="")

        ## build VRT to merge tiles

        system(paste("gdalbuildvrt ",tvrt," ",paste(df$path[df$month==m&df$sensor==s],collapse=" ")))

        ## 

        mod=stack(tvrt)

        names(mod)=c("cf","cfsd","nobs","pobs")

        #################################

        ## correct for orbital artifacts

        ## sample points to fit correction model

        reg=extent(-20016036,20016036, -4295789.46149, 4134293.61707)  #banding region

        cmod=crop(mod,reg)

        names(cmod)=c("cf","cfsd","nobs","pobs")

        modpts=sampleRandom(cmod, size=10000, na.rm=TRUE, xy=T, sp=T)

        modpts=modpts[modpts$nobs>0,]  #drop data from missing tiles

        ### fit popbs correctionmodel (accounting for spatial variation in cf)

        modlm1=bam(cf~s(x,y)+pobs,data=modpts@data)

        summary(modlm1)

        modbeta1=coef(modlm1)["pobs"]

        #modlm2=bam(cf~s(x,y),data=modpts@data)

        #resid=residuals(modlm2)#,newdata=data.frame(x=modpts$x,y=modpts$y,pobs=100))

        #pred=predict(modlm1,newdata=data.frame(x=modpts$x,y=modpts$y,pobs=100))

        #plot(resid~pobs,data=modpts@data);abline(lm(resid~modpts$pobs))

        #plot(modlm1);points(modpts)

        #writeLines(paste(date,"       slope:",round(modbeta1,4)))

        ## mask no data regions (with less than 1 observation per day within that month)

        ## use model above to correct for orbital artifacts

        biasf=function(cf,cfsd,nobs,pobs) {

            ## drop data in areass with nobs<1

            drop=nobs<0|pobs<=50

            cf[drop]=NA

            cfsd[drop]=NA

            nobs[drop]=NA

            pobs[drop]=NA

            bias=round((100-pobs)*modbeta1)

            cfc=cf+bias

            return(c(cf=cf,cfsd=cfsd,bias=bias,cfc=cfc,nobs=nobs,pobs=pobs))}

#        treg=extent(c(-1434564.00523,1784892.95369, 564861.173869, 1880991.3772))

#        treg=extent(c(-2187230.72881, 4017838.07688,  -339907.592509, 3589340.63805))  #all sahara

        mod2=overlay(cmod,fun=biasf,unstack=TRUE,filename=ttif1,format="GTiff",

            dataType="INT1U",overwrite=T,NAflag=255, options=c("COMPRESS=LZW", "BIGTIFF=YES"))

#        system(paste("pksetmask -i ",modvrt," -m ",nObs,

#                     " --operator='<' --msknodata 1 --nodata 255 --operator='>' --msknodata 10 --nodata 10 -o ",modtif1,sep=""))

        ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata 255 -dstnodata 255 -r bilinear -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333",

            "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 500 -wm 500 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5")

#        ops=paste("-t_srs 'EPSG:4326' -srcnodata 255 -dstnodata 255 -multi -r bilinear -tr 0.008333333333333 -0.008333333333333",

#            "-co BIGTIFF=YES -ot Byte --config GDAL_CACHEMAX 300000 -wm 300000 -wo NUM_THREADS:10 -wo SOURCE_EXTRA=5")

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

        ##  convert to netcdf, subset to mean/sd bands

        trans_ops=paste(" -co COMPRESS=DEFLATE -a_nodata 255 -stats -co FORMAT=NC4 -co ZLEVEL=9 -b 4 -b 2")

        system(paste("gdal_translate -of netCDF ",trans_ops," ",ttif2," ",ncfile))

      time:units = \"days since 2000-01-01 ",ifelse(s=="MOD09","10:30:00","13:30:00"),"\" ;

        system(paste("ncks --fl_fmt=netcdf4 -L 9 -A ",tempdir(),"/",date,"_time.nc ",ncfile,sep=""))

        ## build correction factor explanation

                     " -a valid_range,CF,o,ub,\"0,100\" ",

                     " -a _FillValue,CF,o,ub,\"255\" ",

                     " -a missing_value,CF,o,ub,\"255\" ",

                     " -a correction_factor_description,CF,o,c,\"To account for variable observation frequency, CF in each pixel was adjusted by the proportion of days with at least one MODIS observation\" ",

                     " -a correction_factor,CF,o,f,\"",round(modbeta1,4),"\" ",

                     " -a valid_range,CFsd,o,ub,\"0,200\" ",

                     " -a _FillValue,CFsd,o,ub,\"255\" ",

                     " -a missing_value,CFsd,o,ub,\"255\" ",

d=stack("data/mcd09/2014113_combined_20002013_MOD09_1-img-0000000000-0000000000.tif")

#cd=coordinates(d)

#cd$x=as.integer(cd[,"x"]); cd$y=as.integer(cd[,"y"])

#cdr=rasterFromXYZ(cbind(cd,cd))

#names(cdr)

#d=stack(d,cdr)

pts=sampleRandom(d, size=100000, na.rm=TRUE, xy=T, ext=NULL, sp=T)

lm1=bam(cf~s(x,y)+pobs+nobs,data=pts@data)

beta1=coef(lm1)["pobs"]; beta1

beta2=coef(lm1)["nobs"]; beta2

#d2=d

#d2[["pobs"]]=100

#pred1=raster::predict(d,lm1,file=paste("data/bias/pred1.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]]))

#pred2=raster::predict(d2,lm1,file=paste("data/bias/pred2.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)#lm1=bam(as.vector(d[["cf"]]) ~ as.vector(d[["pobs"]]))

#pred3=overlay(pred1,pred2,function(x,y) x-y,file="data/bias/pred3.tif",overwrite=T)

#biasc=overlay(d[["cf"]], pred3,fun=sum,file="data/bias/biasc.tif",overwrite=T)

getbias=function(cf,cfsd,nobs,pobs) return((100-pobs)*beta1+(4-nobs)*beta2)

getbias(cf=c(50,40,30),c(1,1,1),c(3,4,5),c(82,87,87))

bias=overlay(d,fun=getbias, unstack=TRUE,file=paste("data/bias/bias.tif",sep=""),format="GTiff",dataType="INT1U",overwrite=T,NAflag=255)

dc=overlay(d[["cf"]],bias,fun=function(x,y) x+y,

    file=paste("data/bias/biasc.tif",sep=""),

    format="GTiff",dataType="INT1U",overwrite=T,NAflag=255) #,options=c("COMPRESS=LZW","ZLEVEL=9"))