#### Script to facilitate processing of MOD06 data

setwd("/nobackupp1/awilso10/mod06")

library(raster)

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

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

tile="h11v08"  # Venezuela

#tile="h09v04"  # Oregon

tile="h21v09"  #Kenya

### get list of files to process

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

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

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

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

            WHERE east>=",tile_bb$lon_min," AND

                  west<=",tile_bb$lon_max," AND

                  north>=",tile_bb$lat_min," AND

                  south<=",tile_bb$lat_max)

    )

  con=dbDisconnect(con)

  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

## identify which have been completed

done=alldates%in%substr(fdone,14,21)

print(paste(table(done)[2]," out of",length(alldates),

      "(",round(100*table(done)[2]/length(alldates),1),

      "%) dates for tile",tile,

      "have been processed.  Breakdown by year of completed days:"))

print(table(substr(alldates[done],1,4)))

#updatedone=F #update the "done" list using the 

#if(updatedone&exists("fdly")){  #update using table from below

#  done[alldates%in%fdly$dateid[fdly$drop]]=F

#}

## This vector will be used to tell mpiexec which days to include

notdone=alldates[!done]  

#PBS -S /bin/bash

#PBS -m be

#PBS -N mod06

#PBS -q devel

#PBS -V

  source $HDIR/etc/environ.sh

  source /u/awilso10/.bashrc

IDIR=/nobackupp1/awilso10/mod06/

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

LOGSTDERR=$IDIR/log/",tile,"_stderr

### use mpiexec to parallelize across days

Rscript --verbose ",climatescript," --verbose T --tile \"",tile,"\"

### Check the file

# mpirun -np 1 -r ssh R --no-save

## check progress

system("qstat -u awilso10")

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