#### Script to facilitate processing of MOD06 data

library(raster)

verbose=T

system("rm /nobackupp1/awilso10/software/heg/TOOLKIT_MTD/runtime/LogStatus")

#tile="h09v04"  # Oregon

tile="h21v09"  #Kenya

### list of tiles to process

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

datadir="/nobackupp4/datapool/modis/MOD06_L2.005/"

##find swaths in region from sqlite database for the specified date/tile

## path to swath database

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

con=dbConnect("SQLite", dbname = db)

  d=dbGetQuery(con,paste("SELECT * from swath_geo

            WHERE east>=",tile_bb$lon_min[i]," AND

                  west<=",tile_bb$lon_max[i]," AND

                  north>=",tile_bb$lat_min[i]," AND

                  south<=",tile_bb$lat_max[i])

  print(paste("Finished tile",tile_bb$tile[i]))

  return(d)

}))

  fs$id=substr(fs$id,7,19)

### Identify which swaths are available in the datapool

swaths=data.frame(path=list.files(datadir,pattern=paste("hdf$"),recursive=T,full=T),stringsAsFactors=F)  #all swaths in data pool

swaths$id=substr(basename(swaths$path),10,22)

fs$exists=fs$id%in%swaths$id 

fs$path=swaths$path[match(fs$id,swaths$id)]

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

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)

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

## identify which date-tiles have already been run

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

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

#PBS -S /bin/bash

#PBS -m be

#PBS -N mod06

  source /u/awilso10/environ.sh

IDIR=/nobackupp1/awilso10/mod06/

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

LOGSTDERR=$IDIR/log/swath_stderr

system(paste("cat mod06_qsub",sep=""))

system(paste("cat notdone.txt | head",sep=""))

system(paste("cat notdone.txt | wc -l ",sep=""))

system(paste("qsub mod06_qsub",sep=""))

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

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

delay=F

## check running jobs to get JobID of job you want to wait for

system("qstat -u awilso10")

## enter JobID here:

job="881394.pbspl1.nas.nasa.gov"

cat(paste("

#PBS -S /bin/bash

#PBS -l walltime=5:00:00

#PBS -j n

#PBS -m be

#PBS -N mod06_climate

#PBS -q normal

#PBS -V

  source /u/awilso10/environ.sh

IDIR=/nobackupp1/awilso10/mod06/

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

WORKLIST=$IDIR/notdone_climate.txt

EXE=Rscript

LOGSTDOUT=$IDIR/log/climatology_stdout

LOGSTDERR=$IDIR/log/climatology_stderr

### use mpiexec to parallelize across days

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

### Now process the climatologies

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

## check files

system(paste("cat mod06_climatology_qsub",sep=""))        #qsub submission script

system(paste("cat notdone_climate.txt | head",sep=""))    #top of job file

system(paste("cat notdone_climate.txt | wc -l ",sep=""))  #number of jobs to be run

## Submit it

system(paste("qsub mod06_climatology_qsub",sep=""))

## check progress

system("qstat -u awilso10")

# system("qsub -I -l walltime=2:00:00 -lselect=2:ncpus=16:model=san -q devel")

#system(paste("scp ",paste(fs$path[40421:40422],collapse=" ")," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/swaths",sep=""))