[R] climatological standard deviation- (question re-posted)

[Daniel Nordlund]

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of Yogesh Tiwari
> Sent: Wednesday, September 24, 2008 10:55 PM
> To: r-help at stat.math.ethz.ch; sarah.goslee at gmail.com; ssefick at gmail.com;
> milton.ruser at gmail.com
> Subject: Re: [R] climatological standard deviation- (question re-posted)
> 
> Dear R Users,
> Sorry for not explaining the problem clearly.
> Here we go.......
> I am using R on windows OS....
> Here is an example data and script and then the problem what I want to
> solve.
> 
> # (data period 1993-2002)
> #example data
> yr mo co2
> 1993 2 359.543
> 1993 2 358.359
> 1993 2 359.315
> 1993 2 359.293
> 1993 4 362.756
> 1993 4 363.699
> 1993 4 363.505
> 1993 4 363.034
> 1993 6 358.482
> 1993 6 359.283
> 1993 7 356.739
> 1993 7 357.203
> 1993 7 357.257
> 1993 7 357.062
> 1993 8 357.618
> 1993 8 356.965
> 1993 8 356.11
> 1993 8 356.389
> 1993 9 356.418
> 1993 9 356.274
> 1993 9 355.7
> 1993 9 357.389
> 1993 10 351.972
> 1993 10 352.132
> 1993 10 352.181
> 1993 10 352.015
> 1993 10 355.774
> 1993 10 348.42
> 1993 10 348.376
> 1993 10 348.599
> 1993 11 349.954
> 1993 11 350.575
> 1993 11 350.301
> 1993 12 354.179
> 1993 12 353.917
> 1993 12 356.546
> 1994 2 366.965
> 1994 2 367.716
> 1994 2 367.711
> 1994 2 368.595
> 1994 3 367.167
> 1994 3 367.732
> 1994 3 368.427
> 1994 3 368.547
> 1994 3 367.833
> 1994 3 368.423
> 1994 4 361.591
> 1994 4 361.682
> 1994 5 362.952
> 1994 5 362.629
> 1994 5 360.703
> 1994 5 360.813
> 1994 5 366.137
> 1994 5 364.946
> 1994 6 357.072
> 1994 6 357.351
> 1994 7 357.888
> 1994 7 361.456
> 1994 7 357.311
> 1994 7 357.503
> 1994 8 356.242
> 1994 8 356.184
> 1994 8 357.733
> 1994 8 356.729
> 1994 9 357.147
> 1994 9 356.34
> 1994 9 357.626
> 1994 10 354.157
> 1994 10 353.994
> 1994 10 365.504
> 1994 10 365.036
> 1994 11 355.671
> 1994 11 355.209
> 1994 11 352.027
> 1994 11 351.612
> 1994 12 357.318
> 1994 12 357.762
> 1994 12 352.377
> 1994 12 352.719
> 1994 12 355.844
> 1994 12 354.784
> 1995 1 363.532
> 1995 1 364.478
> 1995 1 358.964
> 1995 1 359.402
> 1995 2 363.664
> 1995 2 363.746
> 1995 2 365.041
> 1995 2 364.56
> 1995 3 368.469
> 1995 3 367.798
> 1995 3 362.802
> 1995 3 363.027
> 1995 4 363.698
> 1995 4 363.772
> 1995 4 364.125
> 1995 4 364.408
> 1995 5 362.638
> 1995 5 362.534
> 1995 5 362.701
> 1995 6 366.4
> 1995 6 362.5
> 1995 6 359.634
> 1995 6 360.256
> ....
> ....
> #Script to compute monthly climatology (12 values, Jan-Dec),
> file<-read.csv("CRI-run109-run110.csv")
> names(file)
> attach(file)
> co2.month.mean <- tapply(co2,mo,mean)
> x<- as.integer(row.names(co2.month.mean),co2.month.mean)
> plot(x, co2.month.mean, xlim=c(1,12),
> ylim=c(358,372),xlab=NA,ylab=NA,xaxs="i",yaxs="i",col="grey",type="o",axes=F)
> #
> axis(1,at=1:12,labels=c("J","F","M","A","M","J","J","A","S","O","N","D"))
> axis(2)
> axis(3, at=1:12, labels=FALSE)
> axis(4, at=NULL, labels=FALSE)
> #standard deviation calulation for each month
> #....jan
> time_span <- file[(file$mo==1),]
> sub_co2 <- (time_span$co2)
>  var.jan<-var(sub_co2, na.rm=TRUE)
>  sd.jan<-sqrt(var.jan)
> #....feb
>  time_span <- file[(file$mo==2),]
> sub_co2 <- (time_span$co2)
>  var.feb<-var(sub_co2, na.rm=TRUE)
>  sd.feb<-sqrt(var.feb)
> #....mar
> ...
> 
> #Problem:
> Need to compute climatological standard deviation i.e the 12 values
> calculated from the above data january values, february values, etc.
> and need to plot as an error bar on the ploted curve.
> 
> Is this the best way for such a calculation ?
> 
> Kindly help,
> 
> Best Regards,
> Yogesh
> 
You might find the errbar() function from the Hmisc package useful.  Here is some code to get you started.

#load Hmisc- you may need to install the package
library(Hmisc)

attach(file)
co2.mean <- aggregate(co2,by=list(mo),FUN=mean)
co2.sd <- aggregate(co2,by=list(mo),FUN=sd)
names(co2.mean) <- c('month','mean')
names(co2.sd) <- c('month','sd')
co2.stats <- merge(co2.mean,co2.sd)
with(data=co2.stats, errbar( month, mean, mean+sd, mean-sd ))


Hope this is helpful,

Dan

Daniel Nordlund
Bothell, WA USA