#### Script to facilitate processing of MOD06 data

## working directory

library(raster)

library(RSQLite)

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

load("modlandTiles.Rdata")

## subset to MODLAND tiles

alltiles=system("ls -r MODTILES/ | grep tif$ | cut -c1-6 | sort | uniq - ",intern=T)

tiles=tiles[tiles%in%alltiles]

## subset tile corner matrix to tiles selected above

### get list of files to process

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

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

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

rebuildswathtable=F

if(rebuildswathtable){

  ## path to swath database

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

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

  fs=do.call(rbind.data.frame,lapply(1:nrow(tile_bb),function(i){

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

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

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

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

    d$tile=tile_bb$tile[i]

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

    return(d)

  }))

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

  ## write tile-swath list to disk

  save(fs,swaths,file="swathtile.Rdata")

}

load("swathtile.Rdata")

## get all unique dates

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

#### Generate submission file

stopdate="2011-12-31"

startdate="2009-01-01"

stopdate="2009-12-31"

alldates=format(seq(as.Date(startdate),as.Date(stopdate),1),"%Y%m%d")

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

proclist$avail=paste(proclist$tile,proclist$date,sep="_")%in%paste(avail$tile,avail$date,sep="_")

fdone=data.frame(path=list.files(outdir,pattern="nc$",recursive=T))

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

print(paste(sum(!proclist$done)," out of ",nrow(proclist)," (",round(100*sum(!proclist$done)/nrow(proclist),2),"%) remain"))

### explore tile counts

#x=table(tile=proclist$tile[proclist$done],year=proclist$year[proclist$done])

#x=x[order(rownames(x)),]

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

write.table(paste("--verbose ",script," --date ",proclist$date[tp]," --verbose T --tile ",proclist$tile[tp],sep=""),

test=F

if(test){

  i=2

  time1=system.time(system(paste("Rscript --verbose ",script," --date ",proclist$date[i]," --verbose T --tile ",proclist$tile[i],sep="")))

  hours=round(length(proclist$date[tp])*142/60/60)

  hours=round(length(proclist$date[tp])*time1[3]/60/60,1)

  hours/240

  print(paste("Based on runtime of previous command, it will take",hours," hours to process the full set"))

}

cat(paste("

#PBS -S /bin/bash

##PBS -l walltime=8:00:00

#PBS -m be

#PBS -N mod35

#PBS -q devel

  source /u/awilso10/.bashrc

IDIR=/nobackupp1/awilso10/mod35/

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

WORKLIST=$IDIR/notdone.txt

EXE=Rscript

LOGSTDOUT=$IDIR/log/mod35_stdout

LOGSTDERR=$IDIR/log/mod35_stderr

### use mpiexec to parallelize across days

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

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

### Check the files

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

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

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

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

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

### Now submit the script to generate the climatologies

## this relyies on proclist above so be sure to update above before running

md=table(tile=proclist$tile[!proclist$done],year=proclist$year[!proclist$done])

mdt=names(md[md<10,])

tiles=mdt

                   system("ssh lou 'find MOD35/summary -name \"MOD35_h[0-9][0-9]v[0-9][0-9].nc\"' ",intern=T)),split="_"),function(x) x[2]))

                   system("find summary -name \"MOD35_h[0-9][0-9]v[0-9][0-9].nc\"",intern=T)),split="_"),function(x) x[2]))

## delay start until previous jobs have finished?

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"

### qsub script

cat(paste("

#PBS -S /bin/bash

#PBS -m be

#PBS -N mod35_climate

##PBS -q normal

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

  source /pleiades/u/awilso10/environ.sh

  source /pleiades/u/awilso10/.bashrc

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

WORKLIST=$IDIR/notdone_climate.txt

EXE=Rscript

LOGSTDOUT=$IDIR/log/climatology_stdout

LOGSTDERR=$IDIR/log/climatology_stderr

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

## check files

system(paste("cat mod35_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 mod35_climatology_qsub",sep=""))

## check progress

system("qstat -u awilso10")

## start interactive job on compute node for debugging

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

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

### copy the files back to Yale

#scp `find MOD35/summary -name "MOD35_h[0-9][0-9]v[0-9][0-9].nc"` adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod35/summary/

system("gdalbuildvrt MOD35C6_2009.vrt summary/*2009mean.nc ") 

system("gdal_translate -stats -co \"COMPRESS=LZW\" -of GTiff MOD35C6_2009.vrt MOD35C6_2009.tif ")              

system("scp MOD35C6_2009.tif adamw@acrobates.eeb.24.177.10.190:/Users/adamw/Downloads/")