registerDoMC(7)

setwd("~/acrobates/adamw/projects/cloud")

mpath="/home/adamw/acrobates/adamw/projects/environmental-layers/climate/research/cloud/MOD09/vsnr/"

    dt=.CallOctave("VSNR",d2,epsilon,p,as.matrix(gabor),alpha,maxit+1,prec,C1,argout ="u",verbose=T);

df2=list.files(paste(datadir,"/mcd09tif",sep=""),full=T,pattern="[0-9][.]tif$")

### Build the tiles

modexts=c(xmin=130,xmax=180,ymin=-50,ymax=0)

mydexts=c(xmin=-170,xmax=-140,ymin=-30,ymax=55)

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

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

    file=df2[i]

    outfile=paste(datadir,"/mcd09ctif/",basename(file),sep="")

#    outfile2=paste("data/mcd09tif/",basename(file),sep="")

#    outfile2b=paste("data/mcd09tif/",sub("[.]tif","_sd.tif",basename(file)),sep="")

## set sensor-specific parameters

    sensor=ifelse(grepl("MOD",file),"MOD","MYD")

    if(sensor=="MOD") {

        exts=rbind(exts,modexts)

        scanangle=-15

    }

    if(sensor=="MYD") {

        exts=rbind(exts,mydexts)

        scanangle=15

    }

    foreach(ti=1:nrow(exts)) %dopar% {

        textent=extent(exts$xmin[ti],exts$xmax[ti],exts$ymin[ti],exts$ymax[ti])

        ## extract the tile

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

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

        d=crop(raster(file),textent)

        ## acount for scale of data is 10000*CF

        d=d*.01

        ## skip null tiles - will only have this if tiles are quite small (<10 degrees)

        if(is.null(d@data@values)) return(NULL)

        ## make the gabor kernel

        ## this specifies the 'shape' of the noise we are trying to remove

        psi=fgabor(d,theta=scanangle,x=400,y=4) #3

#        psi=stack(lapply(scanangle,function(a) fgabor(d,theta=a,x=400,y=4))) #3

        ## run the correction function.  

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

        ## write the file

        if(!is.null(outfile)) writeRaster(res*100,file=toutfile,overwrite=T,datatype='INT2S',options=c("COMPRESS=LZW", "PREDICTOR=2"),NAvalue=-32768)

        ## remove temporary files

        rmr(d);rmr(psi);rmr(res)

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

    }

## create VRT of first band of the full image 

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

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

## mosaic the tiles with the original data (keeping the new data when available)

tfiles=paste(c(fvrt,list.files("tmp",pattern=paste(sub("[.]tif","",basename(outfile)),"_[0-9]*[.]tif",sep=""),full=T)),collapse=" ")

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

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

###  calculate monthly means of terra and aqua

