## add tags for distribution

## MOD35

tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD35 Cloud Frequency for 2009 extracted from MOD09GA state_1km bits 0-1. The MOD35 bits encode four categories (with associated confidence that the pixel is actually clear): confidently clear (confidence > 0.99), probably clear (0.99 >= confidence > 0.95), probably cloudy (0.95 >= confidence > 0.66), and confidently cloudy (confidence <= 0.66).  Following the advice of the MODIS science team (Frey, 2010), we binned confidently clear and probably clear together as clear and the other two classes as cloudy.  The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days during 2009'",

  "TIFFTAG_DOCUMENTNAME='Collection 5 MOD35 Cloud Frequency'",

  "TIFFTAG_DATETIME='20090101'",

  "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")

system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py data/MOD35_2009.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))

## MOD09

tags=c("TIFFTAG_IMAGEDESCRIPTION='Collection 5 MOD09 Cloud Frequency for 2009 extracted from MOD09GA \'PGE11\' internal cloud mask algorithm (embedded in MOD09GA \'state_1km\' bit 10. The daily cloud mask time series were summarized to mean cloud frequency (CF) by calculating the proportion of cloudy days during 2009'",

  "TIFFTAG_DOCUMENTNAME='Collection 5 MOD09 Cloud Frequency'",

  "TIFFTAG_DATETIME='20090101'",

  "TIFFTAG_ARTIST='Adam M. Wilson (adam.wilson@yale.edu)'")

system(paste("/usr/local/src/gdal-1.10.0/swig/python/scripts/gdal_edit.py data/MOD09_2009.tif ",paste("-mo ",tags,sep="",collapse=" "),sep=""))

NAvalue(pp_bias)=255

NAvalue(lulc_bias)=255

## read in WWF biome data to summarize by biome

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

  system("wget -O data/teow.zip http://assets.worldwildlife.org/publications/15/files/original/official_teow.zip?1349272619")

  system("unzip data/teow.zip -d data/teow/")

   biome=readOGR("data/teow/wwf_terr_ecos.shp","wwf_terr_ecos")

   biome=biome[biome$BIOME<50,]

   biome2=gUnaryUnion(biome,id=biome$BIOME)

  ## create biome.csv using names in html file   

  biomeid=read.csv("data/teow/biome.csv",stringsAsFactors=F)

  biome2=SpatialPolygonsDataFrame(biome2,data=biomeid)

  writeOGR(biome2,"data/teow","biomes",driver="ESRI Shapefile")

}

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

biome2=extract(pp_bias,biome[biome$Biome==12,],df=T,fun=function(x) data.frame(mean=mean(x,na.rm=T),sd=sd(x,na.rm=T),prop=(sum(!is.na(x))/length(x))))

pat=c(seq(0,100,len=100),254)#seq(0,0.-5,len=2) #,seq(0.05,.1,len=50))

grayr2=colorRampPalette(c("grey","green","red"))#

#grayr2=colorRampPalette(grey(c(.75,.5,.25)))

levelplot(lulc_bias,col.regions=c(grayr2(100),"black"),at=pat,

          colorkey=T,margin=F,maxpixels=1e4)+layer(sp.lines(coast,lwd=.5))

histogram(lulc_bias)

##download 3 days of modis swath data:

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

#logging.basicConfig()

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

os.chdir('/home/adamw/acrobates/adamw/projects/cloud/data/mod09')

#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

MY_SERVICE_ACCOUNT = '205878743334-4mrtqgu0n5rnsv1vanrvv6atqk6vu8am@developer.gserviceaccount.com'

MY_PRIVATE_KEY_FILE = '/home/adamw/EarthEngine_Jeremy-privatekey.p12'

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

#// Date ranges

yearstart=2000

yearstop=2010

monthstart=1

monthstop=12

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

# Loop through months and get monthly % missing data

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

year=2001

month=1

for year in range(yearstart,yearstop+1): {

for month in range(monthstart,monthstop); {

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

## add to plot to confirm it's working

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

## define region for download

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

#  '[[-90, -1], [-90, 20], [-49, 20], [-49, -1]]'  // Northern S. America

#  '[[-180, -90], [-180, 90], [180, 90], [180, -90]]'  //global

#  '[[-180, -60], [-180, 90], [180, 90], [180, -60]]'  // Western Hemisphere

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

path =mod09a.getDownloadUrl({

  'name': 'mod09_'+str(year)+"_"+str(month),  # name the file (otherwise it will be a uninterpretable hash)

  'scale': 926,                               # resolution in meters

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

  'region': region

});

## download with wget

wget.download(path)

}}