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

### Script to download and process the NDP-026D station cloud dataset

### to validate MODIS cloud frequencies

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

library(multicore)

library(doMC)

library(rasterVis)

library(rgdal)

library(reshape)

## Data available here http://cdiac.ornl.gov/epubs/ndp/ndp026d/ndp026d.html

## Get station locations

system("wget -N -nd http://cdiac.ornl.gov/ftp/ndp026d/cat01/01_STID -P data/")

st=read.table("data/01_STID",skip=1)

colnames(st)=c("StaID","LAT","LON","ELEV","ny1","fy1","ly1","ny7","fy7","ly7","SDC","b5c")

st$lat=st$LAT/100

st$lon=st$LON/100

st$lon[st$lon>180]=st$lon[st$lon>180]-360

st=st[,c("StaID","ELEV","lat","lon")]

colnames(st)=c("id","elev","lat","lon")

write.csv(st,"stations.csv",row.names=F)

coordinates(st)=c("lon","lat")

projection(st)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

st@data[,c("lon","lat")]=coordinates(st)

## download data

system("wget -N -nd ftp://cdiac.ornl.gov/pub/ndp026d/cat67_78/* -A '.tc.Z' -P data/")

system("gunzip data/*.Z")

## define FWF widths

f162=c(5,5,4,7,7,7,4) #format 162

c162=c("StaID","YR","Nobs","Amt","Fq","AWP","NC")

## use monthly timeseries

cld=do.call(rbind.data.frame,mclapply(sprintf("%02d",1:12),function(m) {

  d=read.fwf(list.files("data",pattern=paste("MNYDC.",m,".tc$",sep=""),full=T),skip=1,widths=f162)

  colnames(d)=c162

  d$month=as.numeric(m)

  print(m)

  return(d)}

  ))

## add lat/lon

cld[,c("lat","lon")]=coordinates(st)[match(cld$StaID,st$id),]

## drop missing values

cld=cld[,!grepl("Fq|AWP|NC",colnames(cld))]

cld$Amt[cld$Amt<0]=NA

cld$Amt=cld$Amt/100

## calculate means and sds for full record (1970-2009)

Nobsthresh=20 #minimum number of observations to include 

cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){

  data.frame(

      month=x$month[1],

      StaID=x$StaID[1],

      cld_all=mean(x$Amt[x$Nobs>=Nobsthresh],na.rm=T),  # full record

      cldsd_all=sd(x$Amt[x$Nobs>=Nobsthresh],na.rm=T),

      cld=mean(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T), #only MODIS epoch

      cldsd=sd(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T))}))

cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]

## add the MOD09 data to cld

#### Evaluate MOD35 Cloud data

mod09=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc")

mod09std=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonstd.nc")

## overlay the data with 32km diameter (16km radius) buffer

## buffer size from Dybbroe, et al. (2005) doi:10.1175/JAM-2189.1.

buf=16000

bins=cut(st$lat,10)

rerun=F

if(rerun&file.exists("valid.csv")) file.remove("valid.csv")

mod09sta=lapply(levels(bins),function(lb) {

  l=which(bins==lb)

  ## mean

  td=extract(mod09,st[l,],buffer=buf,fun=mean,na.rm=T,df=T)

  td$id=st$id[l]

  td$type="mean"

  ## std

  td2=extract(mod09std,st[l,],buffer=buf,fun=mean,na.rm=T,df=T)

  td2$id=st$id[l]

  td2$type="sd"

  print(lb)#as.vector(c(l,td[,1:4])))

  write.table(rbind(td,td2),"valid.csv",append=T,col.names=F,quote=F,sep=",",row.names=F)

  td

})#,mc.cores=3)

## read it back in

mod09st=read.csv("valid.csv",header=F)[,-c(1)]

colnames(mod09st)=c(names(mod09),"id","type")

mod09stl=melt(mod09st,id.vars=c("id","type"))

mod09stl[,c("year","month")]=do.call(rbind,strsplit(sub("X","",mod09stl$variable),"[.]"))[,1:2]

mod09stl$value[mod09stl$value<0]=NA

mod09stl=cast(mod09stl,id+year+month~type,value="value")

## add it to cld

cldm$mod09=mod09stl$mean[match(paste(cldm$StaID,cldm$month),paste(mod09stl$id,as.numeric(mod09stl$month)))]

cldm$mod09sd=mod09stl$sd[match(paste(cldm$StaID,cldm$month),paste(mod09stl$id,as.numeric(mod09stl$month)))]

## LULC

