Revision 1dc46eb9
Added by Adam Wilson almost 11 years ago
| 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 |
|
Also available in: Unified diff
Updating figures