Project

General

Profile

« Previous | Next » 

Revision 77d8eaf2

Added by Adam Wilson almost 11 years ago

updating figures

View differences:

climate/procedures/MOD09_CloudFigures.R
10 10 
library(texreg)
11 11 
library(reshape)
12 12 
library(caTools)
13 

  		




  
  13
  
    
library(rgeos)


  		




  14
  14
  
    

  		




  15
  15
  
    
## read in data


  		




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


  		




  ......




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


  		




  44
  44
  
    
mod09a=brick("data/cloud_mean.nc",varname="CF");names(mod09a)="Mean Annual Cloud Frequency"


  		




  45
  45
  
    

  		




  46
  
  
    
mod09min=brick("data/cloud_min.nc",varname="CFmin")


  		




  47
  
  
    
mod09max=brick("data/cloud_max.nc",varname="CFmax")


  		




  
  46
  
    
mod09min=brick("data/cloud_min.nc",varname="CF")


  		




  
  47
  
    
mod09max=brick("data/cloud_max.nc",varname="CF")


  		




  48
  48
  
    
mod09sd=brick("data/cloud_std.nc",varname="CFsd")


  		




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


  		




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


  		




  
  49
  
    
mod09metrics=stack(mod09a,mod09min,mod09max,mod09sd)


  		




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


  		




  51
  51
  
    

  		




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


  		




  53
  53
  
    

  		




  ......




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


  		




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


  		




  73
  73
  
    

  		




  74
  
  
    

  		




  75
  
  
    
#### get biome data


  		




  76
  
  
    
##make spatial object


  		




  77
  
  
    
cldms=cldm


  		




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


  		




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


  		




  80
  
  
    

  		




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


  		




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


  		




  83
  
  
    
    teow=teow[teow$BIOME<90,]


  		




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


  		




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


  		




  86
  
  
    
    biome=SpatialPolygonsDataFrame(biome,data=biomeid)


  		




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


  		




  88
  
  
    
}


  		




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


  		




  90
  
  
    
projection(biome)=projection(st)


  		




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


  		




  92
  
  
    

  		




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


  		




  94
  
  
    

  		




  95
  74
  
    
## get stratified sample of points from biomes for illustration


  		




  96
  
  
    
if(!file.exists("output/biomesamplepoints.csv")){


  		




  97
  
  
    
    n_biomesamples=1000


  		




  98
  
  
    
    library(multicore)


  		




  99
  
  
    
    biomesample=do.call(rbind.data.frame,mclapply(1:length(biome),function(i)


  		




  100
  
  
    
        data.frame(biome=biome$BiomeID[i],coordinates(spsample(biome[i,],n=n_biomesamples,type="stratified",nsig=2)))))


  		




  101
  
  
    
    write.csv(biomesample,"output/biomesamplepoints.csv",row.names=F)


  		




  102
  
  
    
}


  		




  103
  
  
    
biomesample=read.csv("output/biomesamplepoints.csv")


  		




  104
  
  
    
coordinates(biomesample)=c("x1","x2")


  		




  
  75
  
    
## if(!file.exists("output/biomesamplepoints.csv")){


  		




  
  76
  
    
##     n_biomesamples=1000


  		




  
  77
  
    
##     library(multicore)


  		




  
  78
  
    
##     biomesample=do.call(rbind.data.frame,mclapply(1:length(biome),function(i)


  		




  
  79
  
    
##         data.frame(biome=biome$BiomeID[i],coordinates(spsample(biome[i,],n=n_biomesamples,type="stratified",nsig=2)))))


  		




  
  80
  
    
##     write.csv(biomesample,"output/biomesamplepoints.csv",row.names=F)


  		




  
  81
  
    
## }


  		




  
  82
  
    
## biomesample=read.csv("output/biomesamplepoints.csv")


  		




  
  83
  
    
## coordinates(biomesample)=c("x1","x2")


  		




  105
  84
  
    

  		




  106
  85
  
    
#biomesample=extract(biomesample,mod09c)


  		




  107
  86
  
    

  		




  ......




  113
  92
  
    

  		




  114
  93
  
    

  		




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


  		




  116
  
  
    
png("output/Figures_%03d.png",width=5000,height=4000,res=600,pointsize=36,bg="transparent")


  		




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


  		




  117
  96
  
    

  		




  118
  
  
    
res=1e5


  		




  
  97
  
    
