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

###  R code to aquire and process MOD06_L2 cloud data from the MODIS platform

##First read in the arguments listed at the command line

args=(commandArgs(TRUE))

##args is now a list of character vectors

## First check to see if arguments are passed.

## if no args were given, print a warning and stop

if(length(args)==0) {stop("No parameters supplied, you must pass parameters")}

## Then cycle through each element of the list and evaluate the expressions.

eval(parse(text=args))

## now there is an i that corresponds to the row in the notdone object that will be processed.

## load libraries

require(sp)

require(rgdal)

require(reshape)

require(ncdf4)

require(geosphere)

require(raster)

require(spgrass6)

## specify some working directories

setwd("/nobackupp1/awilso10/mod06")

tempdir=tempdir()  # to hold intermediate files

outdir="2_daily"   # final daily product

### use MODIS tile as ROI instead

modt=readOGR("modgrid","modis_sinusoidal_grid_world",)

tiles=c("H11V8")

roi=modt[modt$HV%in%tiles,]

## Bounding box of region in lat/lon

roi_ll=spTransform(roi,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

roi_bb=bbox(roi_ll)

## read in list of all files

load("allfiles.Rdata")

load("alldates.Rdata")

load("notdonedates.Rdata")

load("vars.Rdata")

date=notdonedates[i]

#### Function to generate hegtool parameter file for multi-band HDF-EOS file

swath2grid=function(file,vars,outdir,upleft,lowright){

  print(paste("Starting file",basename(file)))

  outfile=paste(outdir,"/",basename(file),sep="")

### First write the parameter file (careful, heg is very finicky!)

  hdr=paste("NUM_RUNS = ",length(vars$varid),"|MULTI_BAND_HDFEOS:",length(vars$varid),sep="")

  grp=paste("

BEGIN

INPUT_FILENAME=",file,"

OBJECT_NAME=mod06

FIELD_NAME=",vars$variable,"|

BAND_NUMBER = 1

OUTPUT_PIXEL_SIZE_X=1000

OUTPUT_PIXEL_SIZE_Y=1000

SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )

SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )

#RESAMPLING_TYPE =",ifelse(grepl("Flag|Mask|Quality",vars),"NN","CUBIC"),"

RESAMPLING_TYPE =NN

OUTPUT_PROJECTION_TYPE = SIN

#UTM_ZONE = 10

 OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )

# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=sn_gctp

ELLIPSOID_CODE = WGS84

OUTPUT_TYPE = HDFEOS

OUTPUT_FILENAME = ",outfile,"

END

",sep="")

  ## if any remnants from previous runs remain, delete them

  if(length(list.files(tempdir(),pattern=basename(file)))>0)

    file.remove(list.files(tempdir(),pattern=basename(file),full=T))

  ## write it to a file

  cat(c(hdr,grp)    , file=paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""))

  ## now run the swath2grid tool

  ## write the tiff file

  log=system(paste("/nobackupp4/pvotava/software/heg/bin/swtif -p ",paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)

      print(paste("Finished ", file))

}

## fix out dir to tempdir

## then move on to grass processing...

#### Run the gridding procedure

lapply(fs$path[fs$dateid==date],function(file){

  swath2grid(file,vars=vars,outdir=outdir,

             upleft=paste(roi_bb[2,2],roi_bb[1,1]),

             lowright=paste(roi_bb[2,1],roi_bb[1,2]))})

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

### Import to GRASS for processing

#fs$grass=paste(fs$month,fs$year,fs$file,sep="_")

