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

# Redirect all warnings to stderr()

#options(warn = -1)

#write("2) write() to stderr", stderr())

#write("2) write() to stdout", stdout())

#warning("2) warning()")

## import commandline arguments

library(getopt)

## load libraries

require(reshape)

require(geosphere)

require(raster)

require(rgdal)

require(spgrass6)

## get options

opta <- getopt(matrix(c(

                        'date', 'd', 1, 'character',

                        'tile', 't', 1, 'character',

                        'verbose','v',1,'logical',

                        'help', 'h', 0, 'logical'

                        ), ncol=4, byrow=TRUE))

if ( !is.null(opta$help) )

  {

       prg <- commandArgs()[1];

          cat(paste("Usage: ", prg,  " --date | -d <file> :: The date to process\n", sep=""));

          q(status=1);

     }

platform="pleiades" 

## default date and tile to play with  (will be overwritten below when running in batch)

  date="20090129"

  tile="h11v08"

  tile="h17v00"

  verbose=T

}

## now update using options if given

  date=opta$date  

  tile=opta$tile 

  verbose=opta$verbose  #print out extensive information for debugging?

}

year=format(as.Date(date,"%Y%m%d"),"%Y")

doy=format(as.Date(date,"%Y%m%d"),"%j")

if(platform=="pleiades"){

  datadir=paste("/nobackupp4/datapool/modis/MOD09GA.005/",year,"/",doy,"/",sep="")

  ncopath="/nasa/sles11/nco/4.0.8/gcc/mpt/bin/"

  ## specify working directory

  basedir="/nobackupp1/awilso10/mod09/" #directory to hold files temporarily before transferring to lou

  setwd(tempdir())

  ## grass database

  ## path to MOD11A1 file for this tile to align grid/extent

  td=raster(gridfile)

}

### print some status messages

if(verbose) print(paste("Processing tile",tile," for date",date))

## load tile information and get bounding box

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

## function to convert binary to decimal to assist in identifying correct values

## this is helpful when defining QA handling below, but isn't used in processing

## b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html

## for example:

## b2d(c(T,T))

  ## set Grass to overwrite

  Sys.setenv(GRASS_OVERWRITE=1)

  Sys.setenv(DEBUG=1)

  Sys.setenv(GRASS_GUI="txt")

### Extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations

## make temporary working directory

  tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")  #temporar

  if(!file.exists(tf)) dir.create(tf)

  ## create output directory if needed

  ncfile=paste(outdir,"MOD09cloud_",tile,"_",date,".nc",sep="")  #this is the 'final' daily output file

  if(verbose) print(paste("Set up temporary grass session in",tf))

  ## Define region by importing one MOD11A1 raster.

  print("Import one MOD11A1 raster to define grid")

    execGRASS("r.in.gdal",input=td@file@name,output="modisgrid",flags=c("quiet","overwrite","o"))

    system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)

  }

if(platform=="litoria"){

  execGRASS("r.in.gdal",input=paste("NETCDF:\"/home/adamw/acrobates/adamw/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc\":CER",sep=""),

            output="modisgrid",flags=c("overwrite","o"))

  system("g.region rast=modisgrid.1 save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)

}

## Identify which files to process

#tfs=list.files(tempdir(),pattern="temp.*hdf")

if(verbose) print(paste(nfs,"swaths available for processing"))

## loop through scenes and process QA flags

  for(i in 1:nfs){

     ## Read in the data from the HDFs

     execGRASS("r.in.gdal",input=paste("QA_",bfile,sep=""),

     execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),

             output=paste("SZ_",i,sep=""),flags=c("overwrite","o")) ; print("")

#     execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))

#     execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))

#     execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))

     ## write out the table of weights for use in the neighborhood analysis to identify bad pixels from swtif

     p=75  #must be odd

     mat=matrix(rep(0,p*p),nrow=p)

     mat[0.5+p/2,]=1

     cat(mat,file="weights.txt")

     execGRASS("r.neighbors",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"clump",sep=""),method="range",size=p,weight="weights.txt")  # too slow!

     system(paste("r.mapcalc \"SZ_",i,"_clump2=SZ_",i,"clump==0\"",sep=""))

#     p=-50:50

#     system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))

#     system(paste("r.mapcalc \"SZ_",i,"_clump=if(min(",paste("SZ_",i,"[0,",min(p,"]",sep="",collapse=","),")==max(",paste("SZ_",i,"[0,",p,"]",sep="",collapse=","),"),1,0)\"",sep=""))

#     vals=do.call(rbind.data.frame,strsplit(execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"),intern=T),split=" "))

#     colnames(vals)=c("value","count")

#     vals$count=as.numeric(as.character(vals$count))

#     vals$value=as.numeric(as.character(vals$value))

#     vals=na.omit(vals)

#     vals$count[vals$value==1&vals$count>10]

                                            #

     #plot(p~value,data=vals)

