## default date and tile to play with  (will be overwritten below when running in batch)

#tile="h11v08"

platform="pleiades" 

verbose=T

  ## path to swath database

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

  ## specify working directory

  setwd(tempdir())

  ## grass database

  gisBase="/u/armichae/pr/grass-6.4.2/"

  ## path to MOD11A1 file for this tile to align grid/extent

  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 "

}

## get bounds of swath to keep and feed into grass when generating tile

## expand a little (0.5 deg) to ensure that there is no clipping of pixels on the edges

## tile will later be aligned with MODLAND tile so the extra will eventually be trimmed

lowright=paste(max(-90,tile_bb$lat_min-0.5),min(180,tile_bb$lon_max+0.5)) #southeast corner

## vars to process

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

  "Cloud_Mask",              "CM",

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

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

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

OBJECT_NAME=mod35

FIELD_NAME=",vars$variable,"|

BAND_NUMBER = 1

SPATIAL_SUBSET_UL_CORNER = ( ",upleft," )

SPATIAL_SUBSET_LR_CORNER = ( ",lowright," )

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

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

  ## now run the swath2grid tool

  ## write the gridded file

}  #end looping over swaths

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

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

  ## create output directory if needed

  ## Identify output file

  ## set up temporary grass instance for this PID

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

  ## Define region by importing one MOD11A1 raster.

  print("Import one MOD11A1 raster to define grid")

  if(platform=="pleiades") {

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

    system("g.region rast=modisgrid save=roi --overwrite",ignore.stdout=F,ignore.stderr=F)

  }

    system(paste("r.mapcalc <<EOF

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

                CM_cloud_",i," =  if((CM1_",i," / 2^0) % 2==1,(CM1_",i," / 2^1) % 2^2,-9999)

EOF",sep=""))

 } #end loop through sub daily files

#### Now generate daily minimum p(clear)

  system(paste("r.mapcalc <<EOF

EOF",sep=""))

## reset null values

  ### Write the files to a netcdf file

  ## create image group to facilitate export as multiband netcdf

  if(file.exists(ncfile)) file.remove(ncfile)  #if it exists already, delete it

  execGRASS("r.out.gdal",input="mod35",output=ncfile,type="Byte",nodata=255,flags=c("quiet"),

system(paste("ncgen -o ",tempdir(),"/time.nc ",tempdir(),"/time.cdl",sep=""))

system(paste(ncopath,"ncks -A ",tempdir(),"/time.nc ",ncfile,sep=""))

## add other attributes

  system(paste(ncopath,"ncatted ",

ncfile,sep=""))

#system(paste(ncopath,"ncatted -a sourcecode,global,o,c,",script," ",ncfile,sep=""))

### qsub script

cat(paste("

#PBS -S /bin/bash

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

#PBS -l walltime=2:00:00

#PBS -j n

#PBS -q devel

#PBS -V

#CORES=160

HDIR=/u/armichae/pr/

  source /u/awilso10/environ.sh

  source /u/awilso10/.bashrc

IDIR=/nobackupp1/awilso10/mod35/

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

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

### Now submit the script to generate the climatologies

tiles

ctiles=c("h10v08","h11v08","h12v08","h10v07","h11v07","h12v07")  # South America

cdone=data.frame(path="",tile="")  #use this if you want to re-run everything

cdone=data.frame(path=sapply(strsplit(basename(

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

cdone$tile=substr(basename(as.character(cdone$path)),1,6)

print(paste(length(ctiles[!ctiles%in%cdone$tile]),"Tiles still need to be processed"))

### qsub script

cat(paste("

#PBS -S /bin/bash

#PBS -j n

#PBS -m be

#PBS -N mod35_climate

##PBS -q ldan

#PBS -V

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

HDIR=/u/armichae/pr/

  source $HDIR/etc/environ.sh

  source /pleiades/u/awilso10/environ.sh

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

### copy the files back to Yale

system("ssh lou")

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

exit