#system(paste("gdalwarp -r near -co \"COMPRESS=LZW\" -tr ",paste(res(mod09),collapse=" ",sep=""),

#             "-tap -multi -t_srs \"",   projection(mod09),"\" /mnt/data/jetzlab/Data/environ/global/landcover/MODIS/MCD12Q1_IGBP_2005_v51.tif ../modis/mod12/MCD12Q1_IGBP_2005_v51.tif"))

lulc=raster("~/acrobates/adamw/projects/interp/data/modis/mod12/MCD12Q1_IGBP_2005_v51.tif")

require(plotKML); data(worldgrids_pal)  #load IGBP palette

IGBP=data.frame(ID=0:16,col=worldgrids_pal$IGBP[-c(18,19)],stringsAsFactors=F)

IGBP$class=rownames(IGBP);rownames(IGBP)=1:nrow(IGBP)

levels(lulc)=list(IGBP)

## function to get modal lulc value

Mode <- function(x) {

      ux <- na.omit(unique(x))

        ux[which.max(tabulate(match(x, ux)))]

      }

lulcst=extract(lulc,st,fun=Mode,buffer=buf,df=T)

colnames(lulcst)=c("id","lulc")

## add it to cld

cldm$lulc=lulcst$lulc[match(cldm$StaID,lulcst$id)]

cldm$lulcc=IGBP$class[match(cldm$lulc,IGBP$ID)]

### Add biome data

if(!file.exists("../teow/biomes.shp")){

    teow=readOGR("/mnt/data/jetzlab/Data/environ/global/teow/official/","wwf_terr_ecos")

    teow=teow[teow$BIOME<90,]

    biome=unionSpatialPolygons(teow,teow$BIOME, threshold=5)

    biomeid=read.csv("/mnt/data/jetzlab/Data/environ/global/teow/official/biome.csv",stringsAsFactors=F)

    biome=SpatialPolygonsDataFrame(biome,data=biomeid[as.numeric(row.names(biome)),])

    writeOGR(biome,"../teow","biomes",driver="ESRI Shapefile",overwrite=T)

}

biome=readOGR("../teow/","biomes")

projection(biome)=projection(st)

#st$biome=over(st,biome,returnList=F)$BiomeID

dists=apply(gDistance(st,biome,byid=T),2,which.min)

st$biome=biome$BiomeID[dists]

cldm$biome=st$biome[match(cldm$StaID,st$id)]

## write out the tables

write.csv(cld,file="cld.csv",row.names=F)

write.csv(cldm,file="cldm.csv",row.names=F)

writeOGR(st,dsn=".",layer="stations",driver="ESRI Shapefile",overwrite_layer=T)

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

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

// MODIS LAND Observation Frequency

// Proportion of days in the month with at least 1 MODIS observation

var p_07 = ee.ImageCollection("MOD09GA").

filter(ee.Filter.calendarRange(2011,2012,"year")).

filter(ee.Filter.calendarRange(7,7,"month")).map(function(img) {

  return img.select(['num_observations_1km']).gte(1)}).mean().multiply(ee.Image(100)).select([0],["p07"]); 

// mean number of observations per day

var n_07 = ee.ImageCollection("MOD09GA").

filter(ee.Filter.calendarRange(2011,2012,"year")).

filter(ee.Filter.calendarRange(7,7,"month")).map(function(img) {

  return img.select(['num_observations_1km'])}).mean().select([0],["n07"]);

var n_12 = ee.ImageCollection("MOD09GA").

filter(ee.Filter.calendarRange(2011,2012,"year")).

filter(ee.Filter.calendarRange(12,12,"month")).map(function(img) {

  return img.select(['num_observations_1km']).gte(1)}).mean().multiply(ee.Image(100)).select([0],["p12"]); 

// mean number of observations per day

var p_12 = ee.ImageCollection("MOD09GA").

filter(ee.Filter.calendarRange(2011,2012,"year")).

filter(ee.Filter.calendarRange(12,12,"month")).map(function(img) {

  return img.select(['num_observations_1km'])}).mean().select([0],["n12"]);

// Combine anything into image for download

var image=n_07.addBands(p_07).addBands(n_12).addBands(p_12)

//var reproj = image.reproject('EPSG:4326', [0.008333333333, 0, -180, 0, -0.008333333333, -90])

//var palette="000000,00FF00,FF0000"

//addToMap(image.select([0]),{min:0,max:100,palette:palette},"July");

//addToMap(n_12,{min:0,max:100,palette:palette},"Decenber");

var driveFolder="EarthEngineOutput"