#     print(sum(vals$p[vals$p>.1]))

     execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),

             output=paste("SoZ_",i,sep=""),flags=c("overwrite","o")) ; print("")

     ## produce the summaries

     system(paste("r.mapcalc <<EOF

                CM_fill_",i," =  if(isnull(CM1_",i,"),1,0)

                QA_useful_",i," =  if((QA_",i," / 2^0) % 2==1,1,0)

                SZ_low_",i," =  if(SZ_",i,"_clump2==0&SZ_",i,"<6000,1,0)

                SoZ_low_",i," =  if(SoZ_",i,"<8500,1,0)

                CM_dayflag_",i," =  if((CM1_",i," / 2^3) % 2==1,1,0)

                CM_cloud_",i," =  if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,null())

                SZday_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==1,SZ_",i,",null())

                SZnight_",i," = if(SZ_",i,"_clump2==0&CM_dayflag_",i,"==0,SZ_",i,",null())

                CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null())

                CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null())

#     CM_dscore_",i," =  if((CM_dayflag_",i,"==0|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,if(SoZ_",i,">=8500,3,4))))

#     CM_nscore_",i," =  if((CM_dayflag_",i,"==1|isnull(CM1_",i,")),0,if(QA_useful_",i,"==0,1,if(SZ_",i,">=6000,2,4)))

     drawplot=F

     if(drawplot){

       d2=stack(

#         raster(readRAST6(paste("QA_useful_",i,sep=""))),

         raster(readRAST6(paste("CM1_",i,sep=""))),

#         raster(readRAST6(paste("CM_cloud_",i,sep=""))),

#         raster(readRAST6(paste("CM_dayflag_",i,sep=""))),

#         raster(readRAST6(paste("CMday_",i,sep=""))),

#         raster(readRAST6(paste("CMnight_",i,sep=""))),

#         raster(readRAST6(paste("CM_fill_",i,sep=""))),

#         raster(readRAST6(paste("SoZ_",i,sep=""))),

         raster(readRAST6(paste("SZ_",i,sep=""))),

         raster(readRAST6(paste("SZ_",i,"_clump",sep=""))),

         raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))

         )

       plot(d2,add=F)

     }

 } #end loop through sub daily files

## select lowest view angle

## use r.series to find minimum

system(paste("r.series input=",paste("SZnight_",1:nfs,sep="",collapse=",")," output=SZnight_min method=min_raster",sep=""))

system(paste("r.series input=",paste("SZday_",1:nfs,sep="",collapse=",")," output=SZday_min method=min_raster",sep=""))

## select cloud observation with lowest sensor zenith for day and night

system(

paste("r.mapcalc <<EOF

              CMday_daily=",paste(paste("if((SZday_min+1)==",1:nfs,",CMday_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse="")),"

              CMnight_daily=",paste(paste("if((SZnight_min+1)==",1:nfs,",CMnight_",1:nfs,",",sep="",collapse=" "),"null()",paste(rep(")",times=nfs),sep="",collapse=""))

))

if(plot){

  ps=1:nfs

  ps=c(12,14,17)

  sz1=brick(lapply(ps,function(i) raster(readRAST6(paste("SZnight_",i,sep="")))))

  sz_clump=brick(lapply(ps,function(i) raster(readRAST6(paste("SZ_",i,"_clump2",sep="")))))

  d=brick(lapply(ps,function(i) raster(readRAST6(paste("CMnight_",i,sep="")))))

  d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))

  library(rasterVis)

  levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100))

  levelplot(sz_clump)

  levelplot(d)

  levelplot(d2)

}

  ### Write the files to a netcdf file

  ## create image group to facilitate export as multiband netcdf

if(file.exists(ncfile)) file.remove(ncfile)  #if it exists already, delete it

  execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),

      createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")

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

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

  cat(paste("

    netcdf time {

      dimensions:

        time = 1 ;

      variables:

        int time(time) ;

      time:units = \"days since 2000-01-01 00:00:00\" ;

      time:calendar = \"gregorian\";

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

    data:

      time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";

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

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

system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))

## add other attributes

  system(paste(ncopath,"ncatted ",

" -a units,CMday,o,c,\"Cloud Flag (0-3)\" ",

" -a missing_value,CMday,o,b,255 ",

" -a _FillValue,CMday,o,b,255 ",

" -a valid_range,CMday,o,b,\"0,3\" ",

" -a long_name,CMday,o,c,\"Cloud Flag from day pixels\" ",

" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ",

" -a missing_value,CMnight,o,b,255 ",

" -a _FillValue,CMnight,o,b,255 ",

" -a valid_range,CMnight,o,b,\"0,3\" ",

" -a long_name,CMnight,o,c,\"Cloud Flag from night pixels\" ",

ncfile,sep=""))

#system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))

## Confirm that the file has the correct attributes, otherwise delete it

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

## confirm it has all 'final variables as specified above"

fvar=all(finalvars%in%strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")[[1]])

  if(ntime!=1|!fvar) {

      print(paste("FILE ERROR:  tile ",tile," and date ",date," was not outputted correctly, deleting... "))

      file.remove(ncfile)

    }

#if(platform=="pleiades"){

#  archivedir=paste("MOD35/",outdir,"/",sep="")  #directory to create on lou

#  system(paste("ssh -q bridge2 \"ssh -q lou mkdir -p ",archivedir,"\"",sep=""))

#  system(paste("ssh -q bridge2 \"scp -q ",ncfile," lou:",archivedir,"\"",sep=""))

#  file.remove(ncfile)

#  file.remove(paste(ncfile,".aux.xml",sep=""))

#}

### delete the temporary files 

#  unlink_.gislock()

#  system(paste("rm -frR ",tempdir(),sep=""))

## delete old files

## quit

q("no",status=0)