### this script generates empty tifs that align with the MODIS land tiles, but are available globally

library(raster)

######

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

maketile=function(tile,type="raster",outdir="MODTILES"){

  ## check tile

  ## list of tile outside valid region

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

    'h11v00','h12v00','h13v00','h22v00','h23v00','h24v00','h25v00','h26v00','h27v00','h28v00','h29v00','h30v00','h31v00',

    'h32v00','h33v00','h34v00','h35v00','h00v01','h01v01','h02v01','h03v01','h04v01','h05v01','h06v01','h07v01','h08v01',

    'h09v01','h10v01','h25v01','h26v01','h27v01','h28v01','h29v01','h30v01','h31v01','h32v01','h33v01','h34v01','h35v01',

    'h00v02','h01v02','h02v02','h03v02','h04v02','h05v02','h06v02','h07v02','h08v02','h27v02','h28v02','h29v02',

    'h30v02','h31v02','h32v02','h33v02','h34v02','h35v02','h00v03','h01v03','h02v03','h03v03','h04v03','h05v03','h30v03',

    'h31v03','h32v03','h33v03','h34v03','h35v03','h00v04','h01v04','h02v04','h03v04','h32v04','h33v04','h34v04','h35v04',

    'h00v05','h01v05','h34v05','h35v05','h00v06','h35v06','h00v11','h35v11','h00v12','h01v12','h34v12','h35v12','h00v13',

    'h01v13','h02v13','h03v13','h32v13','h33v13','h34v13','h35v13','h00v14','h01v14','h02v14','h03v14','h04v14','h05v14',

    'h30v14','h31v14','h32v14','h33v14','h34v14','h35v14','h00v15','h01v15','h02v15','h03v15','h04v15','h05v15','h06v15',

    'h07v15','h08v15','h27v15','h28v15','h29v15','h30v15','h31v15','h32v15','h33v15','h34v15','h35v15','h00v16','h01v16',

    'h02v16','h03v16','h04v16','h05v16','h06v16','h07v16','h08v16','h09v16','h10v16','h25v16','h26v16','h27v16','h28v16',

    'h29v16','h30v16','h31v16','h32v16','h33v16','h34v16','h35v16','h00v17','h01v17','h02v17','h03v17','h04v17','h05v17',

    'h06v17','h07v17','h08v17','h09v17','h10v17','h11v17','h12v17','h13v17','h22v17','h23v17','h24v17','h25v17','h26v17',

    'h27v17','h28v17','h29v17','h30v17','h31v17','h32v17','h33v17','h34v17','h35v17') 

  if(tile%in%nodata) {print(paste(tile," not in region with data, skipping"));return(NULL)}

  ##

  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

  tw=ew/36

  ## number of cells (1200 for 1km data)

  nc=1200

  ##1 km Cell size = tile width / cells

  cs=tw/nc

  ##X coordinate of tile lower-left corner = -20015109.354 + horizontal tile number * tile width

  llx= -(ew/2) + (h * tw)

  lly= -(ew/4) + ((17-v) * tw)

  ## projection

  proj="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" 

  ## make extent

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

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

  names(grid)=tile

  if(type=="raster") {

    writeRaster(grid,paste(outdir,"/",tile,".tif",sep=""),options=c("COMPRESS=LZW"),datatype="INT1U",overwrite=T)

    return(grid)

  }

  if(type=="extent") {

    grid2=data.frame(tile=tile,t(as.vector(extent(grid))))

    names(grid2)=c("tile","xmin", "xmax", "ymin", "ymax")

  return(grid2)

}

}

## run it and save the extents as an Rdata object for easy retrieval

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

## function to confirm that this method results in identical (e same extent, number of rows and columns, projection, resolution, and origin) rasters as MODIS LST data