# Create combined (MOD+MYD) uncorrected mean CF

    foreach(i=1:12) %dopar% {

        ## get files

        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),".tif",sep="")

        ## check if output already exists

        ops=paste("-t_srs 'EPSG:4326' -multi -srcnodata -32768 -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 -wo SOURCE_EXTRA=5")

        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).",

            "The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days. ",

            "Band Descriptions: 1) Mean Monthly Cloud Frequency 2) Four Times the Standard Deviation of Mean Monthly Cloud Frequency"'"),

            "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=""))

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

    }

    foreach(i=1:12) %dopar% {

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

for( i in 1:12){

    f2=list.files(paste(datadir,"/mcd09ctif",sep=""),pattern=paste(".*MCD09_",sprintf("%02d",i),"[.]tif$",sep=""),full=T)

    outfile=paste(datadir,"/mcd09ctif/",basename(file),sep="")

    outfile2=paste("data/mcd09tif/MCD09_",sprintf("%02d",i),".tif",sep="")

    outfile2b=paste("data/mcd09tif/MCD09_",sprintf("%02d",i),"_sd.tif",sep="")

    ops="-stats -scale 0 10000 0 100 -ot Byte -a_nodata 255 -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=yes"

    system(paste("gdal_translate  ",ops," ",paste("-mo ",tags,sep="",collapse=" ")," ",outfile," ",outfile2))

    system(paste("gdal_translate  -b 2 ",ops," ",paste("-mo ",tags,sep="",collapse=" ")," ",file," ",outfile2b))

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

#### stuff below here is old junk.....       

        ##  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))

        ## file.remove(temptffile)

        system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))

        ## create temporary nc file with time information to append to CF data

        cat(paste("

    netcdf time {

      dimensions:

        time = 1 ;

      variables:

        int time(time) ;

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

      time:calendar = \"gregorian\";

      time:long_name = \"time of observation\"; 

    data:

      time=",as.integer(date-as.Date("2000-01-01")),";

    }"),file=paste(tempdir(),"/",date,"_time.cdl",sep=""))

        system(paste("ncgen -o ",tempdir(),"/",date,"_time.nc ",tempdir(),"/",date,"_time.cdl",sep=""))

        ## add time dimension to ncfile and compress (deflate)

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

        ## add other attributes

        system(paste("ncrename -v Band1,CF ",ncfile,sep=""))

        system(paste("ncrename -v Band2,CFsd ",ncfile,sep=""))

        ## build correction factor explanation

        system(paste("ncatted ",

                     ## CF Mean

                     " -a units,CF,o,c,\"%\" ",

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

                                        #               " -a scale_factor,CF,o,b,\"0.1\" ",

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

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

                     " -a long_name,CF,o,c,\"Cloud Frequency (%)\" ",

                     " -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 NETCDF_VARNAME,CF,o,c,\"Cloud Frequency (%)\" ",

                     ## CF Standard Deviation

                     " -a units,CFsd,o,c,\"SD\" ",

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

                     " -a scale_factor,CFsd,o,f,\"0.25\" ",

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

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

                     " -a long_name,CFsd,o,c,\"Cloud Frequency (%) Intra-month (2000-2013) Standard Deviation\" ",

                     " -a NETCDF_VARNAME,CFsd,o,c,\"Cloud Frequency (%) Intra-month (2000-2013) Standard Deviation\" ",

                     ## global

                     " -a title,global,o,c,\"Cloud Climatology from MOD09 Cloud Mask\" ",

                     " -a institution,global,o,c,\"Jetz Lab, EEB, Yale University, New Haven, CT\" ",

                     " -a source,global,o,c,\"Derived from MOD09GA Daily Data\" ",

                     " -a comment,global,o,c,\"Developed by Adam M. Wilson (adam.wilson@yale.edu / http://adamwilson.us)\" ",

                     ncfile,sep=""))

        ## add the fillvalue attribute back (without changing the actual values)

                                        #system(paste("ncatted -a _FillValue,CF,o,b,-32768 ",ncfile,sep=""))

        if(as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))<1) {

            print(paste(ncfile," has no time, deleting"))

            file.remove(ncfile)

        }

        print(paste(basename(ncfile)," Finished"))

    }

