### Figures and tables for MOD09 Cloud Manuscript

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

## libraries

library(rasterVis)

library(latticeExtra)

library(xtable)

library(reshape)

library(caTools)

## read in data

cld=read.csv("data/NDP026D/cld.csv")

cldm=read.csv("data/NDP026D/cldm.csv")

cldy=read.csv("data/NDP026D/cldy.csv")

clda=read.csv("data/NDP026D/clda.csv")

st=read.csv("data/NDP026D/stations.csv")

## add lulc factor information

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)

## month factors

cld$month2=factor(cld$month,labels=month.name)

cldm$month2=factor(cldm$month,labels=month.name)

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

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

##make spatial object

cldms=cldm

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

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

##make spatial object

cldys=cldy

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

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

#### Evaluate MOD35 Cloud data

mod09=brick("data/mod09.nc")

mod09c=brick("data/mod09_clim_mean.nc",varname="CF");names(mod09c)=month.name

mod09a=brick("data/mod09_clim_mac.nc",varname="CF_annual")#;names(mod09c)=month.name

## derivatives

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

  system("cdo -chname,CF,CFmin -timmin data/mod09_clim_mean.nc data/mod09_min.nc")

  system("cdo -chname,CF,CFmax -timmax data/mod09_clim_mean.nc data/mod09_max.nc")

  system("cdo -chname,CF,CFsd -timstd data/mod09_clim_mean.nc data/mod09_std.nc")

  system("cdo -f nc2 merge data/mod09_std.nc data/mod09_min.nc data/mod09_max.nc data/mod09_metrics.nc") 

}

mod09min=raster("data/mod09_metrics.nc",varname="CFmin")

mod09max=raster("data/mod09_metrics.nc",varname="CFmax")

mod09sd=raster("data/mod09_metrics.nc",varname="CFsd")

mod09mean=raster("data/mod09_clim_mac.nc")

names(mod09d)=c("Mean","Minimum","Maximum","Standard Deviation")

plot(mod09a,layers=1,margin=F,maxpixels=100)

## calculated differences

cldm$dif=cldm$mod09-cldm$cld

clda$dif=clda$mod09-clda$cld

## read in global coasts for nice plotting

library(maptools)

data(wrld_simpl)

coast <- unionSpatialPolygons(wrld_simpl, rep("land",nrow(wrld_simpl)), threshold=5)

coast=as(coast,"SpatialLines")

## Figures

n=100

  at=seq(0,100,length=n)

colr=colorRampPalette(c("black","green","red"))

cols=colr(n)

pdf("output/validation.pdf",width=11,height=8.5)

## 4-panel maps

#- Annual average

levelplot(mod09a,col.regions=colr(100),cuts=100,at=seq(0,100,len=100),colorkey=list(space="bottom",adj=1),

  margin=F,maxpixels=1e6,ylab="Latitude",xlab="Longitude",useRaster=T)+

  layer(sp.lines(coast,col="black"),under=F)

#- Monthly minimum

#- Monthly maximum

#- STDEV or Min-Max

### maps of mod09 and NDP

## map of stations

p_mac=levelplot(mod09a,col.regions=colr(100),cuts=100,at=seq(0,100,len=100),colorkey=list(space="bottom",adj=1),

  margin=F,maxpixels=1e6,ylab="Latitude",xlab="Longitude",useRaster=T)+

  layer(panel.xyplot(lon,lat,pch=16,cex=.3,col="black"),data=data.frame(coordinates(st)))+

  layer(sp.lines(coast,col="black"),under=F)

p_mace=xyplot(lat~dif,data=cldm[cldm$lat>-60,],panel=function(x,y,subscripts=T){

  x2=x[order(y)]

  y2=y[order(y)]

  win=8000

  Q50=runquantile(x2,win,probs=c(.25,.5,.75))

  ## polygon

  panel.polygon(c(Q50[,1],rev(Q50[,3])),c(y2,rev(y2)),type="l",col=grey(.8))

### hist

  n=150

  xbins=seq(-70,50,len=n)

  ybins=seq(-60,90,len=n)

  tb=melt(as.matrix(table(

    x=cut(x,xbins,labels=xbins[-1]),

    y=cut(y,ybins,labels=ybins[-1]))))

  qat=unique(tb$value)

  print(qat)

  panel.levelplot(tb$x,tb$y,tb$value,at=qat,col.regions=c("transparent",hmcols(length(qat))),subscripts=1:nrow(tb))

###

#  panel.xyplot(x,y,pch=16,cex=.1,col=)

#  cuts=cut(y,lats,labels=lats[-10])

#  qs=do.call(rbind,tapply(x,cuts,quantile,c(.25,.50,.75),na.rm=T))

  colnames(qs)=c("Q25","Q50","Q75")

  panel.lines(Q50[,1],y2,type="l",col=grey(.5))

  panel.lines(Q50[,2],y2,type="l",col="black")

  panel.lines(Q50[,3],y2,type="l",col=grey(.5))

},asp=1,xlab="Difference (MOD09-Observed)")+layer(panel.abline(v=0,lty="dashed",col="red"))

