|
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 -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)
|