exportImage(image, 

  'PObs',

  {'name': 'PObs',

  'maxPixels':1000000000,

  'driveFolder':driveFolder,

  'crs': 'EPSG:4326', //4326

  'crsTransform':[0.008333333333, 0, -180, 0, -0.008333333333, -90],

//  'scale': 0.008333333333,

  'region': '[[-180, -89], [-180, 89], [180, 89], [180, -89]]'  //global

//  'region': '[[-180, 0], [-180, 10], [180, 10], [180, 0]]'  //

});

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

#// EVI_Cloud_month

#// Make land mask

#// Select the forest classes from the MODIS land cover image and intersect them

#// with elevations above 1000m.

elev = ee.Image('srtm90_v4');

cover = ee.Image('MCD12Q1/MCD12Q1_005_2001_01_01').select('Land_Cover_Type_1');

blank = ee.Image(0);

#// Where (cover == 0) and (elev > 0), set the output to 1.

land = blank.where(

    cover.neq(0).and(cover.neq(15)),//.and(elev.gt(0)),

    1);

palette = ["aec3d4", // water

               "152106", "225129", "369b47", "30eb5b", "387242", // forest

               "6a2325", "c3aa69", "b76031", "d9903d", "91af40",  // shrub, grass, savanah 

               "111149", // wetlands

               "cdb33b", // croplands

               "cc0013", // urban

               "33280d", // crop mosaic

               "d7cdcc", // snow and ice

               "f7e084", // barren

               "6f6f6f"].join(',');// tundra

#// make binary map for forest/nonforest

forest = blank.where(cover.gte(1).and(cover.lte(5)),1);

#//addToMap(forest, {min: 0, max: 1});

#//addToMap(cover, {min: 0, max: 17, palette:palette});

#// MODIS EVI Collection

collection = ee.ImageCollection("MCD43A4_EVI");

#//define reducers and filters

COUNT = ee.call("Reducer.count");

#// Loop through months and get monthly % cloudiness

#//for (var month = 1; month < 2; month += 1) {

#//month=1;

FilterMonth = ee.Filter(ee.call("Filter.calendarRange",

    start=month,end=month,field="month"));

tmonth=collection.filter(FilterMonth)

n=tmonth.getInfo().features.length; #// Get total number of images

print(n+" Layers in the collection for month "+month)

tcloud=tmonth.reduce(COUNT).float()

c=ee.Image(n);  #// make raster with constant value of n

c1=ee.Image(-1);  #// make raster with constant value of -1 to convert to % cloudy

#/////////////////////////////////////////////////

#// Generate the cloud frequency surface:

#// 1 Calculate the number of days with measurements

#// 2 divide by the total number of layers

#// 3 Add -1 and multiply by -1 to invert to % cloudy

#// 4 Rename to "PCloudy_month"

tcloud = tcloud.divide(c).add(c1).multiply(c1).expression("b()*100").int8().select_([0],["PCloudy_"+month]);

#//if(month==1) {var cloud=tcloud}  // if first year, make new object

#//if(month>1)  {var cloud=cloud.addBands(tcloud)}  // if not first year, append to first year

#//} //end loop over months

# // end loop over years

#//print(evi_miss.stats())

#//addToMap(tcloud,{min:0,max:100,palette:"000000,00FF00,FF0000"});

#//addToMap(elev,{min:0,max:2500,palette:"000000,00FF00,FF0000"});

#//centerMap(-122.418, 37.72, 10);

path = tcloud.getDownloadURL({

  'scale': 1000,

  'crs': 'EPSG:4326',

  'region': '[[-90, 0], [-90, 20], [-50, 0], [-50, 20]]'  //h11v08

#//  'region': '[[-180, -90], [-180, 90], [180, -90], [180, 90]]'

});

print('https://earthengine.googleapis.com' + path);

This folder holds scripts related to work on Google Earth Engine, primarily for processing MODIS cloud data.

#!/usr/bin/env python

## Example script that downloads data from Google Earth Engine using the python API

## MODIS MOD09GA data is processed to extract the MOD09 cloud flag and calculate monthly cloud frequency

## import some libraries

import ee

from ee import mapclient

import ee.mapclient 

import datetime

import wget

import os

import sys

import subprocess

import time

import logging

logging.basicConfig(filename='error.log',level=logging.DEBUG)

def Usage():

    print('Usage: ee.MOD9.py -projwin  ulx uly urx ury lrx lry llx lly -year year -month month -regionname 1') 

    sys.exit( 1 )

ulx = float(sys.argv[2])

uly = float(sys.argv[3])

urx = float(sys.argv[4])

ury = float(sys.argv[5])

