Revision f07e0e66
Added by Adam Wilson almost 11 years ago
| climate/procedures/NDP-026D.R | ||
|---|---|---|
| 49 | 49 |
## drop missing values |
| 50 | 50 |
cld=cld[,!grepl("Fq|AWP|NC",colnames(cld))]
|
| 51 | 51 |
cld$Amt[cld$Amt<0]=NA |
| 52 |
#cld$Fq[cld$Fq<0]=NA
|
|
| 53 |
#cld$AWP[cld$AWP<0]=NA |
|
| 54 |
#cld$NC[cld$NC<0]=NA
|
|
| 55 |
#cld=cld[cld$Nobs>0,]
|
|
| 52 |
cld$Amt=cld$Amt/100
|
|
| 53 |
|
|
| 54 |
## calculate means and sds for full record (1970-2009)
|
|
| 55 |
Nobsthresh=20 #minimum number of observations to include
|
|
| 56 | 56 |
|
| 57 |
## calculate means and sds |
|
| 58 | 57 |
cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){
|
| 59 | 58 |
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))})) |
|
| 59 |
month=x$month[1], |
|
| 60 |
StaID=x$StaID[1], |
|
| 61 |
cld_all=mean(x$Amt[x$Nobs>=Nobsthresh],na.rm=T), # full record |
|
| 62 |
cldsd_all=sd(x$Amt[x$Nobs>=Nobsthresh],na.rm=T), |
|
| 63 |
cld=mean(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T), #only MODIS epoch |
|
| 64 |
cldsd=sd(x$Amt[x$YR>=2000&x$Nobs>=Nobsthresh],na.rm=T))})) |
|
| 64 | 65 |
cldm[,c("lat","lon")]=coordinates(st)[match(cldm$StaID,st$id),c("lat","lon")]
|
| 65 | 66 |
|
| 66 | 67 |
|
| 68 |
|
|
| 67 | 69 |
## add the MOD09 data to cld |
| 68 | 70 |
#### Evaluate MOD35 Cloud data |
| 69 | 71 |
mod09=brick("~/acrobates/adamw/projects/cloud/data/cloud_ymonmean.nc")
|
| ... | ... | |
| 91 | 93 |
})#,mc.cores=3) |
| 92 | 94 |
|
| 93 | 95 |
## read it back in |
| 94 |
mod09st=read.csv("valid.csv",header=F)[,-c(1,2)]
|
|
| 95 |
|
|
| 96 |
colnames(mod09st)=c(names(mod09)[-1],"id") |
|
| 97 |
mod09stl=melt(mod09st,id.vars=c("id","sd"))
|
|
| 96 |
mod09st=read.csv("valid.csv",header=F)[,-c(1)]
|
|
| 97 |
colnames(mod09st)=c(names(mod09),"id","type") |
|
| 98 |
mod09stl=melt(mod09st,id.vars=c("id","type"))
|
|
| 98 | 99 |
mod09stl[,c("year","month")]=do.call(rbind,strsplit(sub("X","",mod09stl$variable),"[.]"))[,1:2]
|
| 99 | 100 |
mod09stl$value[mod09stl$value<0]=NA |
| 101 |
mod09stl=cast(mod09stl,id+year+month~type,value="value") |
|
| 100 | 102 |
|
| 101 | 103 |
## add it to cld |
| 102 |
cldm$mod09=mod09stl$value[match(paste(cldm$StaID,cldm$month),paste(mod09stl$id,as.numeric(mod09stl$month)))] |
|
| 104 |
cldm$mod09=mod09stl$mean[match(paste(cldm$StaID,cldm$month),paste(mod09stl$id,as.numeric(mod09stl$month)))] |
|
| 105 |
cldm$mod09sd=mod09stl$sd[match(paste(cldm$StaID,cldm$month),paste(mod09stl$id,as.numeric(mod09stl$month)))] |
|
| 103 | 106 |
|
| 104 | 107 |
|
| 105 | 108 |
## LULC |
| ... | ... | |
| 121 | 124 |
cldm$lulc=lulcst$lulc[match(cldm$StaID,lulcst$id)] |
| 122 | 125 |
cldm$lulcc=IGBP$class[match(cldm$lulc,IGBP$ID)] |
| 123 | 126 |
|
| 124 |
## update cld column names |
|
| 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")]
|
|
| 128 |
|
|
| 129 |
## calculate means and sds |
|
| 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")]
|
|
| 140 |
|
|
| 141 |
## means by year |
|
| 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")]
|
|
| 152 |
|
|
| 153 |
## overall mean |
|
| 154 |
clda=do.call(rbind.data.frame,by(cldm,list(StaID=as.factor(cldm$StaID)),function(x){
|
|
| 155 |
data.frame( |
|
| 156 |
StaID=x$StaID[1], |
|
| 157 |
lulc=x$lulc[1], |
|
| 158 |
mod09=mean(x$mod09,na.rm=T), |
|
| 159 |
mod09sd=sd(x$mod09,na.rm=T), |
|
| 160 |
cld=mean(x$cld,na.rm=T), |
|
| 161 |
cldsd=sd(x$cld,na.rm=T))})) |
|
| 162 |
clda[,c("lat","lon")]=coordinates(st)[match(clda$StaID,st$id),c("lat","lon")]
|
|
| 163 |
|
|
| 164 | 127 |
|
| 165 | 128 |
## write out the tables |
| 166 | 129 |
write.csv(cld,file="cld.csv",row.names=F) |
| 167 |
#write.csv(cldy,file="cldy.csv",row.names=F) |
|
| 168 | 130 |
write.csv(cldm,file="cldm.csv",row.names=F) |
| 169 |
write.csv(clda,file="clda.csv",row.names=F) |
|
| 170 | 131 |
|
| 171 | 132 |
######################################################################### |
| 172 | 133 |
|
Also available in: Unified diff
udpates to validation and addition of initial biodiversity script