Revision d86b0a4a
Added by Adam Wilson almost 11 years ago
| climate/procedures/NDP-026D.R | ||
|---|---|---|
| 25 | 25 |
st=st[,c("StaID","ELEV","lat","lon")]
|
| 26 | 26 |
colnames(st)=c("id","elev","lat","lon")
|
| 27 | 27 |
write.csv(st,"stations.csv",row.names=F) |
| 28 |
|
|
| 28 |
coordinates(st)=c("lon","lat")
|
|
| 29 | 29 |
## download data |
| 30 | 30 |
system("wget -N -nd ftp://cdiac.ornl.gov/pub/ndp026d/cat67_78/* -A '.tc.Z' -P data/")
|
| 31 | 31 |
|
| ... | ... | |
| 45 | 45 |
)) |
| 46 | 46 |
|
| 47 | 47 |
## add lat/lon |
| 48 |
cld[,c("lat","lon")]=st[match(cld$StaID,st$StaID),c("lat","lon")]
|
|
| 48 |
cld[,c("lat","lon")]=st[match(cld$StaID,st$id),c("lat","lon")]
|
|
| 49 | 49 |
|
| 50 | 50 |
## drop missing values |
| 51 |
cld=cld[,!grepl("Fq|AWP|NC",colnames(cld))]
|
|
| 51 | 52 |
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,] |
|
| 53 |
#cld$Fq[cld$Fq<0]=NA
|
|
| 54 |
#cld$AWP[cld$AWP<0]=NA
|
|
| 55 |
#cld$NC[cld$NC<0]=NA
|
|
| 56 |
#cld=cld[cld$Nobs>0,]
|
|
| 56 | 57 |
|
| 57 | 58 |
## add the MOD09 data to cld |
| 58 | 59 |
#### Evaluate MOD35 Cloud data |
| ... | ... | |
| 61 | 62 |
## overlay the data with 32km diameter (16km radius) buffer |
| 62 | 63 |
## buffer size from Dybbroe, et al. (2005) doi:10.1175/JAM-2189.1. |
| 63 | 64 |
buf=16000 |
| 64 |
mod09st=do.call(cbind.data.frame,mclapply(1:nlayers(mod09),function(l) extract(mod09[[l]],st,buffer=buf,fun=mean,na.rm=T,df=T)[,2])) |
|
| 65 |
mod09sta=lapply(1:nlayers(mod09),function(l) {print(l); extract(mod09[[l]],st,buffer=buf,fun=mean,na.rm=T,df=T)[,2]})
|
|
| 66 |
mod09st=do.call(cbind.data.frame,mod09sta) |
|
| 65 | 67 |
#mod09st=mod09st[,!is.na(colnames(mod09st))] |
| 66 | 68 |
colnames(mod09st)=names(mod09) |
| 67 | 69 |
mod09st$id=st$id |
| ... | ... | |
| 82 | 84 |
levels(lulc)=list(IGBP) |
| 83 | 85 |
#lulc=crop(lulc,mod09) |
| 84 | 86 |
Mode <- function(x) {
|
| 85 |
ux <- unique(x)
|
|
| 86 |
ux[which.max(tabulate(match(x, ux)))] |
|
| 87 |
ux <- na.omit(unique(x))
|
|
| 88 |
ux[which.max(tabulate(match(x, ux)),na.rm=T)]
|
|
| 87 | 89 |
} |
| 88 |
lulcst=extract(lulc,st,fun=Mode,buffer=buf,df=T) |
|
| 90 |
lulcst=extract(lulc,st,fun=Mode,na.rm=T,buffer=buf,df=T)
|
|
| 89 | 91 |
colnames(lulcst)=c("id","lulc")
|
| 90 | 92 |
## add it to cld |
| 91 | 93 |
cld$lulc=lulcst$lulc[match(cld$StaID,lulcst$id)] |
| 92 |
cld$lulc=factor(cld$lulc,labels=IGBP$class)
|
|
| 94 |
cld$lulc=factor(as.integer(cld$lulc),labels=IGBP$class)
|
|
| 93 | 95 |
|
| 94 | 96 |
## update cld column names |
| 95 | 97 |
colnames(cld)[grep("Amt",colnames(cld))]="cld"
|
| ... | ... | |
| 135 | 137 |
write.csv(cld,file="cld.csv",row.names=F) |
| 136 | 138 |
write.csv(cldy,file="cldy.csv",row.names=F) |
| 137 | 139 |
write.csv(cldm,file="cldm.csv",row.names=F) |
| 138 |
write.csv(clda,file="clda.csv",row.names=F)
|
|
| 139 |
|
|
| 140 |
write.csv(clda,file="clda.csv",row.names=F |
|
| 141 |
) |
|
| 140 | 142 |
######################################################################### |
| 141 | 143 |
################## |
| 142 | 144 |
### |
Also available in: Unified diff
updates to MOD09 visualizations