### Process a folder of daily MOD35 HDF files to produce a climatology

## import commandline arguments

## get options

opta <- getopt(matrix(c(

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

tile=opta$tile #tile="h11v08"

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

### directory containing daily files

outdir=paste("daily/",tile,"/",sep="")  #directory for separate daily files

  fdly$year=format(fdly$date,"%Y")

nrow(fdly)

## print some summaries

if(verbose) print("Summary of available daily files")

print(table(fdly$year))

print(table(fdly$month))

#print(table(fdly$fvar))

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

## Combine the year-by-year files into a single daily file in the summary directory (for archiving)

if(verbose) print("Merging daily files into single file output")

myear=as.integer(max(fdly$year))  #this year will be used in all dates of monthly climatologies (and day will = 15)

## Monthly means

if(verbose) print("Calculating the monthly means")

system(paste("cdo -O -b I8 -v sorttimestamp -setyear,",myear," -setday,15 -mulc,-1 -subc,100 -ymonmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)

#months=c("01","02","03","04","05","06","07","08","09","10","11","12")

#  month="02"

#system(paste("scp summary/MOD35_",tile,".nc adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/",sep=""))

#system(paste("ncdump -h ",tsdir,"/MOD35_",tile,"_ymonmean.nc ",sep=""))

## Monthly standard deviation

if(verbose) print("Calculating the monthly SD")

system(paste(ncopath,"ncrename -v PClear,PCloud_sd ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,PCloud_sd,o,c,\"Standard Deviation of p(cloud)\" ",

date=opta$date  

tile=opta$tile 

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

## get year and doy from date

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/MOD35_L2.006/",year,"/",doy,"/",sep="")

  ## path to some executables

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

  ## path to swath database

  db="/nobackupp4/pvotava/DB/export/swath_geo.sql.sqlite3.db"

  ## specify working directory

  gisBase="/u/armichae/pr/grass-6.4.2/"

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

  gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]

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

  ## specify working directory

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

  gisBase="/usr/lib/grass64"

   ## location of MOD06 files

  datadir="~/acrobates/adamw/projects/interp/data/modis/mod35"

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

## load tile information and get bounding box

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

upleft=paste(tile_bb$lat_max,tile_bb$lon_min) #northwest corner

lowright=paste(tile_bb$lat_min,tile_bb$lon_max) #southeast corner

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

## Process the gridded files to align exactly with MODLAND tile and produce a daily summary of multiple swaths

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

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

  initGRASS(gisBase=gisBase,gisDbase=tf,SG=as(td,"SpatialGridDataFrame"),override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())

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

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

      file.remove(ncfile)

    }

  ## print out some info

print(paste(" ###################################################################               Finished ",date,

"################################################################"))

## get MODLAND tile information

tb=read.table("http://landweb.nascom.nasa.gov/developers/sn_tiles/sn_bound_10deg.txt",skip=6,nrows=648,header=T)

tb$tile=paste("h",sprintf("%02d",tb$ih),"v",sprintf("%02d",tb$iv),sep="")

save(tb,file="modlandTiles.Rdata")

load("modlandTiles.Rdata")

#tile="h09v04"  # Oregon

tile="h21v09"  #Kenya

### list of tiles to process

tiles=c("h11v08","h21v09","h08v04","h09v04","h08v05","h09v05","h20v11","h31v11")

tile_bb=tb[tb$tile%in%tiles,]

### get list of files to process

datadir="/nobackupp4/datapool/modis/MOD35_L2.006/"

outdir="daily/" #paste("daily/",tile,sep="")

if(verbose) print(paste("###############",nrow(fs)," swath IDs recieved from database"))

## get all unique dates

fs$dateid=format(as.Date(paste(fs$year,fs$day,sep=""),"%Y%j"),"%Y%m%d")

alldates=unique(fs$dateid[fs$exists])

#### Generate submission file

proclist=expand.grid(date=alldates,tile=tiles)

proclist$year=substr(proclist$date,1,4)

## identify which have been completed

fdone$date=substr(basename(as.character(fdone$path)),14,21)

fdone$tile=substr(basename(as.character(fdone$path)),7,12)

proclist$done=paste(proclist$tile,proclist$date,sep="_")%in%substr(basename(as.character(fdone$path)),7,21)

### report on what has already been processed

script="/u/awilso10/environmental-layers/climate/procedures/MOD35_L2_process.r"

## write the table processed by mpiexec

file=paste("notdone.txt",sep=""),row.names=F,col.names=F,quote=F)

### qsub script

### Now submit the script to generate the climatologies

tiles

## write the table processed by mpiexec

write.table(paste("--verbose ",climatescript," --verbose T --tile ",ctiles,sep=""),

### qsub script

cat(paste("

#PBS -S /bin/bash

#PBS -j n

#PBS -m be

#PBS -N mod35_climate

#PBS -V

",if(delay) paste("#PBS -W depend=afterany:",job,sep="")," 

CORES=16

IDIR=/nobackupp1/awilso10/mod35/

##WORKLIST=$HDIR/var/run/pxrRgrs/work.txt

WORKLIST=$IDIR/notdone_climate.txt

EXE=Rscript

LOGSTDOUT=$IDIR/log/climatology_stdout

LOGSTDERR=$IDIR/log/climatology_stderr

mpiexec -np $CORES pxargs -a $WORKLIST -p $EXE -v -v -v --work-analyze 1> $LOGSTDOUT 2> $LOGSTDERR

",sep=""),file=paste("mod35_climatology_qsub",sep=""))