## import commandline arguments

library(getopt)

## get options

opta <- getopt(matrix(c(

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

## default date and tile to play with

tile="h11v08"

platform="pleiades" 

verbose=T

}

if(platform="litoria"){  #if running on local server, use different paths

  ## path to some executables

  swtifpath="sudo MRTDATADIR=\"/usr/local/heg/data\" PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif"

  ## path to swath database

  ## get grid file

  td=raster(paste("~/acrobates/projects/interp/data/modis/mod06/summary/MOD06_",tile,".nc",sep=""),varname="CER")

  projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "

### print some status messages

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

## load tile information and get bounding box

load(file="modlandTiles.Rdata")

tile_bb=tb[tb$tile==tile,] ## identify tile of interest

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

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

  ## create output directory if needed

  ## set up temporary grass instance for this PID

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

  system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""),ignore.stdout=T,ignore.stderr=T)

  ## Define region by importing one MOD11A1 raster.

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

  if(platform=="pleiades") execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),

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

## Identify which files to process

tfs=basename(swaths)

## drop swaths that did not produce an output file (typically due to not overlapping the ROI)

tfs=tfs[tfs%in%list.files(tempdir())]

nfs=length(tfs)

## loop through scenes and process QA flags

  for(i in 1:nfs){

    system(paste("r.mapcalc <<EOF

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

                CM_clear_",i," =  ((CM1_",i," / 2^0) % 2) == 1  &&  ((CM1_",i," / 2^1) % 2^2) >  2 

EOF",sep=""))

     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),

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

   ## QA_CER

   execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")

   ## keep only positive COT values where quality is 'useful' and '>= good' & scale to real units

   system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,",null())\"",sep=""))   

   ## Effective radius ##

   execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),

   execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")

   ## keep only positive CER values where quality is 'useful' and '>= good' & scale to real units

   system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,",null())\"",sep=""))   

   ## Cloud Water Path

#   execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),

#            output=paste("CWP_",i,sep=""),title="cloud_water_path",

#            flags=c("overwrite","o")) ; print("")

#   execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")

#   system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,",null())\"",sep=""))   

   ## set CER to 0 in clear-sky pixels

#   system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))   

#### Now generate daily averages (or maximum in case of cloud flag)

  system(paste("r.mapcalc <<EOF

         COT_daily=int(COT_numer/COT_denom)

         CER_daily=int(CER_numer/CER_denom)

EOF",sep=""))

  ### 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="mod06",output=ncfile,type="Int16",nodata=-32768,flags=c("quiet"),

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 scale_factor,CER,o,d,0.01 -a units,CER,o,c,\"micron\" -a missing_value,CER,o,d,-32768 -a long_name,CER,o,c,\"Cloud Particle Effective Radius\" ",ncfile,sep=""))

  system(paste(ncopath,"ncatted -a scale_factor,CLD,o,d,0.01 -a units,CLD,o,c,\"none\" -a missing_value,CLD,o,d,-32768 -a long_name,CLD,o,c,\"Cloud Mask\" ",ncfile,sep=""))

### delete the temporary files