checkgrid<-function(tile){

  print(tile)

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

  gridfile=list.files("/nobackupp4/datapool/modis/MOD11A2.005/2009.01.01",pattern=paste(tile,".*[.]hdf$",sep=""),recursive=T,full=T)[1]

  if(!file.exists(gridfile)) return(NULL)

  td=raster(paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_8Day_1km_LST:LST_Day_1km",sep=""))

  td2=maketile(tile)

  return(compareRaster(td,td2,res=T,orig=T,rotation=T))

}

testing=F

if(testing){

  ## make the comparison for every tile

  dat=do.call(c,lapply(gs$tile,checkgrid))

  ## summarize the results

  table(dat)

}

### make table of extents

modtiles=do.call(rbind,lapply(1:nrow(gs),function(i) maketile(gs$tile[i],type="extent")))

## get list of tiles with some land

land=read.table("http://modis-land.gsfc.nasa.gov/pdf/sn_bound_10deg.txt",skip=6,nrows=648,header=T)

land$tile=paste("h",sprintf("%02d",land$ih),"v",sprintf("%02d",land$iv),sep="")

land=land[land$lon_min!=-999,]

modtiles$land=modtiles$tile%in%land$tile

## write the table as text

write.csv(modtiles,"MODIS_Tiles.csv",row.names=F,quote=F)

import ee

from ee import mapclient

import logging

logging.basicConfig()

MY_SERVICE_ACCOUNT = '511722844190@developer.gserviceaccount.com'  # replace with your service account

MY_PRIVATE_KEY_FILE = '/Users/adamw/EarthEngine-privatekey.p12'       # replace with you private key file path

ee.Initialize(ee.ServiceAccountCredentials(MY_SERVICE_ACCOUNT, MY_PRIVATE_KEY_FILE))

image = ee.Image('srtm90_v4')

map = ee.mapclient.MapClient()

map.addOverlay(mapclient.MakeOverlay(image.getMapId({'min': 0, 'max': 3000})))

##################    Data preparation for interpolation   #######################################

############################ Extraction of station data ##########################################

#This script perform queries on the Postgres database ghcn for stations matching the             #

#interpolation area. It requires the following inputs:                                           #

# 1)the text file ofGHCND  stations from NCDC matching the database version release              #

# 2)a shape file of the study area with geographic coordinates: lonlat WGS84                     #                                                     #       

# 3)a new coordinate system can be provided as an argument                                       #

# 4)the variable of interest: "TMAX","TMIN" or "PRCP"                                            #

#                                                                                                #

#The outputs are text files and a shape file of a time subset of the database                    #

#AUTHOR: Benoit Parmentier                                                                       #

#DATE: 06/02/212                                                                                 #

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363--                                  #

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

###Loading R library and packages   

library(RPostgreSQL)

library(sp)                                             # Spatial pacakge with class definition by Bivand et al.

library(spdep)                                          # Spatial pacakge with methods and spatial stat. by Bivand et al.

library(rgdal)                                          # GDAL wrapper for R, spatial utilities

library(rgeos)                                          # Polygon buffering and other vector operations

library(reshape)

library(raster)

### Parameters and arguments

db.name <- "ghcn"                #name of the Postgres database

var <- "PRCP"                    #name of the variables to keep: TMIN, TMAX or PRCP

year_start<-"1970"               #starting year for the query (included)

year_end<-"2011"                 #end year for the query (excluded)

buffer=500                      #  size of buffer around shapefile to include in spatial query (in km), set to 0 to use no buffer

infile2<-"ghcnd-stations.txt"                             #This is the textfile of station locations from GHCND

new_proj<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"

#tile="h11v08" 

tile="h09v04"

tile="h21v09"

#out_prefix<-"20121212"                                                 #User defined output prefix

#for Adam

outpath="/home/wilson/data"

path="/home/wilson/data"

setwd(path) 

############ START OF THE SCRIPT #################

#####  Connect to Station database

drv <- dbDriver("PostgreSQL")

db <- dbConnect(drv, dbname=db.name)#,options="statement_timeout = 1m")

##### STEP 1: Select station in the study area

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

#tile_bb=tb[tb$tile==tile,] ## identify tile of interest

#roi_ll=extent(tile_bb$lon_min,tile_bb$lon_max,tile_bb$lat_min,tile_bb$lat_max)

modgrid=readOGR("/home/wilson/data/modisgrid","modis_sinusoidal_grid_world")

roi=modgrid[modgrid$h==as.numeric(substr(tile,2,3))&modgrid$v==as.numeric(substr(tile,5,6)),]

CRS_interp="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

interp_area=spTransform(roi,CRS(CRS_interp))

#####  Buffer shapefile if desired

##     This is done to include stations from outside the region in the interpolation fitting process and reduce edge effects when stiching regions

if(buffer>0){  #only apply buffer if buffer >0

  interp_area=gUnionCascaded(interp_area)  #dissolve any subparts of roi (if there are islands, lakes, etc.)

  interp_areaC=gCentroid(interp_area)       # get centroid of region

  interp_areaB=spTransform(                # buffer roi (transform to azimuthal equidistant with centroid of region for most (?) accurate buffering, add buffer, then transform to WGS84)

    gBuffer(

      spTransform(interp_area,

                  CRS(paste("+proj=aeqd +lat_0=",interp_areaC@coords[2]," +lon_0=",interp_areaC@coords[1]," +ellps=WGS84 +datum=WGS84 +units=m +no_defs ",sep=""))),

      width=buffer*1000),                  # convert buffer (km) to meters

    CRS(CRS_interp))                       # reproject back to original projection

  interp_area=interp_areaB                 # replace original region with buffered region

}

## get bounding box of study area

bbox=bbox(interp_area)

### read in station location information from database

### use the bbox of the region to include only station in rectangular region to speed up overlay

dat_stat=dbGetQuery(db, paste("SELECT id,name,latitude,longitude

                  FROM stations

                  WHERE latitude>=",bbox[2,1]," AND latitude<=",bbox[2,2],"

                  AND longitude>=",bbox[1,1]," AND longitude<=",bbox[1,2],"

                  ;",sep=""))

coordinates(dat_stat)<-c("longitude","latitude")

proj4string(dat_stat)<-CRS_interp

# Spatial query to find relevant stations

inside <- !is.na(over(dat_stat, as(interp_area, "SpatialPolygons")))  #Finding stations contained in the current interpolation area

stat_roi<-dat_stat[inside,]                                            #Finding stations contained in the current interpolation area

stat_roi<-spTransform(stat_roi,CRS(new_proj))         # Project from WGS84 to new coord. system

#Quick visualization of station locations

#plot(interp_area, axes =TRUE)

#plot(stat_roi, pch=1, col="red", cex= 0.7, add=TRUE)

#legend("topleft", pch=1,col="red",bty="n",title= "Stations",cex=1.6)

#### STEP 2: Connecting to the database and query for relevant data 

##  Query to link station location information and observations

##  Concatenate date columns into single field for easy convert to date

##  Divide value by 10 to convert to degrees C and mm

##  Subset to years in year_start -> year_end

##  Drop missing values (-9999)

##  Drop observations that failed quality control (keep only qflag==NA)

system.time(

  dd<<-dbGetQuery(db,  #save object dd (data daily)

                          paste("SELECT station,year||'-'||month||'-'||day AS date,value / 10.0 as value,latitude,longitude,elevation

                                 FROM ghcn, stations

                                 WHERE station = id

                                 AND id IN ('",paste(stat_roi$id,collapse="','"),"')

                                 AND element='",var,"'

                                 AND year>=",year_start,"

                                 AND year<",year_end,"

                                 AND value<>-9999

                                 AND qflag IS NULL

                                 ;",sep=""))

  )  ### print used time in seconds  - only taking ~30 seconds....

dd=dd[!is.na(dd$value),]  #drop any NA values

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

### STEP 3: generate monthly means for climate-aided interpolation

### first extract average daily values by year.  

system.time(

           dm<<-dbGetQuery(db,  # create dm object (data monthly)

                          paste("SELECT station,month,count(value) as count,avg(value)/10.0 as value,latitude,longitude,elevation

                                 FROM ghcn, stations

                                 WHERE station = id

                                 AND id IN ('",paste(stat_roi$id,collapse="','"),"')

                                 AND element='",var,"'

                                 AND year>=",year_start,"

                                 AND year<",year_end,"

                                 AND value<>-9999

                                 AND qflag IS NULL

                                 GROUP BY station, month,latitude,longitude,elevation

                                 ;",sep=""))

            )  ### print used time in seconds  - only taking ~30 seconds....

### drop months with fewer than 75% of the data observations

### is 75% the right number?

#dm=dm[dm$count>=(.75*10*30),]

## add the monthly data to the daily table

dd$month=as.numeric(as.character(format(as.Date(dd$date),"%m")))

dd$avgvalue=dm$value[match(paste(dd$station,dd$month),paste(dm$station,dm$month))]

dd$avgcount=dm$count[match(paste(dd$station,dd$month),paste(dm$station,dm$month))]

### Generate log-transformed ratio of daily:monthly ppt and add it to daily data

if(var=="PRCP"){

  dd$anom=log((1+dd$value)/(1+dd$avgvalue))

  dd$anom[dd$anom==Inf|dd$anom2==-Inf]=NA

}

### Generate temperature anomolies for each daily value

if(var!="PRCP"){

  dd$anom=dd$value-dd$avgvalue

}

###  Transform the subset data frame in a spatial data frame and reproject

dd_sp<-SpatialPointsDataFrame (dd[,c('longitude','latitude')],data=dd,proj=CRS(CRS_interp)) 

dd_sp<-spTransform(dd_sp,CRS(new_proj))         # Project from WGS84 to new coord. system

###  Transform the subset data frame in a spatial data frame and reproject

dm_sp<-SpatialPointsDataFrame (dm[,c('longitude','latitude')],data=dm,proj=CRS(CRS_interp)) 

dm_sp<-spTransform(dm_sp,CRS(new_proj))         # Project from WGS84 to new coord. system

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

### STEP 4: Save results and outuput in textfile and a shape file

##  Save shapefile of station locations

writeOGR(stat_roi,outpath,paste("station_location_",tile,"_",var,sep=""),driver ="ESRI Shapefile")

## save shapefile of daily observations

writeOGR(dd_sp,outpath,paste("/station_daily_",tile,"_",var,sep=""), driver ="ESRI Shapefile") #Note that the layer name is the file name without extension

### write monthly data

writeOGR(dm_sp, outpath, paste("/station_monthly_",tile,"_",var,sep=""), driver ="ESRI Shapefile") #Note that the layer name is the file name without extension

##### END OF SCRIPT ##########

## download cloud products from GEWEX for comparison

setwd("~/acrobates/adamw/projects/cloud/")

library(rasterVis)

## get PATMOS-X 1-deg data

dir1="data/gewex/"

for(y in 2000:2009) system(paste("wget -nc -P ",dir1," http://climserv.ipsl.polytechnique.fr/gewexca/DATA/instruments/PATMOSX/variables/CA/CA_PATMOSX_NOAA_0130PM_",y,".nc.gz",sep=""))

## decompress

lapply(list.files(dir1,pattern="gz",full=T),function(f) system(paste("gzip -dc < ",f," > ",dir1," ",sub(".gz","",basename(f)),sep="")))

## mergetime

system(paste("cdo -mergetime ",list.files(dir1,pattern="nc$",full=T)," ",dir1,"CA_PATMOSX_NOAA.nc",sep=""))

## Get 

# Echo out library paths

.libPaths()

# Load the R MPI package if it is not already loaded

if (!is.loaded("mpi_initialize")) {

  library("Rmpi")

}

# Echo out libraries loaded

library()

# Echo out what is loaded

search()

# Spawn as many slaves as possible

mpi.spawn.Rslaves()

# In case R exits unexpectedly, have it automatically clean up

# resources taken up by Rmpi (slaves, memory, etc...)

.Last <- function(){

  if (is.loaded("mpi_initialize")){

            if (mpi.comm.size(1) > 0){

                              print("Please use mpi.close.Rslaves() to close slaves.")

                                              mpi.close.Rslaves()

                            }

                    print("Please use mpi.quit() to quit R")

                    .Call("mpi_finalize")

          }

}

# Tell all slaves to return a message identifying themselves

mpi.remote.exec(paste("I am",mpi.comm.rank(),"of",mpi.comm.size()))

# Tell all slaves to close down, and exit the program

mpi.close.Rslaves()

mpi.quit()

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"

Copyright (C) 2012 The R Foundation for Statistical Computing

ISBN 3-900051-07-0

Platform: x86_64-unknown-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.

You are welcome to redistribute it under certain conditions.

Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.

Type 'contributors()' for more information and

'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or

'help.start()' for an HTML browser interface to help.

Type 'q()' to quit R.

> # Echo out library paths

> .libPaths()

[1] "/home1/awilso10/R/x86_64-unknown-linux-gnu-library/2.15"

[2] "/nasa/R/2.15.1/lib64/R/library"                         

> 

> # Load the R MPI package if it is not already loaded

> if (!is.loaded("mpi_initialize")) {

+   library("Rmpi")

+ }

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

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

if(verbose) print("Calculating the 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

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

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

if(verbose) print("Calculating the proportion of days with 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=""))

## cld == 0

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

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

system(paste(ncopath,"ncrename -v CLD2,CLD0 ",tsdir,"/MOD06_",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,"/MOD06_",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 -lec,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld01.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD01 ",tsdir,"/MOD06_",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,"/MOD06_",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 -eqc,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld1.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD1 ",tsdir,"/MOD06_",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,"/MOD06_",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 -gtc,1 -setctomiss,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld2.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD23 ",tsdir,"/MOD06_",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,"/MOD06_",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 -gec,1 -selvar,CLD2 ",outdir2,"/MOD06_",tile,"_daily.nc ",tsdir,"/MOD06_",tile,"_ymoncld13.nc",sep=""))

system(paste(ncopath,"ncrename -v CLD2,CLD13 ",tsdir,"/MOD06_",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,"/MOD06_",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 -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,"_ymoncld0.nc  ",tsdir,"/MOD06_",tile,"_ymon.nc",sep=""))

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

system(paste(ncopath,"ncks -A ",tsdir,"/MOD06_",tile,"_ymoncld2.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")

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

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

## connect to server of choice

#system("ssh litoria")

#R

library(sp)

library(rgdal)

library(reshape)

library(ncdf4)

library(geosphere)

library(rgeos)

library(multicore)

library(raster)

library(lattice)

library(rgl)

library(hdf5)

library(spgrass6)

### set options for Raster Package

setOptions(progress="text",timer=T)

X11.options(type="Xlib")

ncores=20  #number of threads to use

setwd("/home/adamw/personal/projects/interp")

setwd("/home/adamw/acrobates/projects/interp")

## read in shapefile as Region of Interest (roi)

roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")

roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))

### use MODIS tile as ROI instead

modt=readOGR("/home/adamw/acrobates/Global/modis_sinusoidal","modis_sinusoidal_grid_world",)

tiles=c("H9V4") #oregon

tiles=c("H11V8") #Venezuela

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)

### Downloading data from LAADSWEB

# subset by geographic area of interest

## download data from ftp site.  Unfortunately the selection has to be selected using the website and orders downloaded via ftp.

#system("wget -r --retr-symlinks ftp://ladsweb.nascom.nasa.gov/orders/500676499/ -P /home/adamw/acrobates/projects/interp/data/modis/MOD06_L2_hdf")

## specify some working directories

gdir="output/"

#datadir="data/modis/MOD06_L2_hdf"

#outdir="data/modis/MOD06_L2_hdf2"

#tifdir="/media/data/MOD06_L2_tif"

#summarydatadir="data/modis/MOD06_climatologies"

datadir="data/modis/Venezuela/MOD06"

outdir="data/modis/Venezuela/MOD06_"

#tifdir="/media/data/MOD06_L2_tif"

summarydatadir="data/modis/MOD06_climatologies"

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)

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

### Installation of hegtool

## needed 32-bit libraries and java for program to install correctly

# system(paste("hegtool -h ",fs$path[1],sep=""))

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

swath2grid=function(i=1,files,vars,outdir,upleft,lowright){

  file=fs$path[i]

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

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

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

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

  grp=paste("

BEGIN

INPUT_FILENAME=",file,"

OBJECT_NAME=mod06

FIELD_NAME=",vars,"|

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 = UTM

UTM_ZONE = 10

# OUTPUT_PROJECTION_PARAMETERS = ( 6371007.181 0.0 0.0 0.0 0.0 0.0 0.0 0.0 86400.0 0.0 1.0 0.0 0.0 0.0 0.0 )

# projection parameters from http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=is_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("sudo MRTDATADIR=/usr/local/heg/data ",

    "PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif -p ",

    paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)

      print(paste("Finished ", file))

}

### update fs with completed files

fs$complete=fs$file%in%list.files(outdir,pattern="hdf$")

table(fs$complete)

## must be run as root to access the data directory!  argh!

## run sudo once here to enter password before continuing

system('sudo ls')

#### Run the gridding procedure

mclapply(which(!fs$complete),function(fi){

  swath2grid(fi,vars=vars$variable,files=fs,

             outdir=outdir,

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

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

         mc.cores=24)

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

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

#######################################################################################33

###  Produce the monthly averages

## get list of daily files

fs2=data.frame(

  path=list.files(tifdir,full=T,recursive=T,pattern="tif$"),

  file=basename(list.files(tifdir,full=F,recursive=T,pattern="tif$")))

fs2$type=substr(fs2$file,1,3)

fs2$date=as.Date(substr(fs2$file,5,12),"%Y%m%d")

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

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

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

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

# Define type/month products

vs=expand.grid(type=unique(fs2$type),month=c("01","02","03","04","05","06","07","08","09","10","11","12"))

## identify which have been completed

done.vs=unique(do.call(rbind,strsplit(list.files(summarydatadir),"_|[.]"))[,c(1,3)])

done.vs=paste(vs$type,vs$month,sep="_")%in%paste(done.vs[,1],done.vs[,2],sep="_")

table(done.vs)

vs=vs[!done.vs,]

#tfs=fs2$path[which(fs2$month==vs$month[i]&fs2$type==vs$type[i])]

#do.call(c,mclapply(2:length(tfs),function(f) #identical(extent(raster(tfs[f-1])),extent(raster(tfs[f])))))

#raster(tfs[23])

## process the monthly summaries using the raster package

mclapply(1:nrow(vs),function(i){

  print(paste("Loading ",vs$type[i]," for month ",vs$month[i]))

  td=stack(fs2$path[which(fs2$month==vs$month[i]&fs2$type==vs$type[i])])

  print(paste("Processing Metric ",vs$type[i]," for month ",vs$month[i]))

  calc(td,mean,na.rm=T,

       filename=paste(summarydatadir,"/",vs$type[i],"_mean_",vs$month[i],".tif",sep=""),

       format="GTiff",overwrite=T)

  calc(td,sd,na.rm=T,

       filename=paste(summarydatadir,"/",vs$type[i],"_sd_",vs$month[i],".tif",sep=""),

       format="GTiff",overwrite=T)

  if(vs$type[i]%in%c("CER")) {

    ## Calculate number of days with effective radius > 20um, found to be linked to precipitating clouds

    ## (Kobayashi 2007)

    calc(td,function(x) mean(ifelse(x<20,0,1),na.rm=T),

      filename=paste(summarydatadir,"/",vs$type[i],"_P20um_",vs$month[i],".tif",sep=""),

      format="GTiff",overwrite=T)

#    td2[td2<20]=0

#    td2[td2>=20]=1

#    calc(td2,mean,na.rm=T,

#         filename=paste(summarydatadir,"/",vs$type[i],"_P20um_",vs$month[i],".tif",sep=""),

#         format="GTiff",overwrite=T)

  }  

  print(paste("Processing missing data for ",vs$type[i]," for month ",vs$month[i]))

  calc(td,function(i)

       sum(!is.na(i)),filename=paste(summarydatadir,"/",vs$type[i],"_count_",vs$month[i],".tif",sep=""),

       format="GTiff",overwrite=T)

  calc(td,function(i) sum(ifelse(i==0,0,1)),

       filename=paste(summarydatadir,"/",vs$type[i],"_clear_",vs$month[i],".tif",sep=""),format="GTiff",overwrite=T)

  gc()

}

)

vs[c(10:18,22:24,34:36),]

#### reproject to Oregon state plane for easy comparison with Benoit's data

dest=raster("data/regions/oregon/lulc/W_Layer1_ClippedTo_OR83M.rst")  # choose one to match (projection, resolution, extent)

projection(dest)=CRS("+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs") #define projection

tifs=list.files(summarydatadir,pattern="*.tif$",full=T) #get list of files to process

mclapply(tifs,function(f) projectRaster(raster(f),dest,filename=paste(dirname(f),"/OR03M_",basename(f),sep=""),overwrite=T))  # warp them

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

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

tile="h11v08"   #can move this to submit script if needed

## read in list of all files

load("allfiles.Rdata")

load("notdone.Rdata")

load("vars.Rdata")

date=notdone[i]

#file=fs$path[fs$dateid==date][1]

## print out some info

print(paste("Processing date ",date," for tile",tile))

## identify tile of interest

tile_bb=tb[tb$tile==tile,]

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

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

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

  outfile=paste(tempdir(),"/",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)

  ## clean up temporary files in working directory

#  file.remove(paste("filetable.temp_",pid,sep=""))

  print(log)

  print(paste("Finished ", file))

}

#### Run the gridding procedure

lapply(fs$path[fs$dateid==date],swath2grid,vars=vars,upleft=paste(tile_bb$lat_max,tile_bb$lon_min),lowright=paste(tile_bb$lat_min,tile_bb$lon_max))

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

gridfile=list.files("/nobackupp4/datapool/modis/MOD11A1.005/2006.01.27/",pattern=paste(tile,".*hdf$",sep=""),full=T)

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 "

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

gisLocation="tmp"

gisMapset="mod06"

## set Grass to overwrite

Sys.setenv(GRASS_OVERWRITE=1)

Sys.setenv(DEBUG=1)

Sys.setenv(GRASS_GUI="txt")

## temporary objects to test function below

# i=1

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

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

#.GRASS_CACHE=spgrass6:::.GRASS_CACHE

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

loadcloud<-function(date,fs){

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

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

  ## set up temporary grass instance for this PID

  initGRASS(gisBase="/nobackupp1/awilso10/software/grass-6.4.3svn",gisDbase=tf,SG=td,override=T,location="mod06",mapset="PERMANENT",home=tf,pid=Sys.getpid())

  system(paste("g.proj -c proj4=\"",projection(td),"\"",sep=""))

  ## Define region by importing one MOD11A1 raster.

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

  execGRASS("r.in.gdal",input=paste("HDF4_EOS:EOS_GRID:\"",gridfile,"\":MODIS_Grid_Daily_1km_LST:Night_view_angl",sep=""),

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

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

  ## Identify which files to process

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

  nfs=length(tfs)

  ## loop through scenes and process QA flags

  for(i in 1:nfs){

     file=paste(tempdir,"/",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(outdir,"/CER_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))

  execGRASS("r.out.gdal",input="COT_daily",output=paste(outdir,"/COT_",date,".tif",sep=""),nodata=-999,flags=c("quiet"))

  execGRASS("r.out.gdal",input="CLD_daily",output=paste(outdir,"/CLD_",date,".tif",sep=""),nodata=99,flags=c("quiet"))

### delete the temporary files 

  unlink_.gislock()

  system("/nobackupp1/awilso10/software/grass-6.4.3svn/etc/clean_temp")

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

### print update

}

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

### Now run it

lapply(date,function(date) loadcloud(date,fs=fs))

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

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

## quit R

q("no")

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

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

## connect to server of choice

#system("ssh litoria")

#R

library(sp)

library(rgdal)

library(reshape)

library(ncdf4)

library(geosphere)

library(rgeos)

library(multicore)

library(raster)

library(lattice)

library(rgl)

library(hdf5)

X11.options(type="Xlib")

ncores=20  #number of threads to use

setwd("/home/adamw/personal/projects/interp")

setwd("/home/adamw/acrobates/projects/interp")

roi=readOGR("data/regions/Test_sites/Oregon.shp","Oregon")

roi=spTransform(roi,CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0"))

### Downloading data from LAADSWEB

# subset by geographic area of interest

# subset: 40-47, -115--125

## download data from ftp site.  Unfortunately the selection has to be selected using the website and orders downloaded via ftp.

system("wget -r --retr-symlinks ftp://ladsweb.nascom.nasa.gov/orders/500676499/ -P /home/adamw/acrobates/projects/interp/data/modis/MOD06_L2_hdf")

gdir="output/"

datadir="data/modis/MOD06_L2_hdf"

outdir="data/modis/MOD06_L2_nc"

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$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$path=as.character(fs$path)

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

## 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=paste(c(

  "Cloud_Effective_Radius",

  "Cloud_Effective_Radius_Uncertainty",

  "Cloud_Optical_Thickness",

  "Cloud_Optical_Thickness_Uncertainty",

  "Cloud_Water_Path",

  "Cloud_Water_Path_Uncertainty",

  "Cloud_Phase_Optical_Properties",

  "Cloud_Multi_Layer_Flag",

  "Cloud_Mask_1km",

  "Quality_Assurance_1km"))

### Installation of hegtool

## needed 32-bit libraries and java for program to install correctly

#system(paste("ncl_convert2nc ",f," -i ",hdfdir," -o ",outdir," -v ",vars," -nc4c -l -cl 1 -B ",sep=""))

#system(paste("h4toh5 ",hdfdir,"/",f," ",outdir,"/",f,sep=""))

# system(paste("hegtool -h ",fs$path[1],sep=""))

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

swath2grid=function(i=1,files,vars=vars,outdir,upleft="47 -125",lowright="41 -115"){

  file=fs$path[i]

  tempfile=paste(tempdir(),"/",fs$file[i],sep="")

  outfile=sub("hdf$","nc",paste(outdir,"/",fs$file[i],sep=""))

  date=fs$date[1]

  origin=as.POSIXct("1970-01-01 00:00:00",tz="GMT")

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

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

  grp=paste("

BEGIN

INPUT_FILENAME=",file,"

OBJECT_NAME=mod06

FIELD_NAME=",vars,"|

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 = NN

OUTPUT_PROJECTION_TYPE = SIN

ELLIPSOID_CODE = WGS84

OUTPUT_PROJECTION_PARAMETERS = ( 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 0.0  )

OUTPUT_TYPE = HDFEOS

OUTPUT_FILENAME = ",tempfile,"

END

",sep="")

  ## write it to a file

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

  ## now run the swath2grid tool  - must be run as root (argh!)!

  if(file.exists(tempfile)) file.remove(tempfile)

  log=system(paste("sudo MRTDATADIR=\"/usr/local/heg/data\" ",

    "PGSHOME=/usr/local/heg/TOOLKIT_MTD PWD=/home/adamw /usr/local/heg/bin/swtif -p ",

    paste(tempdir(),"/",basename(file),"_MODparms.txt -d",sep=""),sep=""),intern=T)

#  system(paste("h5dump -H ",tempdir(),"/",basename(file),sep=""))

#    log=system(paste("swtif -p ",paste(tempdir(),"/",basename(file),"_MODparms.txt",sep=""),sep=""),intern=T)

  ## convert it to netcdf and clean up

  system(paste("NCARG_ROOT=\"/usr/local/ncl_ncarg-6.0.0\" ncl_convert2nc ",

               basename(tempfile)," -i ",tempdir()," -o ",tempdir()," -B ",sep=""))

  ncfile1=paste(tempdir(),"/",sub("hdf","nc",basename(file)),sep="")

    ## add time variable

  print("Adding time dimension")

  system(paste("ncecat -O -u time ",ncfile1,outfile))

  system(paste("ncap2 -s \'time[time]=",

               as.numeric(difftime(fs$datetime[i],origin,units="mins")),"\'  ",outfile,sep=""))

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

  print("Updating netCDF dimensions")

  ## rename dimension variables

  system(paste("ncrename -d YDim_mod06,y -d XDim_mod06,x ",outfile))

  system(paste("ncap2 -s \'x[x]=0;y[y]=0\' ",outfile))

  system(paste("ncap2 -s \'lat[x,y]=0;lon[x,y]=0\' ",outfile))

  system(paste("ncap2 -s \'sinusoidal=0\' ",outfile))

  nc=nc_open(outfile,write=T)

### Get corner coordinates and convert to cell centers

  ncd=system(paste("ncdump -h ",outfile),intern=T)

  UL= as.numeric(do.call(c,strsplit(gsub("[a-z]|[A-Z]|[\\]|[\t]|[\"]|[=]|[(]|[)]|","",

    ncd[grep("UpperLeft",ncd)]),",")))+c(500,-500)

  LR= as.numeric(do.call(c,strsplit(gsub("[a-z]|[A-Z]|[\\]|[\t]|[\"]|[=]|[(]|[)]|","",

    ncd[grep("LowerRight",ncd)]),",")))+c(-500,500)

  xvar=seq(UL[1],LR[1],by=1000)

  yvar=seq(UL[2],LR[2],by=-1000)

  ncvar_put(nc,"x",vals=xvar)

  ncvar_put(nc,"y",vals=yvar)

  ## add lat-lon grid

  grid=expand.grid(x=xvar,y=yvar)

  coordinates(grid)=c("x","y")

  projection(grid)=CRS(" +proj=sinu +lon_0=0 +x_0=0 +y_0=0")

  grid2=spTransform(grid,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

  grid=SpatialPointsDataFrame(grid,data=cbind.data.frame(coordinates(grid),coordinates(grid2)))

  gridded(grid)=T

  fullgrid(grid)=T

  colnames(grid@data)=c("x","y","lon","lat")

  xlon=as.matrix(grid["lon"])

  ylat=as.matrix(grid["lat"])

  ncvar_put(nc,"lon",vals=xlon)

  ncvar_put(nc,"lat",vals=ylat)

  nc_close(nc)

  ## update attributes

  for(v in c(vars[!grepl("Mask|Quality",vars)],paste(vars[grepl("Mask|Quality",vars)],"_0",sep=""))) {

    print(v)

    system(paste("ncatted -a coordinates,",v,",o,c,\"lat lon\" ",outfile,sep=""))

    system(paste("ncatted -a grid_mapping,",v,",o,c,\"sinusoidal\" ",outfile,sep=""))

  }

  system(paste("ncatted -a units,time,o,c,\"Minutes since ",origin,"\" ",outfile))

  system(paste("ncatted -a long_name,time,o,c,\"time\"",outfile))

  system(paste("ncatted -a units,y,o,c,\"m\" ",outfile))

  system(paste("ncatted -a long_name,y,o,c,\"y coordinate of projection\" ",outfile))

  system(paste("ncatted -a standard_name,y,o,c,\"projection_y_coordinate\" ",outfile))

  system(paste("ncatted -a units,x,o,c,\"m\" ",outfile))

  system(paste("ncatted -a long_name,x,o,c,\"x coordinate of projection\" ",outfile))

  system(paste("ncatted -a standard_name,x,o,c,\"projection_x_coordinate\" ",outfile))

  ## grid attributes

  system(paste("ncatted -a units,lat,o,c,\"degrees_north\" ",outfile))

  system(paste("ncatted -a long_name,lat,o,c,\"latitude coordinate\" ",outfile))

  system(paste("ncatted -a standard_name,lat,o,c,\"latitude\" ",outfile))

  system(paste("ncatted -a units,lon,o,c,\"degrees_east\" ",outfile))

  system(paste("ncatted -a long_name,lon,o,c,\"longitude coordinate\" ",outfile))

  system(paste("ncatted -a standard_name,lon,o,c,\"longitude\" ",outfile))

   system(paste("ncatted -a grid_mapping_name,sinusoidal,o,c,\"sinusoidal\" ",outfile))

  system(paste("ncatted -a standard_parallel,sinusoidal,o,c,\"0\" ",outfile))

  system(paste("ncatted -a longitude_of_central_meridian,sinusoidal,o,c,\"0\" ",outfile))

  system(paste("ncatted -a latitude_of_central_meridian,sinusoidal,o,c,\"0\" ",outfile))

  system(paste("cdo griddes ",outfile," > grid.txt"))

  system(paste("cdo mergegrid",outfile,

               "data/modis/MOD06_L2_nc/MOD06_L2.A2006001.2025.051.2010304104117.gscs_000500676714.nc ",

               "data/modis/MOD06_L2_nc/test.nc",sep=" "))

  print(paste("Finished ", file))

}

i=100

#### Run the gridding procedure