######

### proceed by raster using information at http://code.env.duke.edu/projects/mget/wiki/SinusoidalMODIS

  ## check tile

  ## list of tile outside valid region

  nodata=c('h00v00','h01v00','h02v00','h03v00','h04v00','h05v00','h06v00','h07v00','h08v00','h09v00','h10v00',

  ##

  h=as.numeric(substr(tile,2,3))

  v=as.numeric(substr(tile,5,6))

  ## Earth Width (m)

  ew=20015109.354*2

  ## Tile width or height = earth width / 36 = (20015109.354 + 20015109.354) / 36 = 1111950.5196666666 m

  ext=extent(llx,llx+(nc*cs),lly,lly+(nc*cs))

  grid=raster(ext,nrows=nc,ncol=nc,crs=proj)

  names(grid)=tile

  return(grid)

}

gs=expand.grid(h=0:35,v=0:17)

gs$tile=paste("h",sprintf("%02d",gs$h),"v",sprintf("%02d",gs$v),sep="")

modtiles=lapply(1:nrow(gs),function(i) maketile(gs$tile[i]))

names(modtiles)=gs$tile

### path to NCO

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

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

## Get list of all daily files

if(verbose) print(paste("Checking daily output in preparation for generating climatology:",tile))

## 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 title,global,o,c,\"MODIS Cloud Product (MOD35) Summaries\" ",

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

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

## Overall Means

if(verbose) print(paste("Calculating the overall mean:",tile))

system(paste(ncopath,"ncatted ",

## Monthly means

if(verbose) print(paste("Calculating the monthly means:",tile))

system(paste(ncopath,"ncatted ",

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

## Monthly standard deviation

if(verbose) print(paste("Calculating the monthly SD:",tile))

    outdir2,"/MOD35_",tile,"_daily.nc ",

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

system(paste(ncopath,"ncatted ",

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

## frequency of cloud days p(clear<90%)  

## number of observations

if(verbose) print(paste("Calculating the number of missing variables:",tile))

system(paste(ncopath,"ncatted ",

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

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

if(verbose) print(paste("Append all monthly climatologies into a single file:",tile))

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,"_ymonmiss.nc  ",tsdir,"/MOD35_",tile,"_ymon.nc",sep=""))

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

          q(status=1);

     }

platform="pleiades" 

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

if(testing){

  date="20090129"

  tile="h17v00"

  verbose=T

}

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

}

  ## Define region by importing one MOD11A1 raster.

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

  if(platform=="pleiades") {

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

  }

     ## Sensor Zenith      ## extract first bit to keep only "low angle" observations

     execGRASS("r.in.gdal",input=paste("SenZen_",bfile,sep=""),

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

     ## check for interpolation artefacts

#     execGRASS("r.stats",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,".txt",sep=""),flags=c("c"))

#     execGRASS("r.clump",input=paste("SZ_",i,sep=""),output=paste("SZ_",i,"_clump",sep=""))

#     execGRASS("r.stats",input=paste("SZ_",i,"_clump",sep=""),output="-",flags=c("c"))

                                            #

     #plot(p~value,data=vals)

#     print(sum(vals$p[vals$p>.1]))

     ## Solar Zenith      ## extract first bit to keep only "low angle" observations

     execGRASS("r.in.gdal",input=paste("SolZen_",bfile,sep=""),

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

     system(paste("r.mapcalc <<EOF

                CM_fill_",i," =  if(isnull(CM1_",i,"),1,0)

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

                SoZ_low_",i," =  if(SoZ_",i,"<8500,1,0)

                CM_dayflag_",i," =  if((CM1_",i," / 2^3) % 2==1,1,0)

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

                CMday_",i," = if(SoZ_low_",i,"==1&SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==1,CM_cloud_",i,",null())

                CMnight_",i," = if(SZ_low_",i,"==1&QA_useful_",i,"==1&CM_dayflag_",i,"==0,CM_cloud_",i,",null())

EOF",sep=""))

     if(drawplot){

       d2=stack(

#         raster(readRAST6(paste("QA_useful_",i,sep=""))),

#         raster(readRAST6(paste("CM_cloud_",i,sep=""))),

#         raster(readRAST6(paste("CM_dayflag_",i,sep=""))),

#         raster(readRAST6(paste("CMday_",i,sep=""))),

#         raster(readRAST6(paste("CM_fill_",i,sep=""))),

#         raster(readRAST6(paste("SoZ_",i,sep=""))),

         raster(readRAST6(paste("SZ_",i,sep=""))),

         )

       plot(d2,add=F)

     }

))

if(plot){

  d2=brick(list(raster(readRAST6("SZday_min")),raster(readRAST6("SZnight_min")),raster(readRAST6("CMday_daily")),raster(readRAST6("CMnight_daily"))))

  library(rasterVis)

  levelplot(sz1,col.regions=rainbow(100),at=seq(min(sz1@data@min),max(sz1@data@max),len=100))

  ### Write the files to a netcdf file

  ## create image group to facilitate export as multiband netcdf

    execGRASS("i.group",group="mod35",input=c("CMday_daily","CMnight_daily")) ; print("")

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

#      createopt=c("FORMAT=NC4","ZLEVEL=5","COMPRESS=DEFLATE","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")  #for compressed netcdf

      createopt=c("FORMAT=NC","WRITE_GDAL_TAGS=YES","WRITE_LONLAT=NO"),format="netCDF")

" -a missing_value,CMday,o,b,255 ",

" -a _FillValue,CMday,o,b,255 ",

" -a valid_range,CMday,o,b,\"0,3\" ",

" -a units,CMnight,o,c,\"Cloud Flag (0-3)\" ",

" -a missing_value,CMnight,o,b,255 ",

" -a _FillValue,CMnight,o,b,255 ",

" -a valid_range,CMnight,o,b,\"0,3\" ",

ncfile,sep=""))

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

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

      file.remove(ncfile)

    }

############  copy files to lou

#if(platform=="pleiades"){

#  archivedir=paste("MOD35/",outdir,"/",sep="")  #directory to create on lou

  ## print out some info

## delete old files

#system("cleartemp")

save(tb,file="modlandTiles.Rdata")

load("modlandTiles.Rdata")

### list of tiles to process

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

## a northern block of tiles

## subset tile corner matrix to tiles selected above

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

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

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

## write the table processed by mpiexec

tp=((!proclist$done)&proclist$avail)  #date-tiles to process

table(Available=proclist$avail,Completed=proclist$done)

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

### Now submit the script to generate the climatologies

tiles

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

### qsub script

cat(paste("

#PBS -S /bin/bash

#PBS -l walltime=2:00:00

#PBS -j n

#PBS -m be

#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

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/

system("rsync -cavv `find 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/")