res=1e4


  		




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


  		




  120
  99
  
    
    


  		




  121
  100
  
    
## Figure 1: 4-panel summaries


  		




  ......




  134
  113
  
    
  margin=F,maxpixels=res,ylab="",xlab="",useRaster=T,ylim=greg$ylim)+


  		




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


  		




  136
  115
  
    

  		




  
  116
  
    
## four metics


  		




  
  117
  
    
levelplot(mod09metrics,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(space="bottom",adj=2),


  		




  
  118
  
    
  margin=F,maxpixels=res,ylab="",xlab="",useRaster=T,ylim=greg$ylim)+


  		




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


  		




  
  120
  
    

  		




  137
  121
  
    
## Monthly Means


  		




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


  		




  139
  123
  
    
  margin=F,maxpixels=res,ylab="Latitude",xlab="Longitude",useRaster=T,ylim=greg$ylim)+


  		




  ......




  142
  126
  
    
#- Monthly minimum


  		




  143
  127
  
    
#- Monthly maximum


  		




  144
  128
  
    
#- STDEV or Min-Max


  		




  145
  
  
    
#p_mac=levelplot(mod09a,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),xlab="",ylab="",useRaster=T)+


  		




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


  		




  147
  
  
    
#p_min=levelplot(mod09min,col.regions=colr(n),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)


  		




  148
  
  
    
#p_max=levelplot(mod09max,col.regions=colr(n),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)


  		




  149
  
  
    
#p_sd=levelplot(mod09sd,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)


  		




  150
  
  
    
#p3=c("Mean Cloud Frequency (%)"=p_mac,"Max Cloud Frequency (%)"=p_max,"Min Cloud Frequency (%)"=p_min,"Cloud Frequency Variability (SD)"=p_sd,x.same=T,y.same=T,merge.legends=T,layout=c(2,2))


  		




  151
  
  
    
#print(p3)


  		




  
  129
  
    
p_mac=levelplot(mod09a,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,maxpixels=res/10,colorkey=list(space="bottom",height=.75),xlab="",ylab="",useRaster=T)+


  		




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


  		




  
  131
  
    
p_min=levelplot(mod09min,col.regions=colr(n),cuts=99,margin=F,maxpixels=res/10,colorkey=list(space="bottom",height=.75),useRaster=T)+


  		




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


  		




  
  133
  
    
p_max=levelplot(mod09max,col.regions=colr(n),cuts=99,margin=F,maxpixels=res/10,colorkey=list(space="bottom",height=.75),useRaster=T)+


  		




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


  		




  
  135
  
    
p_sd=levelplot(mod09sd,col.regions=colr(n),cuts=99,at=seq(0,100,length=100),margin=F,maxpixels=res/10,colorkey=list(space="bottom",height=.75),useRaster=T)+


  		




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


  		




  
  137
  
    
p3=c("Mean Cloud Frequency (%)"=p_mac,"Max Cloud Frequency (%)"=p_max,"Min Cloud Frequency (%)"=p_min,"Cloud Frequency Variability (SD)"=p_sd,x.same=T,y.same=T,merge.legends=T,layout=c(2,2))


  		




  
  138
  
    
print(p3)


  		




  152
  139
  
    

  		




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


  		




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


  		




  ......




  206
  193
  
    
## }


  		




  207
  194
  
    

  		




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


  		




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


  		




  209
  197
  
    

  		




  210
  198
  
    
#library(BayesFactor)


  		




  211
  199
  
    
#coplot(difm~month2|lulcc,data=cldm[!is.na(cldm$dif)&!is.na(cldm$lulcc),],panel = panel.smooth)


  		




  ......




  213
  201
  
    

  		




  214
  202
  
    
dev.off()


  		




  215
  203
  
    

  		




  
  204
  
    
####################################################################


  		




  
  205
  
    
### Regional Comparisons


  		




  
  206
  
    
## Compare with worldclim and NPP


  		




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


  		




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


  		




  
  209
  
    

  		




  
  210
  
    

  		




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


  		




  
  212
  
    
regs=list(


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  
  220
  
    
  )


  		




  
  221
  
    
for(r in 1:length(regs)){


  		




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


  		




  
  223
  
    
p_map=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=1e5,colorkey=list(space="bottom",height=.75),xlab="",ylab="",main=names(regs)[r],useRaster=T)


  		




  
  224
  
    
