Project

General

Profile

« Previous | Next » 

Revision a6d44f2d

Added by Adam Wilson over 11 years ago

Added script to process NDP-026D station cloud climatologies

View differences:

climate/procedures/NDP-026D.R
6 6 

  		




  7
  7
  
    

  		




  8
  8
  
    
## Get station locations


  		




  9
  
  
    
system("wget -nd http://cdiac.ornl.gov/ftp/ndp026d/cat01/01_STID -P data/")


  		




  
  9
  
    
system("wget -N -nd http://cdiac.ornl.gov/ftp/ndp026d/cat01/01_STID -P data/")


  		




  10
  10
  
    
st=read.table("data/01_STID",skip=1)


  		




  11
  11
  
    
colnames(st)=c("StaID","LAT","LON","ELEV","ny1","fy1","ly1","ny7","fy7","ly7","SDC","b5c")


  		




  12
  12
  
    
st$lat=st$LAT/100


  		




  ......




  35
  35
  
    

  		




  36
  36
  
    
cld[,c("lat","lon")]=st[match(st$StaID,cld$StaID),c("lat","lon")]


  		




  37
  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


  		




  
  38
  
    
## drop missing values


  		




  44
  39
  
    
cld$Amt[cld$Amt<0]=NA


  		




  45
  40
  
    
cld$Fq[cld$Fq<0]=NA


  		




  46
  41
  
    
cld$AWP[cld$AWP<0]=NA


  		




  47
  42
  
    
cld$NC[cld$NC<0]=NA


  		




  48
  43
  
    

  		




  49
  
  
    
cld$col=cut(cld$Amt/100,quantile(cld$Amt/100,seq(0,1,len=5),na.rm=T))


  		




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


  		




  50
  61
  
    

  		




  51
  62
  
    
library(lattice)


  		




  52
  
  
    
xyplot(lat~lon|as.factor(YR)+as.factor(month),groups=col,data=cld,pch=16,cex=.2,auto.key=T)


  		




  
  63
  
    
xyplot(lat~lon|+month,groups=col,data=cldm,pch=16,cex=.2,auto.key=T)


  		




  
  64
  
    

  		












Also available in:  Unified diff