td=readGDAL(paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",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 "

projection(td)="+proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

## fucntion to convert binary to decimal to assist in identifying correct values

b2d=function(x) sum(x * 2^(rev(seq_along(x)) - 1)) #http://tolstoy.newcastle.edu.au/R/e2/help/07/02/10596.html

## for example:

b2d(c(T,T))

### create (or connect to) grass location

gisDbase="/media/data/grassdata"

gisLocation="oregon"

gisMapset="mod06"

## set Grass to overwrite

Sys.setenv(GRASS_OVERWRITE=1)

Sys.setenv(DEBUG=0)

## temporary objects to test function below

 i=1

file=paste(outdir,"/",fs$file[1],sep="")

date=as.Date("2000-05-23")

### Function to extract various SDSs from a single gridded HDF file and use QA data to throw out 'bad' observations

loadcloud<-function(date,fs){

### set up unique grass session

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

  ## set up tempfile for this PID

  initGRASS(gisBase="/usr/lib/grass64",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())

#  system(paste("g.proj -c proj4=\"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +datum=WGS84 +units=m +no_defs\"",sep=""))

    system(paste("g.proj -c proj4=\"+proj=utm +zone=10 +ellps=WGS84 +datum=WGS84 +units=m +no_defs\"",sep=""))

  ## Define region by importing one raster.

  execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",outdir,"/",fs$file[1],"\":mod06:Cloud_Mask_1km_0",sep=""),

            output="modisgrid",flags=c("quiet","overwrite","o"))

  system("g.region rast=modisgrid save=roi --overwrite")

  system("g.region roi")

  system("g.region -p")

  ## Identify which files to process

  tfs=fs$file[fs$date==date]

  nfs=length(tfs)

  ### print some summary info

  print(date)

  ## loop through scenes and process QA flags

  for(i in 1:nfs){

     file=paste(outdir,"/",tfs[i],sep="")

     ## Cloud Mask

     execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Mask_1km_0",sep=""),

              output=paste("CM1_",i,sep=""),flags=c("overwrite","o")) ; print("")

    ## extract cloudy and 'confidently clear' pixels

    system(paste("r.mapcalc <<EOF

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

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

EOF",sep=""))

    ## QA

    execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Quality_Assurance_1km_0",sep=""),

             output=paste("QA_",i,sep=""),flags=c("overwrite","o")) ; print("")

   ## QA_CER

   system(paste("r.mapcalc <<EOF

                 QA_COT_",i,"=   ((QA_",i," / 2^0) % 2^1 )==1

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

                 QA_COT3_",i,"=  ((QA_",i," / 2^3) % 2^2 )==0

                 QA_CER_",i,"=   ((QA_",i," / 2^5) % 2^1 )==1

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

EOF",sep="")) 

#                 QA_CWP_",i,"=   ((QA_",i," / 2^8) % 2^1 )==1

#                 QA_CWP2_",i,"=  ((QA_",i," / 2^9) % 2^2 )==3

   ## Optical Thickness

   execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Optical_Thickness",sep=""),

            output=paste("COT_",i,sep=""),

            title="cloud_effective_radius",

            flags=c("overwrite","o")) ; print("")

   execGRASS("r.null",map=paste("COT_",i,sep=""),setnull="-9999")

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

   system(paste("r.mapcalc \"COT_",i,"=if(QA_COT_",i,"&&QA_COT2_",i,"&&QA_COT3_",i,"&&COT_",i,">=0,COT_",i,"*0.009999999776482582,null())\"",sep=""))   

   ## set COT to 0 in clear-sky pixels

   system(paste("r.mapcalc \"COT2_",i,"=if(CM_clear_",i,"==0,COT_",i,",0)\"",sep=""))   

   ## Effective radius ##

   execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Effective_Radius",sep=""),

            output=paste("CER_",i,sep=""),

            title="cloud_effective_radius",

            flags=c("overwrite","o")) ; print("")

   execGRASS("r.null",map=paste("CER_",i,sep=""),setnull="-9999")

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

   system(paste("r.mapcalc \"CER_",i,"=if(QA_CER_",i,"&&QA_CER2_",i,"&&CER_",i,">=0,CER_",i,"*0.009999999776482582,null())\"",sep=""))   

   ## set CER to 0 in clear-sky pixels

   system(paste("r.mapcalc \"CER2_",i,"=if(CM_clear_",i,"==0,CER_",i,",0)\"",sep=""))   

   ## Cloud Water Path

#   execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",file,"\":mod06:Cloud_Water_Path",sep=""),

#            output=paste("CWP_",i,sep=""),title="cloud_water_path",

#            flags=c("overwrite","o")) ; print("")

#   execGRASS("r.null",map=paste("CWP_",i,sep=""),setnull="-9999")

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

#   system(paste("r.mapcalc \"CWP_",i,"=if(QA_CWP_",i,"&&QA_CWP2_",i,"&&CWP_",i,">=0,CWP_",i,"*0.009999999776482582,null())\"",sep=""))   

   ## set CER to 0 in clear-sky pixels

#   system(paste("r.mapcalc \"CWP2_",i,"=if(CM_clear_",i,"==0,CWP_",i,",0)\"",sep=""))   

 } #end loop through sub daily files

