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

## import commandline arguments

library(getopt)

## get options

opta <- getopt(matrix(c(

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

                        'verbose','v',1,'logical'

                        ), ncol=4, byrow=TRUE))

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

### directory to hold climatology

outdir2="summary" #directory for combined daily files and summarized files

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

### path to NCO

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

### Vector of variables that must be in file or they will be deleted.

###  Formated as output from system(paste("cdo -s showvar ",fdly$path[i]),intern=T)

#finalvars=" CER COT CLD"

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

## Get list of all daily files

if(verbose) print("Checking daily output in preparation for generating climatology")

 fdly=data.frame(path=list.files(outdir,pattern="nc$",full=T),stringsAsFactors=F)

  fdly$file=basename(fdly$path)

  fdly$dateid=substr(fdly$file,14,21)

  fdly$date=as.Date(substr(fdly$file,14,21),"%Y%m%d")

  fdly$month=format(fdly$date,"%m")

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

nrow(fdly)

## check validity (via npar and ntime) of nc files

#for(i in 1:nrow(fdly)){

#  fdly$ntime[i]<-as.numeric(system(paste("cdo -s ntime ",fdly$path[i]),intern=T))

#  fdly$npar[i]<-as.numeric(system(paste("cdo -s npar ",fdly$path[i]),intern=T))

#  fdly$fyear[i]<-as.numeric(system(paste("cdo -s showyear ",fdly$path[i]),intern=T))

#  fdly$fmonth[i]<-as.numeric(system(paste("cdo -s showmon ",fdly$path[i]),intern=T))

#  fdly$fvar[i]<-system(paste("cdo -s showvar ",fdly$path[i]),intern=T)

#  print(paste(i," out of ",nrow(fdly)," for year ",  fdly$fyear[i]))

#}

## print some summaries

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

print(table(fdly$year))

print(table(fdly$month))

#print(table(fdly$fvar))

## Identify which files failed test

#fdly$drop=is.na(fdly$npar)|fdly$fvar!=finalvars

## delete files that fail check?

delete=F

if(delete) {

  print(paste(sum(fdly$drop),"files will be deleted"))

  file.remove(as.character(fdly$path[fdly$drop]))

}

## remove dropped files from list

#fdly=fdly[!fdly$drop,]

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

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

## create temporary directory to put intermediate files (will be deleted when R quits)

tsdir=paste(tempdir(),"/summary",sep="")

if(!file.exists(tsdir)) dir.create(tsdir,recursive=T)

## merge all daily files to create a single file with all dates

system(paste(ncopath,"ncrcat -O ",outdir,"/*nc ",outdir2,"/MOD35_",tile,"_daily.nc",sep=""))

## Update attributes

system(paste(ncopath,"ncatted ",

" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",

" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Daily Timeseries\" ",

" -a institution,global,o,c,\"Yale University\" ",

" -a source,global,o,c,\"MODIS Cloud Mask (MOD35)\" ",

" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",

outdir2,"/MOD35_",tile,"_daily.nc",sep=""))

### produce a monthly timeseries?

#system(paste("cdo -O monmean ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_monmean.nc",sep=""))

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

##  Generate the Climatologies

if(verbose) print("Generate monthly climatologies")

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 sorttimestamp -setyear,",myear," -setday,15 -ymonmean -mulc,1.0 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmean.nc",sep=""),wait=T)

## Monthly standard deviation

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

system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonstd ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))

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

system(paste(ncopath,"ncatted ",

#" -a long_name,CER_sd,o,c,\"Cloud Particle Effective Radius (standard deviation of daily observations)\" ",

" -a long_name,CLD_sd,o,c,\"Cloud Mask (standard deviation of daily observations)\" ",

#" -a long_name,COT_sd,o,c,\"Cloud Optical Thickness (standard deviation of daily observations)\" ",

tsdir,"/MOD35_",tile,"_ymonstd.nc",sep=""))

## cld == 0

if(verbose) print("Calculating the proportion of cloudy days")

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean -mulc,1.0 -eqc,0 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD0 ",tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CLD0,o,c,\"Proportion of Days with Cloud Mask == 0\" ",

" -a units,CLD0,o,c,\"Proportion\" ",

tsdir,"/MOD35_",tile,"_ymoncld0.nc",sep=""))

## cld == 0|1

if(verbose) print("Calculating the proportion of cloudy days")

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -lec,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD01 ",tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CLD01,o,c,\"Proportion of Days with Cloud Mask == 0|1\" ",

" -a units,CLD01,o,c,\"Proportion\" ",

tsdir,"/MOD35_",tile,"_ymoncld01.nc",sep=""))

## cld == 1

if(verbose) print("Calculating the proportion of uncertain days")

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -eqc,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD1 ",tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CLD1,o,c,\"Proportion of Days with Cloud Mask == 1 (uncertain)\" ",

" -a units,CLD1,o,c,\"Proportion\" ",

tsdir,"/MOD35_",tile,"_ymoncld1.nc",sep=""))

## cld >= 2 (setting cld==01 to missing because 'uncertain')

if(verbose) print("Calculating the proportion of clear days")

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -gtc,1 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD23 ",tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CLD23,o,c,\"Proportion of Days with Cloud Mask >= 2 (Probably Clear or Certainly Clear)\" ",

" -a units,CLD23,o,c,\"Proportion\" ",

tsdir,"/MOD35_",tile,"_ymoncld2.nc",sep=""))

## cld >= 1

if(verbose) print("Calculating the proportion of clear days")

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -gec,1 -selvar,CLD2 ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD13 ",tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CLD13,o,c,\"Proportion of Days with Cloud Mask >= 1\" ",

