### Process a folder of daily MOD06 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

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

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

  file=list.files(outdir,pattern="nc$"))

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

print(table(fdly$fmonth))

print(table(fdly$fvar))

## delete files that fail check?

delete=T

if(delete) {

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

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

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

  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,"/MOD06_",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 (MOD06) Daily Timeseries\" ",

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

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

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

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

### produce a monthly timeseries?

#system(paste("cdo -O monmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",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

system(paste("cdo -O sorttimestamp -setyear,",myear," -setday,15 -ymonmean ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymonmean.nc",sep=""),wait=T)

## Monthly standard deviation

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

system(paste(ncopath,"ncrename -v CER,CER_sd -v CLD,CLD_sd -v COT,COT_sd ",tsdir,"/MOD06_",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,"/MOD06_",tile,"_ymonstd.nc",sep=""))

## cer > 20

system(paste("cdo -O  sorttimestamp -setyear,",myear," -setday,15 -ymonmean -gtc,20 -selvar,CER ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))

system(paste(ncopath,"ncrename -v CER,CER20 ",tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))

system(paste(ncopath,"ncatted ",

" -a long_name,CER20,o,c,\"Proportion of Days with Cloud Particle Effective Radius > 20um\" ",

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

tsdir,"/MOD06_",tile,"_ymoncer20.nc",sep=""))

## number of observations

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

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

system(paste(ncopath,"ncatted ",

" -a long_name,CER_pmiss,o,c,\"Proportion of Days with missing data for CER\" ",

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

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

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

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

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

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

## append variables to a single file

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

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

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

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

system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymonmiss.nc  ",tsdir,"/MOD06_",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,"/MOD06_",tile,"_ymon.nc",sep=""))

## invert latitude so it plays nicely with gdal

system(paste(ncopath,"ncpdq -O -a -y ",tsdir,"/MOD06_",tile,"_ymon.nc ",outdir2,"/MOD06_",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 (MOD06) Climatology\" ",

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

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

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

outdir2,"/MOD06_",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")

## location of MOD06 files

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

### print some status messages

require(rgdal)

## path to swath database

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

## path to NCO

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

## specify working directory

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

td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""))

projection(td)="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs +datum=WGS84 +ellps=WGS84 "

## vector of variables expected to be in final netcdf file.  If these are not present, the file will be deleted at the end.

finalvars=c("CER","COT","CLD")

loadcloud<-function(date,swaths,ncfile){

  ## make temporary working directory

  tf=paste(tempdir(),"/grass", Sys.getpid(),"/", sep="")  #temporar

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

  ## create output directory if needed

  if(!file.exists(dirname(ncfile))) dir.create(dirname(ncfile))

  if(verbose) print(paste("Set up temporary grass session in",tf))

  tfs=basename(swaths)

    ## extract cloudy and 'probably/confidently clear' pixels

                CM_clear_",i," =  ((CM1_",i," / 2^0) % 2) == 1  &&  ((CM1_",i," / 2^1) % 2^2) >  2 

                 QA_COT2_",i,"=  ((QA_",i," / 2^1) % 2^2 )>=2

                 QA_CER2_",i,"=  ((QA_",i," / 2^6) % 2^2 )>=2

   ## keep only positive COT values where quality is 'useful' and '>= good' & scale to real units

   ## keep only positive CER values where quality is 'useful' and '>= good' & scale to real units

      time=",as.integer(as.Date(date,"%Y%m%d")-as.Date("2000-01-01")),";

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

  if(verbose) print("Accessing swath ID's from database")

  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," AND

                  year==",format(as.Date(date,"%Y%m%d"),"%Y")," AND

                  day==",as.numeric(format(as.Date(date,"%Y%m%d"),"%j"))

    ))

  con=dbDisconnect(con)

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

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

  ## find the swaths on disk (using datadir)

  swaths=list.files(datadir,pattern=paste(fs$id,collapse="|"),recursive=T,full=T)

  ## Run the gridding procedure

  lapply(swaths,swath2grid,vars=vars,upleft=paste(tile_bb$lat_max,tile_bb$lon_min),lowright=paste(tile_bb$lat_min,tile_bb$lon_max))

    outfiles=paste(tempdir(),"/",basename(swaths),sep="")

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

  ## run the mod06 processing for this date

  loadcloud(date,swaths=swaths,ncfile=ncfile)

  ## confirm it has all 'final variables as specified above"

  fvar=all(finalvars%in%do.call(c,strsplit(system(paste("cdo -s showvar ",ncfile),intern=T)," ")))

  if(!ntime==1&fvar) {

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

      file.remove(ncfile)

    }

system("cleartemp")

#tile="h11v07"  # Venezuela coast

#tile="h09v04"  # Oregon

tile="h21v09"  #Kenya

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

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

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

### report on what has already been processed

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

#}

## Identify which dates still need to be processed

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

notdone=alldates[!done]  

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

## write the table processed by mpiexec

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

CORES=400

LOGSTDOUT=$IDIR/log/",tile,"_stdout

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

### use mpiexec to parallelize across days

### Now process the climatologies

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

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

summarydir="summary"

system(paste("scp ",summarydir,"/MOD06_",tile,".nc 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=""))