if(ramdisk) {

    ## unmount the ram disk

    system(paste("sudo umount ",tmpfs)

}

### merge all the tiles to a single global composite

#system(paste("ncdump -h ",list.files(tempdir,pattern="mod09.*.nc$",full=T)[10]))

file.remove("tmp/mod09_2000-01-15.nc")

system(paste("cdo -O  mergetime -setrtomiss,-32768,-1 ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/cloud_monthly.nc"))

#  Overall mean

system(paste("cdo -O  timmean data/cloud_monthly.nc  data/cloud_mean.nc"))

### generate the monthly mean and sd

#system(paste("cdo -P 10 -O merge -ymonmean data/mod09.nc -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim.nc"))

system(paste("cdo  -f nc4c -O -ymonmean data/cloud_monthly.nc data/cloud_ymonmean.nc"))

## Seasonal Means

system(paste("cdo  -f nc4c -O -yseasmean data/cloud_monthly.nc data/cloud_yseasmean.nc"))

system(paste("cdo  -f nc4c -O -yseasstd data/cloud_monthly.nc data/cloud_yseasstd.nc"))

## standard deviations, had to break to limit memory usage

system(paste("cdo  -f nc4c -z zip -O -chname,CF,CF_sd -ymonstd -selmon,1,2,3,4,5,6 data/cloud_monthly.nc data/cloud_ymonsd_1-6.nc"))

system(paste("cdo  -f nc4c -z zip -O -chname,CF,CF_sd -ymonstd -selmon,7,8,9,10,11,12 data/cloud_monthly.nc data/cloud_ymonsd_7-12.nc"))

system(paste("cdo  -f nc4c -z zip -O mergetime  data/cloud_ymonsd_1-6.nc  data/cloud_ymonsd_7-12.nc data/cloud_ymonstd.nc"))

system("cdo -f nc4c -z zip  timmin data/cloud_ymonmean.nc data/cloud_min.nc")

system("cdo -f nc4c -z zip  timmax data/cloud_ymonmean.nc data/cloud_max.nc")

## standard deviation of mean monthly values give intra-annual variability

system("cdo -f nc4c -z zip -chname,CF,CFsd -timstd data/cloud_ymonmean.nc data/cloud_std_intra.nc")

## mean of monthly standard deviations give inter-annual variability 

system("cdo -f nc4c -z zip -chname,CF,CFsd -timmean data/cloud_ymonstd.nc data/cloud_std_inter.nc")

# Regressions through time by season

s=c("DJF","MAM","JJA","SON")

system(paste("cdo  -f nc4c -O regres -selseas,",s[1]," data/cloud_monthly.nc data/slope_",s[1],".nc",sep=""))

system(paste("cdo  -f nc4c -O regres -selseas,",s[2]," data/cloud_monthly.nc data/slope_",s[2],".nc",sep=""))

system(paste("cdo  -f nc4c -O regres -selseas,",s[3]," data/cloud_monthly.nc data/slope_",s[3],".nc",sep=""))

system(paste("cdo  -f nc4c -O regres -selseas,",s[4]," data/cloud_monthly.nc data/slope_",s[4],".nc",sep=""))

## Daily animations

regs=list(

    Venezuela=extent(c(-69,-59,0,7)),

    Cascades=extent(c(-122.8,-118,44.9,47)),

    Hawaii=extent(c(-156.5,-154,18.75,20.5)),

    Boliva=extent(c(-71,-63,-20,-15)),

    CFR=extent(c(17.75,22.5,-34.8,-32.6)),

    Madagascar=extent(c(46,52,-17,-12))

    )

r=1

system(paste("cdo  -f nc4c -O inttime,2012-01-15,12:00:00,7day  -sellonlatbox,",

             paste(regs[[r]]@xmin,regs[[r]]@xmax,regs[[r]]@ymin,regs[[r]]@ymax,sep=","),

             "  data/cloud_monthly.nc data/daily_",names(regs[r]),".nc",sep=""))

### Long term summaries

seasconc <- function(x,return.Pc=T,return.thetat=F) {

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

          ## Precipitation Concentration function

          ## This function calculates Precipitation Concentration based on Markham's (1970) technique as described in Schulze (1997)

          ## South Africa Atlas of Agrohydology and Climatology - R E Schulze, M Maharaj, S D Lynch, B J Howe, and B Melvile-Thomson

          ## Pages 37-38

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

          ## x is a vector of precipitation quantities - the mean for each factor in "months" will be taken,

          ## so it does not matter if the data are daily or monthly, as long as the "months" factor correctly

          ## identifies them into 12 monthly bins, collapse indicates whether the data are already summarized as monthly means.

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

          theta=seq(30,360,30)*(pi/180)                                       # set up angles for each month & convert to radians

                  if(sum(is.na(x))==12) { return(cbind(Pc=NA,thetat=NA)) ; stop}

                  if(return.Pc) {

                              rt=sqrt(sum(x * cos(theta))^2 + sum(x * sin(theta))^2)    # the magnitude of the summation

                                        Pc=as.integer(round((rt/sum(x))*100))}

                  if(return.thetat){

                              s1=sum(x*sin(theta),na.rm=T); s2=sum(x*cos(theta),na.rm=T)

                                        if(s1>=0 & s2>=0)  {thetat=abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                                        if(s1>0 & s2<0)  {thetat=180-abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                                        if(s1<0 & s2<0)  {thetat=180+abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                                        if(s1<0 & s2>0)  {thetat=360-abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                             thetat=as.integer(round(thetat))

                            }

                  if(return.thetat&return.Pc) return(c(conc=Pc,theta=thetat))

                  if(return.Pc)          return(Pc)

                  if(return.thetat)  return(thetat)

        }

## read in monthly dataset

mod09=brick("data/cloud_ymonmean.nc",varname="CF")

plot(mod09[1])

mod09_seas=calc(mod09,seasconc,return.Pc=T,return.thetat=F,overwrite=T,filename="data/mod09_seas.nc",NAflag=255,datatype="INT1U")

mod09_seas2=calc(mod09,seasconc,return.Pc=F,return.thetat=T,overwrite=T,filename="data/mod09_seas_theta.nc",datatype="INT1U")

plot(mod09_seas)