Validation

----------------------

```{r,echo=F}

## working directory

setwd("/mnt/data/personal/adamw/projects/cloud/")

## libraries

library(rasterVis)

library(latticeExtra)

library(xtable)

library(texreg)

library(reshape)

library(caTools)

library(rgeos)

library(raster)

## read in global coasts for nice plotting

library(maptools)

library(rgdal)

```

```{r, loaddata,echo=F}

# read in data

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

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

## month factors

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

### Drop valitation station-months with fewer than 20 years of data for full record or less than 10 years for MODIS-era record

cldm$cld_all[cldm$cldn_all<20]=NA

cldm$cldsd_all[cldm$cldn_all<20]=NA

cldm$cld[cldm$cldn<10]=NA

cldm$cldsd[cldm$cldn<10]=NA

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

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

#coast=getRgshhsMap("/mnt/data/jetzlab/Data/environ/global/gshhg/gshhs_h.b", xlim = NULL, ylim = NULL, level = 4) 

land=readShapePoly("/mnt/data/jetzlab/Data/environ/global/gshhg/GSHHS_shp/c/GSHHS_c_L1.shp",force_ring=TRUE)

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

CP <- as(extent(-180, 180, -60, 84), "SpatialPolygons")

proj4string(CP) <- CRS(proj4string(land))

coast=as(land[land$area>50,],"SpatialLines")

## Clip the map

land <- gIntersection(land, CP, byid=F)

coast <- gIntersection(coast, CP, byid=F)

#### Evaluate MOD35 Cloud data

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

```

