Project

General

Profile

Download (1.6 KB) Statistics
| Branch: | Revision:
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)