tmac=crop(mod09a,regs[[r]])


  		




  
  225
  
    
p_mac=levelplot(tmac,col.regions=grey(seq(0,1,len=100)),cuts=100,at=seq(tmac@data@min,tmac@data@max,len=100),margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)


  		




  
  226
  
    
#p_npp=levelplot(crop(npp,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.5),zscaleLog=T,useRaster=T)  #"NPP"=p_npp,


  		




  
  227
  
    
p3=c("MOD09 Cloud Frequency (%)"=p_mac,"WorldClim Mean Annual Precip (mm)"=p_map,x.same=T,y.same=T,merge.legends=T,layout=c(2,1))


  		




  
  228
  
    
print(p3)


  		




  
  229
  
    
}


  		




  
  230
  
    
dev.off()


  		




  
  231
  
    

  		




  
  232
  
    

  		




  
  233
  
    
## reduced resolution


  		




  
  234
  
    

  		




  
  235
  
    
## read in GEWEX 1-degree data


  		




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


  		




  
  237
  
    

  		




  
  238
  
    
mod09_8km=aggregate(mod09_mac,8)


  		




  
  239
  
    

  		




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


  		




  
  241
  
    
p1=levelplot(mod09_mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  242
  
    
#p2=levelplot(mod09_8km,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  243
  
    
p3=levelplot(gewex,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  244
  
    
print(c(p1,p3,x.same=T,y.same=T,merge.legends=F))


  		




  
  245
  
    

  		




  
  246
  
    
p1=levelplot(crop(mac,regs[["Venezuela"]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  247
  
    
#p2=levelplot(crop(mod09_8km,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  248
  
    
p3=levelplot(crop(gewex,regs[["Venezuela"]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  249
  
    
print(c(MOD09=p1,GEWEX=p3,x.same=T,y.same=T,merge.legends=F))


  		




  
  250
  
    

  		




  
  251
  
    
p1=levelplot(crop(mod09_mac,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  252
  
    
#p2=levelplot(crop(mod09_8km,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  253
  
    
p3=levelplot(crop(mod09_1deg,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  
  254
  
    
print(c(p1,p3,x.same=T,y.same=T,merge.legends=F))


  		




  
  255
  
    

  		




  
  256
  
    
dev.off()


  		




  
  257
  
    

  		




  216
  258
  
    

  		




  217
  259
  
    

  		




  218
  260
  
    
## Validation table construction


  		




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


  		




  219
  262
  
    

  		




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


  		




  221
  264
  
    

  		




  ......




  230
  273
  
    
mod_bs=lapply(bs,function(bs) lm(cld_all~mod09,data=cldm[cldm$biome==bs,]))


  		




  231
  274
  
    
names(mod_bs)=biome$Biome


  		




  232
  275
  
    

  		




  233
  
  
    
screenreg(mod_bs,digits=2,single.row=F,custom.model.names=names(mod_bs))


  		




  
  276
  
    
bt=do.call(rbind,tapply(cldm$difm,list(cldm$month,cldm$biome),function(x) c(n=length(x),mean=round(mean(x,na.rm=T),1),sd=round(sd(x,na.rm=T),1))))


  		




  
  277
  
    
bt


  		




  
  278
  
    
cast(cldm,biome~month,fun=function(x) paste(round(mean(x,na.rm=T),1)," (",round(sd(x,na.rm=T),1),")",sep=""),value="difm")


  		




  
  279
  
    
melt(cast(cldm,biome~month,fun=function(x) c(mean=mean(x,na.rm=T),sd=sd(x,na.rm=T)),value="difm"))


  		




  234
  280
  
    

  		




  
  281
  
    
bt=melt(cast(cldm,biome~month2,fun=function(x) c(mean=mean(x,na.rm=T),sd=sd(x,na.rm=T)),value="mod09"))


  		




  
  282
  
    
xyplot(value~month2,data=bt)


  		




  
  283
  
    

  		




  
  284
  
    
print(xtable(bt),type="html")


  		




  235
  285
  
    

  		




  236
  286
  
    
#lm_all=lm(cld_all~mod09,data=cldm[!is.na(cldm$cld),])


  		




  237
  287
  
    
lm_mod=lm(cld~mod09+biome,data=cldm)


  		




  ......




  291
  341
  
    
cov.model <- "exponential"


  		




  292
  342
  
    

  		




  293
  343
  
    
## Run the model


  		




  294
  
  
    
n.samples <- 100


  		




  
  344
  
    
n.samples <- 500


  		




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


  		




  296
  346
  
    

  		




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


  		




  297
  348
  
    
## summarize


  		




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


  		




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


  		




  ......




  308
  359
  
    

  		




  309
  360
  
    

  		




  310
  361
  
    

  		




  
  362
  
    

  		




  311
  363
  
    
library(texreg)


  		




  312
  364
  
    
 extract.lm <- function(model) {


  		




  313
  365
  
    
     s <- summary(model)


  		




  ......




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


  		




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


  		




  362
  414
  
    

  		




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


  		




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


  		




  365
  
  
    
print(xtable(lulctl[,c("lulcc","count","mean","rmse")],digits=1),"html")


  		




  366
  
  
    

  		




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


  		




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


  		




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


  		




  
  418
  
    
    


  		




  367
  419
  
    

  		




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


  		




  369
  421
  
    
lulcrmse


  		












  

    
      climate/procedures/MOD09_Visualize.R
    
  





  7
  7
  
    
## read in global coasts for nice plotting


  		




  8
  8
  
    
library(maptools)


  		




  9
  9
  
    

  		




  10
  
  
    
data(wrld_simpl)


  		




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


  		




  12
  
  
    
coast=as(coast,"SpatialLines")


  		




  13
  10
  
    
#coast=spTransform(coast,CRS(projection(mod35)))


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  14
  16
  
    

  		




  15
  17
  
    

  		




  16
  18
  
    
#### Evaluate MOD35 Cloud data


  		




  ......




  28
  30
  
    
}


  		




  29
  31
  
    

  		




  30
  32
  
    
#### Evaluate MOD35 Cloud data


  		




  31
  
  
    
mmc=brick("~/acrobates/adamw/projects/cloud/data/mod09_clim_mean.nc",varname="CF")


  		




  
  33
  
    
mmc=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc",varname="CF")


  		




  32
  34
  
    
names(mmc)=month.name


  		




  33
  
  
    
NAvalue(mmc)=-1


  		




  34
  35
  
    

  		




  35
  36
  
    
cols=colorRampPalette(c("#000000","#00FF00","#FF0000"))#"black","blue","red"))


  		




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


  		




  36
  38
  
    
for(i in 1:12){


  		




  37
  
  
    
png(paste("output/mod09_animation_",i,".png",sep=""),width=2000,height=1000)


  		




  38
  
  
    
  print(i)


  		




  39
  
  
    
  r=mmc[[i]]


  		




  40
  
  
    
  print(levelplot(r,col.regions=cols(100),at=seq(0,100,len=100),margin=F,maxpixels=1e7,ylim=c(-60,70),


  		




  41
  
  
    
                  main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99))+


  		




  42
  
  
    
    layer(sp.lines(coast))


  		




  43
  
  
    
dev.off()


  		




  44
  
  
    
}


  		




  45
  
  
    

  		




  46
  
  
    

  		




  47
  
  
    
## climatologies


  		




  48
  
  
    
mac=brick("~/acrobates/adamw/projects/cloud/data/cloud_mean.nc",varname="CF_annual")


  		




  49
  
  
    

  		




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


  		




  51
  
  
    
levelplot(mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)+


  		




  52
  
  
    
    layer(sp.lines(coast,lwd=.5,col="black"))


  		




  53
  
  
    
dev.off()


  		




  54
  
  
    

  		




  55
  
  
    

  		




  56
  
  
    
## Compare with worldclim and NPP


  		




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


  		




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


  		




  59
  
  
    
npp=raster("/mnt/data/jetzlab/Data/environ/global/MODIS/MOD17A3/MOD17A3_Science_NPP_mean_00_12.tif",sep="")


  		




  60
  
  
    

  		




  61
  
  
    

  		




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


  		




  63
  
  
    
regs=list(


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  71
  
  
    
  )


  		




  72
  
  
    
for(r in 1:length(regs)){


  		




  73
  
  
    
p_map=levelplot(crop(wc_map,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),xlab="",ylab="",main=names(regs)[r],useRaster=T)


  		




  74
  
  
    
p_mac=levelplot(crop(mac,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)


  		




  75
  
  
    
#p_npp=levelplot(crop(npp,regs[[r]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.5),zscaleLog=T,useRaster=T)  #"NPP"=p_npp,


  		




  76
  
  
    
p3=c("WorldClim Mean Annual Precip (mm)"=p_map,"MOD09 Cloud Frequency (%)"=p_mac,x.same=T,y.same=T,merge.legends=T,layout=c(2,1))


  		




  77
  
  
    
print(p3)


  		




  
  39
  
    
    print(i)


  		




  
  40
  
    
    r=mmc[[i]]


  		




  
  41
  
    
    print(levelplot(r,col.regions=cols(100),at=seq(1,100,len=100),margin=F,maxpixels=1e6,ylim=c(-60,73),


  		




  
  42
  
    
                    main=paste(month.name[i]),cex.main=3,scales=list(draw=F),cuts=99,ylab="",xlab="")+


  		




  
  43
  
    
                        layer(panel.polygon(x=c(-180,-180,180,180),y=c(-90,90,90,-90),col="black"),under=T)+


  		




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


  		




  78
  45
  
    
}


  		




  79
  46
  
    
dev.off()


  		




  80
  47
  
    

  		




  81
  48
  
    

  		




  82
  
  
    
## reduced resolution


  		




  83
  
  
    

  		




  84
  
  
    
## read in GEWEX 1-degree data


  		




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


  		




  86
  
  
    

  		




  87
  
  
    
mod09_8km=aggregate(mod09_mac,8)


  		




  88
  
  
    

  		




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


  		




  90
  
  
    
p1=levelplot(mod09_mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  91
  
  
    
#p2=levelplot(mod09_8km,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  92
  
  
    
p3=levelplot(gewex,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  93
  
  
    
print(c(p1,p3,x.same=T,y.same=T,merge.legends=F))


  		




  94
  
  
    

  		




  95
  
  
    
p1=levelplot(crop(mac,regs[["Venezuela"]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  96
  
  
    
#p2=levelplot(crop(mod09_8km,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  97
  
  
    
p3=levelplot(crop(gewex,regs[["Venezuela"]]),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  98
  
  
    
print(c(MOD09=p1,GEWEX=p3,x.same=T,y.same=T,merge.legends=F))


  		




  99
  
  
    

  		




  100
  
  
    
p1=levelplot(crop(mod09_mac,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  101
  
  
    
#p2=levelplot(crop(mod09_8km,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  102
  
  
    
p3=levelplot(crop(mod09_1deg,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)


  		




  103
  
  
    
print(c(p1,p3,x.same=T,y.same=T,merge.legends=F))


  		




  104
  
  
    

  		




  105
  
  
    
dev.off()


  		




  106
  
  
    

  		












  

    
      climate/procedures/NDP-026D.R
    
  





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


  		




  
  2
  
    
### to validate MODIS cloud frequencies


  		




  2
  3
  
    

  		




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


  		




  4
  5
  
    

  		




  ......




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


  		




  126
  127
  
    

  		




  127
  128
  
    

  		




  
  129
  
    
### Add biome data


  		




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


  		




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


  		




  
  132
  
    
    teow=teow[teow$BIOME<90,]


  		




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


  		




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


  		




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


  		




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


  		




  
  137
  
    
}


  		




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


  		




  
  139
  
    
projection(biome)=projection(st)


  		




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


  		




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


  		




  
  142
  
    
st$biome=biome$BiomeID[dists]


  		




  
  143
  
    

  		




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


  		




  
  145
  
    

  		




  
  146
  
    

  		




  128
  147
  
    
## write out the tables


  		




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


  		




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


  		




  131
  
  
    

  		




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


  		




  132
  151
  
    
#########################################################################


  		




  133
  152
  
    

  		












  

    
      climate/procedures/ee_compile.R
    
  





  139
  139
  
    
# Regressions through time by season


  		




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


  		




  141
  141
  
    

  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  146
  146
  
    

  		




  147
  147
  
    

  		




  148
  148
  
    

  		




  
  149
  
    
## Daily animations


  		




  
  150
  
    
regs=list(


  		




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


  		




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


  		




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


  		




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


  		




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


  		




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


  		




  
  157
  
    
    )


  		




  
  158
  
    

  		




  
  159
  
    
r=1


  		




  
  160
  
    

  		




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


  		




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


  		




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


  		




  
  164
  
    

  		




  
  165
  
    

  		




  149
  166
  
    

  		




  150
  167
  
    

  		




  151
  168
  
    
### Long term summaries


  		












Also available in:  Unified diff