print(p_mac,position=c(0,0,.75,1),more=T)

print(p_mace,position=c(0.75,0,1,1))

p1=c("MODIS Cloud Frequency and NDP-026D Validation Stations"=p_mac,"Difference (MOD09-NDP026D)"=p_mace,x.same=F,y.same=T,layout=c(2,1))

resizePanels(p1,w=c(.75,.25))

quantile(cldm$dif,seq(0,1,len=6),na.rm=T)

at=c(-70,-50,-25,-10,-5,0,5,10,25,50,70)

bwr=colorRampPalette(c("blue","grey","red"))

xyplot(lat~lon|month2,data=cldm,groups=cut(cldm$dif,at),

       par.settings=list(superpose.symbol=list(col=bwr(length(at)-1))),pch=16,cex=.25,

       auto.key=list(space="right",title="Difference\n(MOD09-NDP026D)",cex.title=1),asp=1,

       main="NDP-026D Cloud Climatology Stations",ylab="Latitude",xlab="Longitude")+

  layer(sp.lines(coast,col="black",lwd=.1),under=F)

### heatmap of mod09 vs. NDP for all months

hmcols=colorRampPalette(c("grey","blue","red"))

#hmcols=colorRampPalette(c(grey(.8),grey(.3),grey(.2)))

tr=c(0,120)

colkey <- draw.colorkey(list(col = hmcols(tr[2]), at = tr[1]:tr[2],height=.25))

xyplot(cld~mod09,data=cld[cld$Nobs>10,],panel=function(x,y,subscripts){

  n=150

  bins=seq(0,100,len=n)

  tb=melt(as.matrix(table(

    x=cut(x,bins,labels=bins[-1]),

    y=cut(y,bins,labels=bins[-1]))))

  qat=tr[1]:tr[2]#unique(tb$value)

  print(qat)

  panel.levelplot(tb$x,tb$y,tb$value,at=qat,col.regions=c("transparent",hmcols(length(qat))),subscripts=subscripts)

  },asp=1,scales=list(at=seq(0,100,len=6)),ylab="NDP Mean Cloud Amount (%)",xlab="MOD09 Cloud Frequency (%)",

       legend= list(right = list(fun = colkey,title="Station Count")))+

  layer(panel.abline(0,1,col="black",lwd=2))

#  layer(panel.ablineq(lm(y ~ x), r.sq = TRUE,at = 0.6,pos=1, offset=22,digits=2,col="blue"), style = 1)

xyplot(cld~mod09|month2,data=cld[cld$Nobs>50,],panel=function(x,y,subscripts){

  n=50

  bins=seq(0,100,len=n)

  tb=melt(as.matrix(table(

    x=cut(x,bins,labels=bins[-1]),

    y=cut(y,bins,labels=bins[-1]))))

  qat=unique(tb$value)

  qat=0:78

  qat=tr[1]:tr[2]#unique(tb$value)

  panel.levelplot(tb$x,tb$y,tb$value,at=qat,col.regions=c("transparent",hmcols(length(qat))),subscripts=1:nrow(tb))

  panel.abline(0,1,col="black",lwd=2)

#  panel.ablineq(lm(y ~ x), r.sq = TRUE,at = 0.6,pos=1, offset=0,digits=2,col="blue")

  panel.text(70,10,bquote(paste(R^2,"=",.(round(summary(lm(y ~ x))$r.squared,2)))),cex=1.2)

},asp=1,scales=list(at=seq(0,100,len=6),useRaster=T,colorkey=list(width=.5,title="Number of Stations")),

          ylab="NDP Mean Cloud Amount (%)",xlab="MOD09 Cloud Frequency (%)",

              legend= list(right = list(fun = colkey)))+ layer(panel.abline(0,1,col="black",lwd=2))

## Monthly Climatologies

for(i in 1:2){

 p1=xyplot(cld~mod09|month2,data=cldm,cex=.2,pch=16,subscripts=T,ylab="NDP Mean Cloud Amount",xlab="MOD09 Cloud Frequency (%)")+

  layer(panel.lines(1:100,predict(lm(y~x),newdata=data.frame(x=1:100)),col="green"))+

  layer(panel.lines(1:100,predict(lm(y~x+I(x^2)),newdata=data.frame(x=1:100)),col="blue"))+

  layer(panel.abline(0,1,col="red"))

    if(i==2){

     p1=p1+layer(panel.segments(mod09[subscripts],cld[subscripts]-cldsd[subscripts],mod09[subscripts],cld[subscripts]+cldsd[subscripts],subscripts=subscripts,col="grey"),data=cldm,under=T,magicdots=T)

     p1=p1+layer(panel.segments(mod09[subscripts]-mod09sd[subscripts],cld[subscripts],mod09[subscripts]+mod09sd[subscripts],cld[subscripts],subscripts=subscripts,col="grey"),data=cldm,under=T,magicdots=T) 

       }

print(p1)

}

