|
1
|
#! /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 -N -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
|
## drop missing values
|
|
39
|
cld$Amt[cld$Amt<0]=NA
|
|
40
|
cld$Fq[cld$Fq<0]=NA
|
|
41
|
cld$AWP[cld$AWP<0]=NA
|
|
42
|
cld$NC[cld$NC<0]=NA
|
|
43
|
|
|
44
|
## calculate means
|
|
45
|
cldm=do.call(rbind.data.frame,by(cld,list(month=as.factor(cld$month),StaID=as.factor(cld$StaID)),function(x){
|
|
46
|
data.frame(month=x$month[1],StaID=x$StaID[1],Amt=mean(x$Amt[x$Nobs>20],na.rm=T))}))
|
|
47
|
cldm[,c("lat","lon")]=st[match(st$StaID,cldm$StaID),c("lat","lon")]
|
|
48
|
|
|
49
|
|
|
50
|
|
|
51
|
## write out the table
|
|
52
|
write.csv(cldm,file="cldm.csv")
|
|
53
|
|
|
54
|
|
|
55
|
##################
|
|
56
|
###
|
|
57
|
cldm=read.csv("cldm.csv")
|
|
58
|
|
|
59
|
## add a color key
|
|
60
|
cldm$col=cut(cldm$Amt/100,quantile(cldm$Amt/100,seq(0,1,len=5),na.rm=T))
|
|
61
|
|
|
62
|
library(lattice)
|
|
63
|
xyplot(lat~lon|+month,groups=col,data=cldm,pch=16,cex=.2,auto.key=T)
|
|
64
|
|