#### Now generate daily averages (or maximum in case of cloud flag)

  system(paste("r.mapcalc <<EOF

         COT_denom=",paste("!isnull(COT2_",1:nfs,")",sep="",collapse="+"),"

         COT_numer=",paste("if(isnull(COT2_",1:nfs,"),0,COT2_",1:nfs,")",sep="",collapse="+"),"

         COT_daily=COT_numer/COT_denom

         CER_denom=",paste("!isnull(CER2_",1:nfs,")",sep="",collapse="+"),"

         CER_numer=",paste("if(isnull(CER2_",1:nfs,"),0,CER2_",1:nfs,")",sep="",collapse="+"),"

         CER_daily=CER_numer/CER_denom

         CLD_daily=max(",paste("if(isnull(CM_cloud_",1:nfs,"),0,CM_cloud_",1:nfs,")",sep="",collapse=","),") 

EOF",sep=""))

  #### Write the files to a geotiff

  execGRASS("r.out.gdal",input="CER_daily",output=paste(tifdir,"/CER_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))

  execGRASS("r.out.gdal",input="COT_daily",output=paste(tifdir,"/COT_",format(date,"%Y%m%d"),".tif",sep=""),nodata=-999,flags=c("quiet"))

  execGRASS("r.out.gdal",input="CLD_daily",output=paste(tifdir,"/CLD_",format(date,"%Y%m%d"),".tif",sep=""),nodata=99,flags=c("quiet"))

### delete the temporary files 

  unlink_.gislock()

  system("/usr/lib/grass64/etc/clean_temp")

 system(paste("rm -R ",tf,sep=""))

### print update

  print(paste(" ###################################################################               Finished ",date,"

################################################################"))

}

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

### Now run it

tdates=sort(unique(fs$date))

done=tdates%in%as.Date(substr(list.files(tifdir),5,12),"%Y%m%d")

table(done)

tdates=tdates[!done]

mclapply(tdates,function(date) loadcloud(date,fs=fs))

#### Script to facilitate processing of MOD06 data

setwd("/nobackupp1/awilso10/mod06")

### get list of files to process

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

fs=data.frame(

  path=list.files(datadir,full=T,recursive=T,pattern="hdf"),

  file=basename(list.files(datadir,full=F,recursive=T,pattern="hdf")))

fs$date=as.Date(substr(fs$file,11,17),"%Y%j")

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

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

fs$time=substr(fs$file,19,22)

fs$datetime=as.POSIXct(strptime(paste(substr(fs$file,11,17),substr(fs$file,19,22)), '%Y%j %H%M'))

fs$dateid=format(fs$date,"%Y%m%d")

fs$path=as.character(fs$path)

fs$file=as.character(fs$file)

## get all unique dates

alldates=unique(fs$dateid)

## write it out

save(fs,file="allfiles.Rdata")

save(alldates,file="alldates.Rdata")

notdonedates=alldates

save(notdonedates,file="notdonedates.Rdata")

## output ROI

#get bounding box of region in m

#ge=SpatialPoints(data.frame(lon=c(-125,-115),lat=c(40,47)))

#projection(ge)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

#ge2=spTransform(ge, CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))

## vars

vars=as.data.frame(matrix(c(

  "Cloud_Effective_Radius",              "CER",

  "Cloud_Effective_Radius_Uncertainty",  "CERU",

  "Cloud_Optical_Thickness",             "COT",

  "Cloud_Optical_Thickness_Uncertainty", "COTU",

  "Cloud_Water_Path",                    "CWP",

  "Cloud_Water_Path_Uncertainty",        "CWPU",

  "Cloud_Phase_Optical_Properties",      "CPOP",

  "Cloud_Multi_Layer_Flag",              "CMLF",

  "Cloud_Mask_1km",                      "CM1",

  "Quality_Assurance_1km",               "QA"),

  byrow=T,ncol=2,dimnames=list(1:10,c("variable","varid"))),stringsAsFactors=F)

save(vars,file="vars.Rdata")

### Submission script

cat("

#PBS -S /bin/csh

#PBS -N cfd

# This example uses the Harpertown nodes

# User job can access ~7.6 GB of memory per Harpertown node.

# A memory intensive job that needs more than ~0.9 GB

# per process should use less than 8 cores per node

# to allow more memory per MPI process. This example

# asks for 64 nodes and 4 MPI processes per node.

# This request implies 64x4 = 256 MPI processes for the job.

#PBS -l select=64:ncpus=8:mpiprocs=4:model=har

#PBS -l walltime=4:00:00

#PBS -j oe

#PBS -W group_list=a0801

#PBS -m e

# Load some modules

module load gcc

module load hdf5

module load netcdf/4.1.3/gcc/mpt

module load mpi

module load tcl-tk/8.5.11

module load udunits/2.1.19

module load szip/2.1/gcc

module load R

module load git

# By default, PBS executes your job from your home directory.

# However, you can use the environment variable

# PBS_O_WORKDIR to change to the directory where

# you submitted your job.

cd $PBS_O_WORKDIR

# use of dplace to pin processes to processors may improve performance

# Here you request to pin processes to processors 2, 3, 6, 7 of each node.

# This helps for using the Harpertown nodes, but not for Nehalem-EP or

# Westmere-EP nodes

# The resource request of select=64 and mpiprocs=4 implies

# that you want to have 256 MPI processes in total.

# If this is correct, you can omit the -np 256 for mpiexec

# that you might have used before.

mpiexec dplace -s1 -c2,3,6,7 ./grinder < run_input > output

# It is a good practice to write stderr and stdout to a file (ex: output)

# Otherwise, they will be written to the PBS stderr and stdout in /PBS/spool,

# which has limited amount  of space. When /PBS/spool is filled up, any job

# that tries to write to /PBS/spool will die.

# -end of script-