| 1 |
4ef959c2
|
Adam M. Wilson
|
#! /bin/R
|
| 2 |
|
|
### Script to download and process the NDP-026D station cloud dataset
|
| 3 |
|
|
setwd("~/acrobates/projects/interp/data/NDP026D")
|
| 4 |
|
|
|
| 5 |
|
|
## available here http://cdiac.ornl.gov/epubs/ndp/ndp026d/ndp026d.html
|
| 6 |
|
|
|
| 7 |
|
|
|
| 8 |
|
|
## Get station locations
|
| 9 |
|
|
system("wget -nd http://cdiac.ornl.gov/ftp/ndp026d/cat01/01_STID -P data/")
|
| 10 |
|
|
st=read.table("data/01_STID",skip=1)
|
| 11 |
|
|
colnames(st)=c("StaID","LAT","LON","ELEV","ny1","fy1","ly1","ny7","fy7","ly7","SDC","b5c")
|
| 12 |
|
|
st$lat=st$LAT/100
|
| 13 |
|
|
st$lon=st$LON/100
|
| 14 |
|
|
st$lon[st$lon>180]=st$lon[st$lon>180]-360
|
| 15 |
|
|
|
| 16 |
|
|
## check a plot
|
| 17 |
|
|
plot(lat~lon,data=st,pch=16,cex=.5)
|
| 18 |
|
|
|
| 19 |
|
|
|
| 20 |
|
|
## get monthly mean cloud amount MMCA
|
| 21 |
|
|
system("wget -N -nd ftp://cdiac.ornl.gov/pub/ndp026d/cat67_78/* -A '.tc.Z' -P data/")
|
| 22 |
|
|
system("gunzip data/*.Z")
|
| 23 |
|
|
|
| 24 |
|
|
#f121=c(6,6,6,7,6,7,6,2) #format 121
|
| 25 |
|
|
#c121=c("StaID","NobD","AvgDy","NobN","AvgNt","NobDN","AvgDN","Acode")
|
| 26 |
|
|
f162=c(5,5,4,7,7,7,4) #format 121
|
| 27 |
|
|
c162=c("StaID","YR","Nobs","Amt","Fq","AWP","NC")
|
| 28 |
|
|
|
| 29 |
|
|
cld=do.call(rbind.data.frame,lapply(sprintf("%02d",1:12),function(m) {
|
| 30 |
|
|
d=read.fwf(list.files("data",pattern=paste("MNYDC.",m,".tc",sep=""),full=T),skip=1,widths=f162)
|
| 31 |
|
|
colnames(d)=c162
|
| 32 |
|
|
d$month=as.numeric(m)
|
| 33 |
|
|
return(d)}
|
| 34 |
|
|
))
|
| 35 |
|
|
|
| 36 |
|
|
cld[,c("lat","lon")]=st[match(st$StaID,cld$StaID),c("lat","lon")]
|
| 37 |
|
|
|
| 38 |
|
|
## calculate means
|
| 39 |
|
|
cldm=by(cld,list(as.factor(cld$month),as.factor(cld$StaID)),function(x){
|
| 40 |
|
|
data.frame(Amt=mean(x$Amt[x$Nobs>20],na.rm=T))})
|
| 41 |
|
|
|
| 42 |
|
|
|
| 43 |
|
|
## add color val
|
| 44 |
|
|
cld$Amt[cld$Amt<0]=NA
|
| 45 |
|
|
cld$Fq[cld$Fq<0]=NA
|
| 46 |
|
|
cld$AWP[cld$AWP<0]=NA
|
| 47 |
|
|
cld$NC[cld$NC<0]=NA
|
| 48 |
|
|
|
| 49 |
|
|
cld$col=cut(cld$Amt/100,quantile(cld$Amt/100,seq(0,1,len=5),na.rm=T))
|
| 50 |
|
|
|
| 51 |
|
|
library(lattice)
|
| 52 |
|
|
xyplot(lat~lon|as.factor(YR)+as.factor(month),groups=col,data=cld,pch=16,cex=.2,auto.key=T)
|