```{r}

## Figures

n=100

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

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

cols=colr(n)

## set plotting parameters

my.theme = trellis.par.get()

my.theme$strip.background=list(col="transparent")

trellis.par.set(my.theme)

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

png("output/CF_Figures_%03d.png",width=5000,height=4000,res=600,pointsize=42,bg="white")

## set plotting parameters

my.theme = trellis.par.get()

my.theme$strip.background=list(col="transparent")

trellis.par.set(my.theme)

greg=list(ylim=c(-60,84),xlim=c(-180,180))

bgr=function(x,n=100,br=0,c1=c("darkblue","blue","grey"),c2=c("grey","red","purple")){

    at=unique(c(seq(min(x,na.rm=T),max(x,na.rm=T),len=n)))

    bg=colorRampPalette(c1)

    gr=colorRampPalette(c2)

    return(list(at=at,col=c(bg(sum(at<br)),gr(sum(at>=br)))))

}

cldm$resid=NA

# get residuals of simple linear model

cldm$resid[!is.na(cldm$cld_all)&!is.na(cldm$mod09)]=residuals(lm(mod09~cld_all,data=cldm))

colat=bgr(cldm$resid)

phist=histogram(cldm$resid,breaks=colat$at,border=NA,col=colat$col,xlim=c(-30,30),type="count",xlab="MODCF Residuals")#,seq(0,1,len=6),na.rm=T)

pmap=xyplot(lat~lon|month2,data=cldm,groups=cut(cldm$resid,rev(colat$at)),

       par.settings=list(superpose.symbol=list(col=colat$col)),pch=16,cex=.25,

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

       ylab="Latitude",xlab="Longitude")+

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

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

print(pmap,position=c(0,0,1,.78))

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

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

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

tr=c(0,66)

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

xyplot(cld_all~mod09|month2,data=cldm,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)

  print(max(qat))

  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.abline(lm(y ~ x),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)

},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="MODCF Cloud Frequency (%)",

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

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

dev.off()

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

### Regional Comparisons

## Compare with worldclim and NPP

#wc=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/prec_",1:12,".bil",sep="")))

wc_map=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/bio_12.bil",sep="")))

wc_dem=stack(as.list(paste("/mnt/data/jetzlab/Data/environ/global/worldclim/alt.bil",sep="")))

regs=list(

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

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

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

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

  #reg2=extent(c(-81,-70,-4,10))

  )

## read in GEWEX 1-degree data

gewex=mean(brick("data/gewex/CA_PATMOSX_NOAA.nc",varname="a_CA"))

names(gewex)="PATMOS-x GEWEX AVHRR"

## calculate 1-degree means of MODCF data

#MOD_gewex=gewex

#MOD_gewex@data@values=1:length(MOD_gewex@data@values)

#MOD_gewex2=zonal(mod09a,MOD_gewex,fun='mean')

png("output/Resolution_Figures_%03d.png",width=5500,height=4000,res=600,pointsize=36,bg="white")

trellis.par.set(my.theme)

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

r="Venezuela"

# ylab.right = "Cloud Frequency (%)",par.settings = list(layout.widths = list(axis.key.padding = 0.1,axis.left=0.6,ylab.right = 3,right.padding=2)),

pars=list(layout.heights=list(key.bottom=2,key.top=1),layout.widths = list(axis.key.padding = 3,axis.left=0.6))

p1=levelplot(crop(mod09a,regs[[r]]),col.regions=grey(seq(0,1,len=100)),at=seq(45,100,len=99),

    colorkey=list(space="top",width=1,height=.75,labels=list(labels=c(50,75,100),at=c(50,75,100))),

    cuts=99,margin=F,max.pixels=1e6,par.settings = pars)

p2=levelplot(crop(gewex,regs[[r]]),col.regions=grey(seq(0,1,len=100)),at=seq(.45,1,len=99),cuts=99,margin=F,max.pixels=1e6,

    colorkey=list(space="top",width=1,height=.75,labels=list(labels=c(50,75,100),at=c(.50,.75,1))),

    par.settings = pars)

tmap=crop(wc_map,regs[[r]])

p3=levelplot(tmap,col.regions=grey(seq(0,1,len=100)),cuts=100,at=seq(tmap@data@min,tmap@data@max,len=100),margin=F,maxpixels=1e6,

    colorkey=list(space="bottom",height=.75,width=1),xlab="",ylab="",main=names(regs)[r],useRaster=T,

    par.settings = pars)

p4=levelplot(crop(wc_dem,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6,

    colorkey=list(space="bottom",height=.75,width=1),

    par.settings = pars)#,labels=list(labels=c(1000,4000),at=c(1000,4000))))

print(c("MODCF (%)"=p1,"PATMOS-x GEWEX (%)"=p2,"WorldClim Precip (mm)"=p3,"Elevation (m)"=p4,x.same=T,y.same=T,merge.legends=T,layout=c(2,2)))

dev.off()

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

### Some stats for paper

## number of stations retained

length(unique(cldm$StaID[!is.na(cldm$cld_all)]))

length(unique(cldm$StaID[!is.na(cldm$cld)]))

# approximate size of mod09ga archive - get total size for one day from the USGS website

size=scan("http://e4ftl01.cr.usgs.gov/MOLT/MOD09GA.005/2000.04.30/",what="char")

## extract all filesizes in MB (all the HDFs) and sum them and covert to TB for the length of the full record

sum(as.numeric(sub("M","",grep("[0-9]*M$",size,value=T))))/1024/1024*as.numeric(as.Date("2013-12-31")-as.Date("2000-02-24"))

## seasonal variability

cellStats(mod09sd,"mean")

## Validation table construction

quantile(cldm$difm,na.rm=T)

summary(lm(cld_all~mod09+lat,data=cldm))

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

### summary by biome

biomep$id=1:nrow(biomep)

biomepl=melt(biomep@data,id.vars=c("id","code","biome","realm"))

colnames(biomepl)[grep("variable",colnames(biomepl))]="month"

biomepl$month=factor(biomepl$month,ordered=T,levels=month.name)

biomepl$realm=factor(biomepl$realm,ordered=T,levels=c("Antarctic","Australasia","Oceania", "IndoMalay", "Neotropics","Palearctic","Nearctic" ))

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

png("output/Biome_Figures_%03d.png",width=5500,height=4000,res=600,pointsize=36,bg="white")

trellis.par.set(my.theme)

#[biomepl$id%in%sample(biomep$id,10000),]

p1=useOuterStrips(xyplot(value~month|realm+biome,groups=id,data=biomepl,panel=function(x,y,groups = groups, subscripts = subscripts){

    panel.xyplot(x,y,col=grey(0.6,0.2),type="l",lwd=.5,groups=groups,subscripts=subscripts)

    panel.smoother(y ~ s(x), method = "gam",lwd=2,subscripts=subscripts,n=24)

},scales=list(y=list(at=c(0,100),lim=c(-20,120),cex=.75,alternating=2,tck=c(0,1)),x=list(at=c(1,7),rot=45,alternating=1)),ylab="Biome",xlab.top="Geographic Realm",ylab.right="MODCF (%)", xlab="Month"),

    strip=strip.custom(par.strip.text = list(cex = .7)),strip.left=strip.custom(horizontal=TRUE,par.strip.text = list(cex = .75)))

p1$par.settings$layout.widths$strip.left[1]=13

p1$par.strip.text$lines=.65

print(p1)

dev.off()

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

## assess temporal stability

## spatialy subset data to stations at least 10km apart

st2=remove.duplicates(st,zero=10)

## Subset data

## drop missing observations

cldm.t=cldm[!is.na(cldm$cld_all)&!is.na(cldm$mod09)&!is.na(cldm$biome),]

cldm.t=cldm.t[cldm.t$lat>=-60,]

#  make sure all stations have all mod09 data

stdrop=names(which(tapply(cldm.t$month,cldm.t$StaID,length)!=12))

cldm.t=cldm.t[!cldm.t$StaID%in%stdrop,]

# Keep only stations at least 10km apart 

cldm.t=cldm.t[cldm.t$StaID%in%st2$id,]

## Subset to only some months, if desired

#cldm.t=cldm.t[cldm.t$month%in%1:3,]

## Select Knots

knots=spsample(land,500,type="regular")

                                        #  reshape data

m.cld=cast(cldm.t,StaID+lat+lon+biome~month,value="cld_all");colnames(m.cld)[-(1:4)]=paste("cld.",colnames(m.cld)[-(1:4)],sep="")

m.mod09=cast(cldm.t,StaID~month,value="mod09");colnames(m.mod09)[-1]=paste("mod09.",colnames(m.mod09)[-1],sep="")

mdata=cbind(m.cld,m.mod09)

## cast to 

coords <- as.matrix(m.cld[,c("lon","lat")])#as.matrix(ne.temp[,c("UTMX", "UTMY")]/1000)

max.d <- max(iDist(coords))

##make symbolic model formula statement for each month

mods <- lapply(paste(paste(paste("cld.",1:N.t,sep=''),paste("mod09.",1:N.t,sep=''),sep='~'),"",sep=""), as.formula)

tlm=model.matrix(lm(mods[[1]],data=mdata))

N.t <- ncol(m.mod09)-1 ##number of months

n <- nrow(m.cld) ##number of observation per months

p <- ncol(tlm) #number of regression parameters in each month

starting <- list("beta"=rep(0,N.t*p), "phi"=rep(3/(0.5*max.d), N.t),

                 "sigma.sq"=rep(2,N.t), "tau.sq"=rep(1, N.t),

                 "sigma.eta"=diag(rep(0.01, p)))

tuning <- list("phi"=rep(5, N.t))

priors <- list("beta.0.Norm"=list(rep(0,p), diag(1000,p)),

               "phi.Unif"=list(rep(3/(0.9*max.d), N.t), rep(3/(0.05*max.d), N.t)),

               "sigma.sq.IG"=list(rep(2,N.t), rep(10,N.t)),

               "tau.sq.IG"=list(rep(2,N.t), rep(5,N.t)),

               "sigma.eta.IW"=list(2, diag(0.001,p)))

cov.model <- "exponential"

## Run the model

n.samples <- 500

m.1=spDynLM(mods,data=mdata,coords=coords,knots=coordinates(knots),n.samples=n.samples,starting=starting,tuning=tuning,priors=priors,cov.model=cov.model,get.fitted=T,n.report=25)

save(m.1,file="output/m.1.Rdata")

## summarize

burn.in <- floor(0.75*n.samples)

quant <- function(x){quantile(x, prob=c(0.5, 0.025, 0.975))}

beta <- apply(m.1$p.beta.samples[burn.in:n.samples,], 2, quant)

beta.0 <- beta[,grep("Intercept", colnames(beta))]

beta.1 <- beta[,grep("mod09", colnames(beta))]

### Compare time periods

library(texreg)

 extract.lm <- function(model) {

     s <- summary(model)

     names <- rownames(s$coef)

     co <- s$coef[, 1]

     se <- s$coef[, 2]

     pval <- s$coef[, 4]

     rs <- s$r.squared

     n <- as.integer(nobs(model))

     rmse=sqrt(mean((residuals(s)^2)))

     gof <- c(rs, rmse, n)

     gof.names <- c("R-Squared","RMSE","n")

     tr <- createTexreg(coef.names = names, coef = co, se = se, 

                        pvalues = pval, gof.names = gof.names, gof = gof)

     return(tr)

 }

setMethod("extract", signature = className("lm", "stats"),definition = extract.lm)

forms=c("cld~mod09+month2+lat")

lm_all=lm(cld_all~mod09+lat,data=cldm[!is.na(cldm$cld),])

### Compare two time periods

lm_all1=lm(cld_all~mod09,data=cldm[!is.na(cldm$cld)&cldm$cldn_all>=10,])

lm_mod=lm(cld~mod09,data=cldm[cldm$cldn==10,])

mods=list("1970-2009"=lm_all1,"2000-2009"=lm_mod)

screenreg(mods,digits=2,single.row=T,custom.model.names=names(mods),custom.coef.names = c("Intercept", "MODCF"))

htmlreg(mods,file = "output/tempstab.doc",

        custom.model.names = names(mods),

        single.row = T, inline.css = FALSE,doctype = TRUE, html.tag = TRUE, head.tag = TRUE, body.tag = TRUE)

## assess latitude bias

cldm$abslat=abs(cldm$lat)

cldm$absdif=abs(cldm$difm)

abslm=lm(absdif~abslat*I(abslat^2),data=cldm[cldm$cldn_all>30,])

xyplot(absdif~abslat|month2,type=c("p","smooth"),data=cldm,cex=.25,pch=16)

plot(absdif~abslat,data=cldm[cldm$cldn_all>30,],cex=.25,pch=16)

lines(0:90,predict(abslm,newdata=data.frame(abslat=0:90),type="response"),col="red")

bf=anovaBF(dif~lulcc+month2,data=cldm[!is.na(cldm$dif)&!is.na(cldm$lulcc),])

ch=posterior(bf, iterations = 1000)

summary(bf)

plot(bf)

## explore validation error

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

## Table of RMSE's by lulc by month

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

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

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

lulctl=ddply(cldm,c("lulc"),function(x) c(count=nrow(x),mean=paste(round(mean(x$difm,na.rm=T),2)," (",round(sd(x$difm,na.rm=T),2),")",sep=""),rmse=round(sqrt(mean((x$difm)^2,na.rm=T)),2)))

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

    print(xtable(lulctl[order(lulctl$rmse),c("lulcc","count","mean","rmse")],digits=1),type="html",include.rownames=F,file="output/lulcc.doc",row.names=F)

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

lulcrmse

lulcrmse.q=round(do.call(rbind,apply(lulcrmse,1,function(x) data.frame(Min=min(x,na.rm=T),Mean=mean(x,na.rm=T),Max=max(x,na.rm=T),SD=sd(x,na.rm=T)))),1)#quantile,c(0.025,0.5,.975),na.rm=T)),1)

lulcrmse.q=lulcrmse.q[order(lulcrmse.q$Mean,decreasing=T),]

lulcrmse.q

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

bgyr=colorRampPalette(c("blue","green","yellow","red"))

levelplot(rmse~month*lulcc,data=lulcrmsel,col.regions=bgyr(1000),at=quantile(lulcrmsel$rmse,seq(0,1,len=100),na.rm=T))

### Linear models

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

mod09_mean=stack(list.files("data/mcd09tif/",pattern="MCD09_mean_[0-9]*[.]tif",full=T))

NAvalue(mod09_mean)=255

names(mod09_mean)=month.name

mod09_sd=stack(list.files("data/mcd09tif/",pattern="MCD09_sd_[0-9]*[.]tif",full=T))

NAvalue(mod09_sd)=255

names(mod09_sd)=month.name

beginCluster(12)

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

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

endCluster()

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

colnames(mod09stl)[grep("variable",colnames(mod09stl))]="month"

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

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

cldm$monthname=month.name[cldm$month]

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

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

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

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

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

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

writeOGR(st,dsn="data/validation/",layer="stations",driver="ESRI Shapefile",overwrite_layer=T)