Revision 1dc46eb9
Added by Adam Wilson almost 11 years ago
| climate/procedures/MOD09_CloudFigures.R | ||
|---|---|---|
| 41 | 41 |
projection(cldys)=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
|
| 42 | 42 |
|
| 43 | 43 |
#### Evaluate MOD35 Cloud data |
| 44 |
mod09=brick("data/cloud_daily.nc")
|
|
| 44 |
mod09=brick("data/cloud_monthly.nc")
|
|
| 45 | 45 |
mod09c=brick("data/cloud_ymonmean.nc",varname="CF");names(mod09c)=month.name
|
| 46 | 46 |
mod09a=brick("data/cloud_mean.nc",varname="CF_annual")#;names(mod09c)=month.name
|
| 47 | 47 |
|
| ... | ... | |
| 76 | 76 |
|
| 77 | 77 |
## Figure 1: 4-panel summaries |
| 78 | 78 |
#- Annual average |
| 79 |
levelplot(mod09a,col.regions=colr(100),cuts=100,at=seq(0,100,len=100),colorkey=list(space="bottom",adj=1),
|
|
| 79 |
levelplot(mod09a,col.regions=colr(n),cuts=100,at=seq(0,100,len=100),colorkey=list(space="bottom",adj=1),
|
|
| 80 | 80 |
margin=F,maxpixels=1e6,ylab="Latitude",xlab="Longitude",useRaster=T)+ |
| 81 | 81 |
layer(sp.lines(coast,col="black"),under=F) |
| 82 | 82 |
#- Monthly minimum |
| 83 | 83 |
#- Monthly maximum |
| 84 | 84 |
#- STDEV or Min-Max |
| 85 |
p_mac=levelplot(mac,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)
|
|
| 86 |
p_min=levelplot(mod09min,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)
|
|
| 87 |
p_max=levelplot(mod09max,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)
|
|
| 88 |
p_sd=levelplot(mod09sd,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)
|
|
| 85 |
p_mac=levelplot(mod09a,col.regions=colr(n),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),xlab="",ylab="",main=names(regs)[r],useRaster=T)
|
|
| 86 |
p_min=levelplot(mod09min,col.regions=colr(n),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)
|
|
| 87 |
p_max=levelplot(mod09max,col.regions=colr(n),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)
|
|
| 88 |
p_sd=levelplot(mod09sd,col.regions=colr(n),cuts=99,margin=F,maxpixels=1e5,colorkey=list(space="bottom",height=.75),useRaster=T)
|
|
| 89 | 89 |
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))
|
| 90 | 90 |
print(p3) |
| 91 | 91 |
|
| climate/procedures/MOD09_Visualize.R | ||
|---|---|---|
| 45 | 45 |
|
| 46 | 46 |
|
| 47 | 47 |
## climatologies |
| 48 |
mac=brick("~/acrobates/adamw/projects/cloud/data/mod09_clim_mac.nc",varname="CF_annual")
|
|
| 48 |
mac=brick("~/acrobates/adamw/projects/cloud/data/cloud_mean.nc",varname="CF_annual")
|
|
| 49 | 49 |
|
| 50 | 50 |
pdf("output/mod09_climatology.pdf",width=11,height=8.5)
|
| 51 | 51 |
levelplot(mac,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e6)+ |
| ... | ... | |
| 82 | 82 |
## reduced resolution |
| 83 | 83 |
|
| 84 | 84 |
## read in GEWEX 1-degree data |
| 85 |
gewex=raster("data/gewex/CA_PATMOSX_NOAA.nc")
|
|
| 85 |
gewex=mean(brick("data/gewex/CA_PATMOSX_NOAA.nc",varname="a_CA"))
|
|
| 86 | 86 |
|
| 87 | 87 |
mod09_8km=aggregate(mod09_mac,8) |
| 88 | 88 |
|
| 89 | 89 |
pdf("output/mod09_resolution.pdf",width=11,height=8.5)
|
| 90 | 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(mod09_1deg,col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5)
|
|
| 93 |
print(c(p1,p2,p3,x.same=T,y.same=T,merge.legends=F))
|
|
| 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 | 94 |
|
| 95 |
|
|
| 96 |
p1=levelplot(crop(mod09_mac,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5) |
|
| 97 |
p2=levelplot(crop(mod09_8km,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5) |
|
| 98 |
p3=levelplot(crop(mod09_1deg,reg2),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5) |
|
| 99 |
print(c(p1,p2,p3,x.same=T,y.same=T,merge.legends=F)) |
|
| 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)) |
|
| 100 | 99 |
|
| 101 | 100 |
p1=levelplot(crop(mod09_mac,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5) |
| 102 |
p2=levelplot(crop(mod09_8km,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)
|
|
| 103 | 102 |
p3=levelplot(crop(mod09_1deg,reg3),col.regions=grey(seq(0,1,len=100)),cuts=99,margin=F,max.pixels=1e5) |
| 104 |
print(c(p1,p2,p3,x.same=T,y.same=T,merge.legends=F))
|
|
| 103 |
print(c(p1,p3,x.same=T,y.same=T,merge.legends=F)) |
|
| 105 | 104 |
|
| 106 | 105 |
dev.off() |
| 107 | 106 |
|
| climate/procedures/NDP-026D.R | ||
|---|---|---|
| 23 | 23 |
write.csv(st,"stations.csv",row.names=F) |
| 24 | 24 |
coordinates(st)=c("lon","lat")
|
| 25 | 25 |
projection(st)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" |
| 26 |
st@data[,c("lon","lat")]=coordinates(st)
|
|
| 26 | 27 |
|
| 27 | 28 |
## download data |
| 28 | 29 |
system("wget -N -nd ftp://cdiac.ornl.gov/pub/ndp026d/cat67_78/* -A '.tc.Z' -P data/")
|
| ... | ... | |
| 53 | 54 |
#cld$NC[cld$NC<0]=NA |
| 54 | 55 |
#cld=cld[cld$Nobs>0,] |
| 55 | 56 |
|
| 57 |
## calculate means and sds |
|
| 58 |
cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){
|
|
| 59 |
data.frame( |
|
| 60 |
month=x$month[1], |
|
| 61 |
StaID=x$StaID[1], |
|
| 62 |
cld=mean(x$cld[x$Nobs>60],na.rm=T), |
|
| 63 |
cldsd=sd(x$cld[x$Nobs>60],na.rm=T))})) |
|
| 64 |
cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
|
| 65 |
|
|
| 66 |
|
|
| 56 | 67 |
## add the MOD09 data to cld |
| 57 | 68 |
#### Evaluate MOD35 Cloud data |
| 58 | 69 |
mod09=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc")
|
| 70 |
mod09std=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonstd.nc")
|
|
| 59 | 71 |
|
| 60 | 72 |
## overlay the data with 32km diameter (16km radius) buffer |
| 61 | 73 |
## buffer size from Dybbroe, et al. (2005) doi:10.1175/JAM-2189.1. |
| 62 | 74 |
buf=16000 |
| 63 |
bins=cut(1:nrow(st),100) |
|
| 64 |
if(file.exists("valid.csv")) file.remove("valid.csv")
|
|
| 75 |
bins=cut(st$lat,10) |
|
| 76 |
rerun=F |
|
| 77 |
if(rerun&file.exists("valid.csv")) file.remove("valid.csv")
|
|
| 65 | 78 |
mod09sta=lapply(levels(bins),function(lb) {
|
| 66 | 79 |
l=which(bins==lb) |
| 80 |
## mean |
|
| 67 | 81 |
td=extract(mod09,st[l,],buffer=buf,fun=mean,na.rm=T,df=T) |
| 68 | 82 |
td$id=st$id[l] |
| 83 |
td$type="mean" |
|
| 84 |
## std |
|
| 85 |
td2=extract(mod09std,st[l,],buffer=buf,fun=mean,na.rm=T,df=T) |
|
| 86 |
td2$id=st$id[l] |
|
| 87 |
td2$type="sd" |
|
| 69 | 88 |
print(lb)#as.vector(c(l,td[,1:4]))) |
| 70 |
write.table(td,"valid.csv",append=T,col.names=F,quote=F,sep=",",row.names=F)
|
|
| 89 |
write.table(rbind(td,td2),"valid.csv",append=T,col.names=F,quote=F,sep=",",row.names=F)
|
|
| 71 | 90 |
td |
| 72 | 91 |
})#,mc.cores=3) |
| 73 | 92 |
|
| ... | ... | |
| 75 | 94 |
mod09st=read.csv("valid.csv",header=F)[,-c(1,2)]
|
| 76 | 95 |
|
| 77 | 96 |
colnames(mod09st)=c(names(mod09)[-1],"id") |
| 78 |
mod09stl=melt(mod09st,id.vars="id")
|
|
| 97 |
mod09stl=melt(mod09st,id.vars=c("id","sd"))
|
|
| 79 | 98 |
mod09stl[,c("year","month")]=do.call(rbind,strsplit(sub("X","",mod09stl$variable),"[.]"))[,1:2]
|
| 99 |
mod09stl$value[mod09stl$value<0]=NA |
|
| 80 | 100 |
|
| 81 | 101 |
## add it to cld |
| 82 |
cld$mod09=mod09stl$value[match(paste(cld$StaID,cld$YR,cld$month),paste(mod09stl$id,mod09stl$year,as.numeric(mod09stl$month)))]
|
|
| 102 |
cldm$mod09=mod09stl$value[match(paste(cldm$StaID,cldm$month),paste(mod09stl$id,as.numeric(mod09stl$month)))]
|
|
| 83 | 103 |
|
| 84 | 104 |
|
| 85 | 105 |
## LULC |
| ... | ... | |
| 98 | 118 |
lulcst=extract(lulc,st,fun=Mode,buffer=buf,df=T) |
| 99 | 119 |
colnames(lulcst)=c("id","lulc")
|
| 100 | 120 |
## add it to cld |
| 101 |
cld$lulc=lulcst$lulc[match(cld$StaID,lulcst$id)]
|
|
| 102 |
cld$lulcc=IGBP$class[match(cld$lulc,IGBP$ID)]
|
|
| 121 |
cldm$lulc=lulcst$lulc[match(cldm$StaID,lulcst$id)]
|
|
| 122 |
cldm$lulcc=IGBP$class[match(cldm$lulc,IGBP$ID)]
|
|
| 103 | 123 |
|
| 104 | 124 |
## update cld column names |
| 105 |
colnames(cld)[grep("Amt",colnames(cld))]="cld"
|
|
| 106 |
cld$cld=cld$cld/100
|
|
| 107 |
cld[,c("lat","lon")]=coordinates(st)[match(cld$StaID,st$id),c("lat","lon")]
|
|
| 125 |
colnames(cldm)[grep("Amt",colnames(cldm))]="cld"
|
|
| 126 |
cldm$cld=cldm$cld/100
|
|
| 127 |
cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
|
| 108 | 128 |
|
| 109 | 129 |
## calculate means and sds |
| 110 |
cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){
|
|
| 111 |
data.frame( |
|
| 112 |
month=x$month[1], |
|
| 113 |
lulc=x$lulc[1], |
|
| 114 |
StaID=x$StaID[1], |
|
| 115 |
mod09=mean(x$mod09,na.rm=T), |
|
| 116 |
mod09sd=sd(x$mod09,na.rm=T), |
|
| 117 |
cld=mean(x$cld[x$Nobs>50],na.rm=T), |
|
| 118 |
cldsd=sd(x$cld[x$Nobs>50],na.rm=T))})) |
|
| 119 |
cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
|
| 130 |
#cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){
|
|
| 131 |
# data.frame(
|
|
| 132 |
# month=x$month[1],
|
|
| 133 |
# lulc=x$lulc[1],
|
|
| 134 |
# StaID=x$StaID[1],
|
|
| 135 |
# mod09=mean(x$mod09,na.rm=T),
|
|
| 136 |
# mod09sd=sd(x$mod09,na.rm=T),
|
|
| 137 |
# cld=mean(x$cld[x$Nobs>50],na.rm=T),
|
|
| 138 |
# cldsd=sd(x$cld[x$Nobs>50],na.rm=T))}))
|
|
| 139 |
#cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
|
| 120 | 140 |
|
| 121 | 141 |
## means by year |
| 122 |
cldy=do.call(rbind.data.frame,by(cld,list(year=as.factor(cld$YR),StaID=as.factor(cld$StaID)),function(x){
|
|
| 123 |
data.frame( |
|
| 124 |
year=x$YR[1], |
|
| 125 |
StaID=x$StaID[1], |
|
| 126 |
lulc=x$lulc[1], |
|
| 127 |
mod09=mean(x$mod09,na.rm=T), |
|
| 128 |
mod09sd=sd(x$mod09,na.rm=T), |
|
| 129 |
cld=mean(x$cld[x$Nobs>50]/100,na.rm=T), |
|
| 130 |
cldsd=sd(x$cld[x$Nobs>50]/100,na.rm=T))})) |
|
| 131 |
cldy[,c("lat","lon")]=coordinates(st)[match(cldy$StaID,st$id),c("lat","lon")]
|
|
| 142 |
#cldy=do.call(rbind.data.frame,by(cld,list(year=as.factor(cld$YR),StaID=as.factor(cld$StaID)),function(x){
|
|
| 143 |
# data.frame(
|
|
| 144 |
# year=x$YR[1],
|
|
| 145 |
# StaID=x$StaID[1],
|
|
| 146 |
# lulc=x$lulc[1],
|
|
| 147 |
# mod09=mean(x$mod09,na.rm=T),
|
|
| 148 |
# mod09sd=sd(x$mod09,na.rm=T),
|
|
| 149 |
# cld=mean(x$cld[x$Nobs>50]/100,na.rm=T),
|
|
| 150 |
# cldsd=sd(x$cld[x$Nobs>50]/100,na.rm=T))}))
|
|
| 151 |
#cldy[,c("lat","lon")]=coordinates(st)[match(cldy$StaID,st$id),c("lat","lon")]
|
|
| 132 | 152 |
|
| 133 | 153 |
## overall mean |
| 134 |
clda=do.call(rbind.data.frame,by(cld,list(StaID=as.factor(cld$StaID)),function(x){
|
|
| 154 |
clda=do.call(rbind.data.frame,by(cldm,list(StaID=as.factor(cldm$StaID)),function(x){
|
|
| 135 | 155 |
data.frame( |
| 136 | 156 |
StaID=x$StaID[1], |
| 137 | 157 |
lulc=x$lulc[1], |
| 138 | 158 |
mod09=mean(x$mod09,na.rm=T), |
| 139 | 159 |
mod09sd=sd(x$mod09,na.rm=T), |
| 140 |
cld=mean(x$cld[x$Nobs>10],na.rm=T),
|
|
| 141 |
cldsd=sd(x$cld[x$Nobs>10],na.rm=T))}))
|
|
| 160 |
cld=mean(x$cld,na.rm=T), |
|
| 161 |
cldsd=sd(x$cld,na.rm=T))})) |
|
| 142 | 162 |
clda[,c("lat","lon")]=coordinates(st)[match(clda$StaID,st$id),c("lat","lon")]
|
| 143 | 163 |
|
| 144 | 164 |
|
| 145 | 165 |
## write out the tables |
| 146 | 166 |
write.csv(cld,file="cld.csv",row.names=F) |
| 147 |
write.csv(cldy,file="cldy.csv",row.names=F) |
|
| 167 |
#write.csv(cldy,file="cldy.csv",row.names=F)
|
|
| 148 | 168 |
write.csv(cldm,file="cldm.csv",row.names=F) |
| 149 | 169 |
write.csv(clda,file="clda.csv",row.names=F) |
| 150 | 170 |
|
| climate/procedures/ee_compile.R | ||
|---|---|---|
| 114 | 114 |
### generate the monthly mean and sd |
| 115 | 115 |
#system(paste("cdo -P 10 -O merge -ymonmean data/mod09.nc -chname,CF,CF_sd -ymonstd data/mod09.nc data/mod09_clim.nc"))
|
| 116 | 116 |
system(paste("cdo -f nc4c -O -ymonmean data/cloud_monthly.nc data/cloud_ymonmean.nc"))
|
| 117 |
system(paste("cdo -f nc4c -O -chname,CF,CF_sd -ymonstd data/cloud_monthly.nc data/cloud_ymonsd.nc"))
|
|
| 117 |
## standard deviations, had to break to limit memory usage |
|
| 118 |
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"))
|
|
| 119 |
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"))
|
|
| 120 |
system(paste("cdo -f nc4c -O mergetime data/cloud_ymonsd_1-6.nc data/cloud_ymonsd_7-12.nc data/cloud_ymonstd.nc"))
|
|
| 121 |
|
|
| 118 | 122 |
|
| 119 | 123 |
#if(!file.exists("data/mod09_metrics.nc")) {
|
| 120 | 124 |
system("cdo -f nc4c -chname,CF,CFmin -timmin data/cloud_ymonmean.nc data/cloud_min.nc")
|
Also available in: Unified diff
Updating figures