" -a units,CLD13,o,c,\"Proportion\" ",

tsdir,"/MOD35_",tile,"_ymoncld13.nc",sep=""))

## number of observations

if(verbose) print("Calculating the number of missing variables")

system(paste("cdo -O sorttimestamp  -setyear,",myear," -setday,15 -nint -mulc,100 -ymonmean  -mulc,1.0 -eqc,9999 -setmisstoc,9999   -selvar,CLD ",outdir2,"/MOD35_",tile,"_daily.nc ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD,CLD_pmiss ",tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CLD_pmiss,o,c,\"Proportion of Days with missing data for CLD\" ",

" -a scale_factor,CLD_pmiss,o,d,0.01 ",

" -a units,CLD_pmiss,o,c,\"Proportion\" ",

tsdir,"/MOD35_",tile,"_ymonmiss.nc",sep=""))

## TODO: fix projection information so GDAL can read it correctly.

## clean up variables?

## append variables to a single file

if(verbose) print("Append all monthly climatologies into a single file")

system(paste(ncopath,"ncks -O ",tsdir,"/MOD35_",tile,"_ymonmean.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonstd.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld0.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld01.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld1.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymoncld2.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

system(paste(ncopath,"ncks -A ",tsdir,"/MOD35_",tile,"_ymonmiss.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

## append sinusoidal grid from one of input files as CDO doesn't transfer all attributes

if(verbose) print("Clean up file (update attributes, flip latitudes, add grid description")

#system(paste(ncopath,"ncea -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1],"  ",tsdir,"/sinusoidal.nc",sep=""))

#system(paste(ncopath,"ncks -A -d time,0,1 -v sinusoidal ",list.files(outdir,full=T,pattern="[.]nc$")[1],"  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

## invert latitude so it plays nicely with gdal

system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD35_",tile,"_ymon.nc ",outdir2,"/MOD35_",tile,".nc",sep=""))

## update attributes

system(paste(ncopath,"ncatted ",

#" -a standard_parallel,sinusoidal,o,c,\"0\" ",

#" -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",

#" -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",

" -a units,time,o,c,\"days since 2000-1-1 0:0:0\" ",

" -a title,global,o,c,\"MODIS Cloud Product (MOD35) Climatology\" ",

" -a institution,global,o,c,\"Yale University\" ",

" -a source,global,o,c,\"MODIS Cloud Product (MOD35)\" ",

" -a comment,global,o,c,\"Compiled by Adam M. Wilson (adam.wilson@yale.edu)\" ",

outdir2,"/MOD35_",tile,".nc",sep=""))

print(paste("###############################  Processed ",nrow(fdly),"days for tile:",tile," ###################################################"))

print("Years:")

print(table(fdly$fyear))

print("Months:")

print(table(fdly$fmonth))

## quit R

q("no")

  datadir=paste("/nobackupp4/datapool/modis/MOD35_L2.006/",year,"/",doy,"/",sep="")

  setwd("/nobackupp1/awilso10/mod35")

  ncfile=paste(outdir,"/MOD35_",tile,"_",date,".nc",sep="")  #this is the 'final' daily output file

#### Script to facilitate processing of MOD06 data

  setwd("/nobackupp1/awilso10/mod35")

library(rgdal)

library(raster)

library(RSQLite)

verbose=T

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

## delete temporary log file that can grow to GB

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

tile="h11v08"  # Venezuela

#tile="h11v07"  # Venezuela coast

#tile="h09v04"  # Oregon

tile="h21v09"  #Kenya

### list of tiles to process

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

tiles=tiles[c(1,4)]

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

### get list of files to process

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

#datadir="/nobackupp1/awilso10/mod06/data"   #for data downloaded from 

outdir="daily/" #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=do.call(rbind.data.frame,lapply(1:nrow(tile_bb),function(i){

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

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

    )

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

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

alldates=format(seq(as.Date("2000-03-01"),as.Date("2011-12-31"),1),"%Y%m%d")

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

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

## identify which have been completed

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

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

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

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

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

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

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

## write the table processed by mpiexec

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

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

### qsub script

cat(paste("

#PBS -S /bin/bash

#PBS -l select=50:ncpus=8:mpiprocs=8

##PBS -l select=2:ncpus=8:mpiprocs=8

##PBS -l select=2:ncpus=4:mpiprocs=4

#PBS -l walltime=5:00:00

#PBS -j n

#PBS -m be

#PBS -N mod35

#PBS -q normal

#PBS -V

CORES=400

HDIR=/u/armichae/pr/

#  source $HDIR/etc/environ.sh

  source /u/awilso10/environ.sh

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

## Submit it

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

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

### Now submit the script to generate the climatologies

tiles

ctiles=tiles[c(2)]  #subset to only some tiles (for example if some aren't finished yet)?

climatescript="/u/awilso10/environmental-layers/climate/procedures/MOD35_Climatology.r"

## 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 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 -l select=50:ncpus=8:mpiprocs=8

#PBS -l select=4:ncpus=4:mpiprocs=4

#PBS -l walltime=2:00:00

#PBS -j n

#PBS -m be

#PBS -N mod35_climate

#PBS -q devel

#PBS -V

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

CORES=16

HDIR=/u/armichae/pr/

#  source $HDIR/etc/environ.sh

  source /u/awilso10/environ.sh

  source /u/awilso10/.bashrc

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

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

summarydir="summary"

sumfiles=list.files("summary",pattern="^MOD06_.*[0-9][.]nc",full=T)

system(paste("scp ",paste(sumfiles,collapse=" ")," adamw@acrobates.eeb.yale.edu:/data/personal/adamw/projects/interp/data/modis/mod06/summary",sep=""))

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

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