lrx = float(sys.argv[6])

lry = float(sys.argv[7])

llx = float(sys.argv[8])

lly = float(sys.argv[9])

year = int(sys.argv[11])

month = int(sys.argv[13])

regionname = str(sys.argv[15])

#```

#ulx=-159

#uly=20

#lrx=-154.5

#lry=18.5

#year=2001

#month=6

#```

## Define output filename

output=regionname+'_'+str(year)+'_'+str(month)

## set working directory (where files will be downloaded)

cwd='/mnt/data2/projects/cloud/mod09/'

## Wget always starts with a temporary file called "download", so move to a subdirectory to

## prevent overwrting when parallel downloading 

if not os.path.exists(cwd+output):

    os.makedirs(cwd+output)

os.chdir(cwd+output)

      # output filename

unzippedfilename=output+".mod09.tif"

      # Check if file already exists and continue if so...

if(os.path.exists(unzippedfilename)):

    sys.exit("File exists:"+output)    

## initialize GEE

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

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

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

## set map center to speed up viewing

#ee.mapclient.centerMap(-121.767, 46.852, 11)

#///////////////////////////////////

#// Function to extract cloud flags

def getmod09(img): return(img.select(['state_1km']).expression("((b(0)/1024)%2)").gte(1)); 

#////////////////////////////////////////////////////

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

# Processing Function

# MOD09 internal cloud flag for this year-month

      # to filter by a date range:  filterDate(datetime.datetime(yearstart,monthstart,1),datetime.datetime(yearstop,monthstop,31))