bwplot(lulcc~dif,data=cldm,horiz=T,xlab="Difference (MOD09-Observed)",varwidth=T,notch=T)+layer(panel.abline(v=0))

dev.off()

summary(lm(cld~mod09,data=cld))

## explore validation error

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

## Table of RMSE's by lulc by month

lulcrmsel=ddply(cldm,c("month","lulc"),function(x) c(count=nrow(x),rmse=sqrt(mean((x$mod09-x$cld)^2,na.rm=T))))

lulcrmsel=lulcrmsel[!is.na(lulcrmsel$lulc),]

lulcrmsel$lulcc=as.factor(IGBP$class[match(lulcrmsel$lulc,IGBP$ID)])

lulcrmse=cast(lulcrmsel,lulcc~month,value="rmse")

lulcrmse

print(xtable(lulcrmse,digits=1),"html")

levelplot(rmse~lulc*month,data=lulcrmsel,col.regions=heat.colors(20))

### Linear models

summary(lm(dif~as.factor(lulc)+lat+month2,data=cldm))

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

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

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

## read it back in

colnames(mod09st)=c(names(mod09)[-1],"id")

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

## function to get modal lulc value

Mode <- function(x) {

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

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

             cld=mean(x$cld[x$Nobs>50],na.rm=T),

             cldsd=sd(x$cld[x$Nobs>50],na.rm=T))}))

             cld=mean(x$cld[x$Nobs>50]/100,na.rm=T),

             cldsd=sd(x$cld[x$Nobs>50]/100,na.rm=T))}))

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

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

#!/usr/bin/env python

import sys

import logging

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

def Usage():

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

    sys.exit( 1 )

ulx = float(sys.argv[2])

uly = float(sys.argv[3])

lrx = float(sys.argv[4])

lry = float(sys.argv[5])

year = int(sys.argv[7])

month = int(sys.argv[9])

regionname = str(sys.argv[11])

#```

#ulx=-159

#uly=20

#lrx=-154.5

#lry=18.5

#year=2001

#month=6

#```

## Define output filename

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

unzippedfilename=output+".mod09.tif"

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

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

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

# Processing Function

# MOD09 internal cloud flag for this year-month

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

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=[ulx,lry], [ulx, uly], [lrx, uly], [lrx, lry]  #h11v08

strregion=str(list(region))

#data.getInfo()

## print a status update

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

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

#```

# TODO:  

#  use MODIS projection

path =mod09a.getDownloadUrl({

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

        'scale': 926,                              # resolution in meters

        'crs': 'EPSG:4326',                         #  projection

        'region': strregion                        # region defined above

        });

#if(not(os.path.exists(output+".tif"))):

    # download with wget

print("Downloading "+output) 

wget.download(path)

#call(["wget"+path,shell=T])

print("Unzipping "+output)

zipstatus=call("unzip "+output+".zip",shell=True)

if(zipstatus==9):

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

#        print("ERROR: "+output+" unzip-able")

#        os.remove(output+".zip")

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

os.remove(output+".zip")

print(output+' Finished!')

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

## Load list of tiles

tiles=read.table("tile_lat_long_10d.txt",header=T)

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

jobs[,c("ULX","ULY","LRX","LRY")]=tiles[match(jobs$tile,tiles$Tile),c("ULX","ULY","LRX","LRY")]

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

testtiles=c("h02v07","h02v06","h02v08","h03v07","h03v06")

todo=which(jobs$tile%in%testtiles)

#todo=todo[1:3]

#todo=1:nrow(jobs)

lapply(todo,function(i)

         system(paste("~/acrobates/adamw/projects/environmental-layers/climate/procedures/ee.MOD09.py -projwin ",jobs$ULX[i]," ",jobs$ULY[i]," ",jobs$LRX[i]," ",jobs$LRY[i],

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

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

## subset to testtiles?

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

table(df$year,df$month)

  system(paste("gdal_merge.py -o ",ncfile," -n -32768 -of netCDF -ot Int16 ",paste(df$path[df$date==date],collapse=" ")))

               " -a units,CF,o,c,\"Proportion Days Cloudy\" ",

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

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

               " -a long_name,CF,o,c,\"Proportion cloudy days (%)\" ",

               ncfile,sep=""))

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

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