mod09 = ee.ImageCollection("MOD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);

#      myd09 = ee.ImageCollection("MYD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);

      # calculate mean cloudiness (%), rename band, multiply by 100, and convert to integer

mod09a=mod09.mean().select([0], ['mod09']).multiply(ee.Image(1000)).int16();

#      myd09a=myd09.mean().select([0], ['MYD09_'+str(year)+'_'+str(month)]).multiply(ee.Image(100)).int8();

## Set data equal to whatver you want downloaded

data=mod09a

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

## define region for download

region=[llx, lly],[lrx, lry],[urx, ury],[ulx,uly]  #h11v08

strregion=str(list(region))

# add to plot to confirm it's working

#ee.mapclient.addToMap(data, {'range': '0,100'}, 'MOD09')

      # build the URL and name the object (so that when it's unzipped we know what it is!)

path =mod09a.getDownloadUrl({

        'crs': 'SR-ORG:6974',          # MODIS Sinusoidal

        'scale': '926.625433055833',   # MODIS ~1km

        'name': output,  # name the file (otherwise it will be a uninterpretable hash)

        'region': strregion                        # region defined above

        });

# print info to confirm there is data

print(data.getInfo())

print(' Processing.... '+output+'     Coords:'+strregion)

print(path)

#test=wget.download(path)

test=subprocess.call(['-c','-q','--timeout=0','--ignore-length','--no-http-keep-alive','-O','mod09.zip',path],executable='wget')

print('download sucess for tile '+output+':   '+str(test))

## Sometimes EE will serve a corrupt zipped file with no error

# try to unzip it

zipstatus=subprocess.call("unzip -o mod09.zip",shell=True)

if zipstatus==0:  #if sucessful, quit

    os.remove("mod09.zip")

    sys.exit("Finished:  "+output)

# if file doesn't exists or it didn't unzip, try again      

if zipstatus!=0:

    print("Zip Error for:  "+output+"...         Trying again (#2)")    

    time.sleep(15)

    test=subprocess.call(['-c','-q','--timeout=0','--ignore-length','--no-http-keep-alive','-O','mod09.zip',path],executable='wget')

zipstatus=subprocess.call("unzip -o mod09.zip",shell=True)

if zipstatus==0:  #if sucessful, quit

    os.remove("mod09.zip")

    sys.exit("Finished:  "+output)

# if file doesn't exists or it didn't unzip, try again      

if zipstatus!=0:

    print("Zip Error for:  "+output+"...         Trying again (#3)")    

    time.sleep(30)

    test=subprocess.call(['-c','-q','--timeout=0','--ignore-length','--no-http-keep-alive','-O','mod09.zip',path],executable='wget')

# try again #4

zipstatus=subprocess.call("unzip -o mod09.zip",shell=True)

if zipstatus==0:  #if sucessful, quit

    os.remove("mod09.zip")

    sys.exit("Finished:  "+output)

# if file doesn't exists or it didn't unzip, try again      

if zipstatus!=0:

    print("Zip Error for:  "+output+"...         Trying again (#4)")    

    time.sleep(30)

    test=subprocess.call(['-c','-q','--timeout=0','--ignore-length','--no-http-keep-alive','-O','mod09.zip',path],executable='wget')

zipstatus=subprocess.call("unzip -o mod09.zip",shell=True)

if zipstatus==0:  #if sucessful, quit

    os.remove("mod09.zip")

    sys.exit("Finished:  "+output)

## delete the zipped file (the unzipped version is kept)

os.remove("mod09.zip")

sys.exit("Zip Error for:  "+output)    

#!/usr/bin/env python

## Example script that downloads data from Google Earth Engine using the python API

## MODIS MOD09GA data is processed to extract the MOD09 cloud flag and calculate monthly cloud frequency

## import some libraries

import ee

from ee import mapclient

import ee.mapclient 

import datetime

import wget

import os

import sys

def Usage():

    print('Usage: ee.MOD9.py -projwin  ulx uly lrx lry -o output   ') 

    sys.exit( 1 )

ulx = float(sys.argv[2])

uly = float(sys.argv[3])

lrx = float(sys.argv[4])

lry = float(sys.argv[5])

output = sys.argv[7]

print ( output , ulx ,uly ,lrx , lry )

#import logging

MY_SERVICE_ACCOUNT = '364044830827-ubb6ja607b8j7t8m9uooi4c01vgah4ms@developer.gserviceaccount.com'  # replace with your service account

MY_PRIVATE_KEY_FILE = '/home/selv/GEE/fe3f13d90031e3eedaa9974baa6994e467b828f7-privatekey.p12'      # replace with you private key file path

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

#///////////////////////////////////

#// Function to extract cloud flags

def getmod09(img): return(img.select(['state_1km']).expression("((b(0)/1024)%2)"));

## set a year-month if you don't want to run the loop (for testing)

year=2001

month=1

print('Processing '+str(year)+'_'+str(month))

## MOD09 internal cloud flag for this year-month

## to filter by a date range:  filterDate(datetime.datetime(yearstart,monthstart,1),datetime.datetime(yearstop,monthstop,31))

mod09 = ee.ImageCollection("MOD09GA").filter(ee.Filter.calendarRange(year,year,"year")).filter(ee.Filter.calendarRange(month,month,"month")).map(getmod09);

## calculate mean cloudiness (%), rename band, multiply by 100, and convert to integer

mod09a=mod09.mean().select([0], ['MOD09_'+str(year)+'_'+str(month)]).multiply(ee.Image(100)).byte();

## print info to confirm there is data

# mod09a.getInfo()

## define region for download

region=[ulx,lry ], [ulx, uly], [lrx, uly], [lrx, lry]  #h11v08

strregion=str(list(region))

## Define tiles

region='[[-72, -1], [-72, 11], [-59, 11], [-59, -1]]'

## build the URL and name the object (so that when it's unzipped we know what it is!)

path =mod09a.getDownloadUrl({

'name': output,  # name the file (otherwise it will be a uninterpretable hash)

'scale': 926,                               # resolution in meters

'crs': 'EPSG:4326',                         # MODIS sinusoidal

'region': strregion                  # region defined above

});

## download with wget

wget.download(path)

###  Script to compile the monthly cloud data from earth engine into a netcdf file for further processing

library(raster)

library(doMC)

library(multicore)

library(foreach)

#library(doMPI)

registerDoMC(4)

#beginCluster(4)

wd="~/acrobates/adamw/projects/cloud"

setwd(wd)

tempdir="tmp"

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

##  Get list of available files

df=data.frame(path=list.files("/mnt/data2/projects/cloud/mod09",pattern="*.tif$",full=T,recur=T),stringsAsFactors=F)

df[,c("region","year","month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]"))[,c(1,2,3)]

df$date=as.Date(paste(df$year,"_",df$month,"_15",sep=""),"%Y_%m_%d")

## add stats to test for missing data

addstats=F

if(addstats){

    df[,c("max","min","mean","sd")]=do.call(rbind.data.frame,mclapply(1:nrow(df),function(i) as.numeric(sub("^.*[=]","",grep("STATISTICS",system(paste("gdalinfo -stats",df$path[i]),inter=T),value=T)))))

    table(df$sd==0)

}

## subset to testtiles?

#df=df[df$region%in%testtiles,]

#df=df[df$month==1,]

table(df$year,df$month)

writeLines(paste("Tiling options will produce",nrow(tiles),"tiles and ",nrow(jobs),"tile-months.  Current todo list is ",length(todo)))

## drop some if not complete

#df=df[df$month%in%1:9&df$year%in%c(2001:2012),]

rerun=F  # set to true to recalculate all dates even if file already exists

## Loop over existing months to build composite netcdf files

foreach(date=unique(df$date)) %dopar% {

    ## get date

  print(date)

  ## Define output and check if it already exists

  vrtfile=paste(tempdir,"/mod09_",date,".vrt",sep="")

  ncfile=paste(tempdir,"/mod09_",date,".nc",sep="")

  tffile=paste(tempdir,"/mod09_",date,".tif",sep="")

  if(!rerun&file.exists(ncfile)) return(NA)

  ## merge regions to a new netcdf file

#  system(paste("gdal_merge.py -o ",tffile," -init -32768  -n -32768.000 -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))

  system(paste("gdalbuildvrt -overwrite -srcnodata -32768 -vrtnodata -32768 ",vrtfile," ",paste(df$path[df$date==date],collapse=" ")))

  ## Warp to WGS84 grid and convert to netcdf

  ##ops="-t_srs 'EPSG:4326' -multi -r cubic -te -180 0 180 10 -tr 0.008333333333333 -0.008333333333333 -wo SOURCE_EXTRA=50" #-wo SAMPLE_GRID=YES -wo SAMPLE_STEPS=100

  ops="-t_srs 'EPSG:4326' -multi -r cubic -te -180 -90 180 90 -tr 0.008333333333333 -0.008333333333333  -wo SOURCE_EXTRA=50"

  system(paste("gdalwarp -overwrite ",ops," -et 0 -srcnodata -32768 -dstnodata -32768 ",vrtfile," ",tffile," -ot Int16"))

  system(paste("gdal_translate -of netCDF ",tffile," ",ncfile))

                                        #  system(paste("gdalwarp -overwrite ",ops," -srcnodata -32768 -dstnodata -32768 -of netCDF ",vrtfile," ",tffile," -ot Int16"))

  file.remove(tffile)

  setwd(wd)

  system(paste("ncecat -O -u time ",ncfile," ",ncfile,sep=""))

## create temporary nc file with time information to append to MOD06 data

  cat(paste("

    netcdf time {

      dimensions:

        time = 1 ;

      variables:

        int time(time) ;

      time:units = \"days since 2000-01-01 00:00:00\" ;

      time:calendar = \"gregorian\";

      time:long_name = \"time of observation\"; 

    data:

      time=",as.integer(date-as.Date("2000-01-01")),";

    }"),file=paste(tempdir,"/",date,"_time.cdl",sep=""))

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

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

## add other attributes

  system(paste("ncrename -v Band1,CF ",ncfile,sep=""))

  system(paste("ncatted ",

               " -a units,CF,o,c,\"%\" ",

#               " -a valid_range,CF,o,b,\"0,100\" ",

               " -a scale_factor,CF,o,f,\"0.1\" ",

               " -a _FillValue,CF,o,f,\"-32768\" ",

               " -a long_name,CF,o,c,\"Cloud Frequency(%)\" ",

               " -a NETCDF_VARNAME,CF,o,c,\"Cloud Frequency(%)\" ",

               " -a title,global,o,c,\"Cloud Climatology from MOD09 Cloud Mask\" ",

               " -a institution,global,o,c,\"Jetz Lab, EEB, Yale University, New Haven, CT\" ",

               " -a source,global,o,c,\"Derived from MOD09GA Daily Data\" ",

               " -a comment,global,o,c,\"Developed by Adam M. Wilson (adam.wilson@yale.edu / http://adamwilson.us)\" ",

               ncfile,sep=""))

               ## add the fillvalue attribute back (without changing the actual values)

#system(paste("ncatted -a _FillValue,CF,o,b,-32768 ",ncfile,sep=""))

if(as.numeric(system(paste("cdo -s ntime ",ncfile),intern=T))<1) {

  print(paste(ncfile," has no time, deleting"))

  file.remove(ncfile)

}

  print(paste(basename(ncfile)," Finished"))

}

### merge all the tiles to a single global composite

#system(paste("ncdump -h ",list.files(tempdir,pattern="mod09.*.nc$",full=T)[10]))

file.remove("tmp/mod09_2000-01-15.nc")

system(paste("cdo -O  mergetime -setrtomiss,-32768,-1 ",paste(list.files(tempdir,pattern="mod09.*.nc$",full=T),collapse=" ")," data/cloud_monthly.nc"))

#  Overall mean

system(paste("cdo -O  timmean data/cloud_monthly.nc  data/cloud_mean.nc"))

### generate the monthly mean and sd

#system(paste("cdo -P 10 -O merge -ymonmean data/mod09.nc -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim.nc"))

system(paste("cdo  -f nc4c -O -ymonmean data/cloud_monthly.nc data/cloud_ymonmean.nc"))

## Seasonal Means

system(paste("cdo  -f nc4c -O -yseasmean data/cloud_monthly.nc data/cloud_yseasmean.nc"))

system(paste("cdo  -f nc4c -O -yseasstd data/cloud_monthly.nc data/cloud_yseasstd.nc"))

## standard deviations, had to break to limit memory usage

system(paste("cdo  -f nc4c -O -chname,CF,CF_sd -ymonstd -selmon,1,2,3,4,5,6 data/cloud_monthly.nc data/cloud_ymonsd_1-6.nc"))

system(paste("cdo  -f nc4c -O -chname,CF,CF_sd -ymonstd -selmon,7,8,9,10,11,12 data/cloud_monthly.nc data/cloud_ymonsd_7-12.nc"))

system(paste("cdo  -f nc4c -O mergetime  data/cloud_ymonsd_1-6.nc  data/cloud_ymonsd_7-12.nc data/cloud_ymonstd.nc"))

#if(!file.exists("data/mod09_metrics.nc")) {

    system("cdo -f nc4c  timmin data/cloud_ymonmean.nc data/cloud_min.nc")

    system("cdo -f nc4c  timmax data/cloud_ymonmean.nc data/cloud_max.nc")

    system("cdo -f nc4c -chname,CF,CFsd -timstd data/cloud_ymonmean.nc data/cloud_std.nc")

#    system("cdo -f nc2 merge data/mod09_std.nc data/mod09_min.nc data/cloud_max.nc data/cloud_metrics.nc")

#    system("cdo merge -chname,CF,CFmin -timmin data/cloud_ymonmean.nc -chname,CF,CFmax -timmax data/cloud_ymonmean.nc  -chname,CF,CFsd -timstd data/cloud_ymonmean.nc  data/cloud_metrics.nc")

#}

# Regressions through time by season

s=c("DJF","MAM","JJA","SON")

system(paste("cdo  -f nc4c -O regres -selseas,",s[1]," data/cloud_monthly.nc data/slope_",s[1],".nc",sep=""))

system(paste("cdo  -f nc4c -O regres -selseas,",s[2]," data/cloud_monthly.nc data/slope_",s[2],".nc",sep=""))

system(paste("cdo  -f nc4c -O regres -selseas,",s[3]," data/cloud_monthly.nc data/slope_",s[3],".nc",sep=""))

system(paste("cdo  -f nc4c -O regres -selseas,",s[4]," data/cloud_monthly.nc data/slope_",s[4],".nc",sep=""))

## Daily animations

regs=list(

    Venezuela=extent(c(-69,-59,0,7)),

    Cascades=extent(c(-122.8,-118,44.9,47)),

    Hawaii=extent(c(-156.5,-154,18.75,20.5)),

    Boliva=extent(c(-71,-63,-20,-15)),

    CFR=extent(c(17.75,22.5,-34.8,-32.6)),

    Madagascar=extent(c(46,52,-17,-12))

    )

r=1

system(paste("cdo  -f nc4c -O inttime,2012-01-15,12:00:00,1day  -sellonlatbox,",

             paste(regs[[r]]@xmin,regs[[r]]@xmax,regs[[r]]@ymin,regs[[r]]@ymax,sep=","),

             "  data/cloud_monthly.nc data/daily_",names(regs[r]),".nc",sep=""))

### Long term summaries

seasconc <- function(x,return.Pc=T,return.thetat=F) {

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

          ## Precipitation Concentration function

          ## This function calculates Precipitation Concentration based on Markham's (1970) technique as described in Schulze (1997)

          ## South Africa Atlas of Agrohydology and Climatology - R E Schulze, M Maharaj, S D Lynch, B J Howe, and B Melvile-Thomson

          ## Pages 37-38

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

          ## x is a vector of precipitation quantities - the mean for each factor in "months" will be taken,

          ## so it does not matter if the data are daily or monthly, as long as the "months" factor correctly

          ## identifies them into 12 monthly bins, collapse indicates whether the data are already summarized as monthly means.

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

          theta=seq(30,360,30)*(pi/180)                                       # set up angles for each month & convert to radians

                  if(sum(is.na(x))==12) { return(cbind(Pc=NA,thetat=NA)) ; stop}

                  if(return.Pc) {

                              rt=sqrt(sum(x * cos(theta))^2 + sum(x * sin(theta))^2)    # the magnitude of the summation

                                        Pc=as.integer(round((rt/sum(x))*100))}

                  if(return.thetat){

                              s1=sum(x*sin(theta),na.rm=T); s2=sum(x*cos(theta),na.rm=T)

                                        if(s1>=0 & s2>=0)  {thetat=abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                                        if(s1>0 & s2<0)  {thetat=180-abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                                        if(s1<0 & s2<0)  {thetat=180+abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                                        if(s1<0 & s2>0)  {thetat=360-abs((180/pi)*(atan(sum(x*sin(theta),na.rm=T)/sum(x*cos(theta),na.rm=T))))}

                             thetat=as.integer(round(thetat))

                            }

                  if(return.thetat&return.Pc) return(c(conc=Pc,theta=thetat))

                  if(return.Pc)          return(Pc)

                  if(return.thetat)  return(thetat)

        }

## read in monthly dataset

mod09=brick("data/cloud_ymonmean.nc",varname="CF")

plot(mod09[1])

mod09_seas=calc(mod09,seasconc,return.Pc=T,return.thetat=F,overwrite=T,filename="data/mod09_seas.nc",NAflag=255,datatype="INT1U")

mod09_seas2=calc(mod09,seasconc,return.Pc=F,return.thetat=T,overwrite=T,filename="data/mod09_seas_theta.nc",datatype="INT1U")

plot(mod09_seas)

library(doMC)

library(foreach)

registerDoMC(4)

wd="~/acrobates/adamw/projects/cloud"

setwd(wd)

### Build Tiles

## bin sizes in degrees

ybin=30

xbin=30

tiles=expand.grid(ulx=seq(-180,180-xbin,by=xbin),uly=seq(90,-90+ybin,by=-ybin))

tiles$h=factor(tiles$ulx,labels=paste("h",sprintf("%02d",1:length(unique(tiles$ulx))),sep=""))

tiles$v=factor(tiles$uly,labels=paste("v",sprintf("%02d",1:length(unique(tiles$uly))),sep=""))

tiles$tile=paste(tiles$h,tiles$v,sep="")

tiles$urx=tiles$ulx+xbin

tiles$ury=tiles$uly

tiles$lrx=tiles$ulx+xbin

tiles$lry=tiles$uly-ybin

tiles$llx=tiles$ulx

tiles$lly=tiles$uly-ybin

tiles$cy=(tiles$uly+tiles$lry)/2

tiles$cx=(tiles$ulx+tiles$urx)/2

tiles=tiles[,c("tile","h","v","ulx","uly","urx","ury","lrx","lry","llx","lly","cx","cy")]

jobs=expand.grid(tile=tiles$tile,year=2000:2012,month=1:12)

jobs[,c("ulx","uly","urx","ury","lrx","lry","llx","lly")]=tiles[match(jobs$tile,tiles$tile),c("ulx","uly","urx","ury","lrx","lry","llx","lly")]

## drop Janurary 2000 from list (pre-modis)

jobs=jobs[!(jobs$year==2000&jobs$month==1),]

## Run the python downloading script

#system("~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin -159 20 -154.5 18.5 -year 2001 -month 6 -region test")   

i=1

todo=1:nrow(jobs)

##  Get list of available files

df=data.frame(path=list.files("/mnt/data2/projects/cloud/mod09",pattern="*.tif$",full=T,recur=T),stringsAsFactors=F)

df[,c("region","year","month")]=do.call(rbind,strsplit(basename(df$path),"_|[.]"))[,c(1,2,3)]

df$date=as.Date(paste(df$year,"_",df$month,"_15",sep=""),"%Y_%m_%d")

table(df$year,df$month)

checkcomplete=T

if(checkcomplete&exists("df")){  #if desired (and "df" exists from below) drop complete date-tiles

todo=which(!paste(jobs$tile,jobs$year,jobs$month)%in%paste(df$region,df$year,df$month))

}

writeLines(paste("Tiling options will produce",nrow(tiles),"tiles and ",nrow(jobs),"tile-months.  Current todo list is ",length(todo)))

t=foreach(i=todo,.inorder=FALSE,.verbose=F) %dopar%{

    system(paste("python ~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin ",

                      jobs$ulx[i]," ",jobs$uly[i]," ",jobs$urx[i]," ",jobs$ury[i]," ",jobs$lrx[i]," ",jobs$lry[i]," ",jobs$llx[i]," ",jobs$lly[i]," ",

                      "  -year ",jobs$year[i]," -month ",jobs$month[i]," -region ",jobs$tile[i],sep=""))

     }