### R code for my Seasonal Trends graph

I had a request for the code I used to generate the graphs in my Seasonal Trends post.

While it looks complex, the R code for it is very simple IF you have the data ready. I'm assuming that you already have the temperature dataset you want as an R object (I have several datasets in an object I simply call "monthly": GISS, Berkeley Earth, Cowtan-Way, HadCRUT4, UAH, and RSS, along with the year/decimal month Time, Year, and numeric Month). The code I used to generate the graph is as follows:

While it looks complex, the R code for it is very simple IF you have the data ready. I'm assuming that you already have the temperature dataset you want as an R object (I have several datasets in an object I simply call "monthly": GISS, Berkeley Earth, Cowtan-Way, HadCRUT4, UAH, and RSS, along with the year/decimal month Time, Year, and numeric Month). The code I used to generate the graph is as follows:

#Create separate datasets for each season

monthly.79<-subset(monthly, Time>=1978.92 & Time<2013.91)

DJF<-subset(monthly.79, Month=="12" | Month =="1" | Month=="2")

DJF$Year_2 <- numeric (length (DJF$Year))

for (i in 1:length (DJF$Year) ) {

if ( DJF$Month [i] == 12) {

DJF$Year_2[i] <- DJF$Year [i] + 1

}

else {

DJF$Year_2[i] <- DJF$Year [i]

}

}

MAM<-subset(monthly.79, Month=="3" | Month =="4" | Month=="5")

JJA<-subset(monthly.79, Month=="6" | Month =="7" | Month=="8")

SON<-subset(monthly.79, Month=="9" | Month=="10" | Month=="11")

#Calculate the seasonal average for each year

DJF<-aggregate(DJF$BEST, by=list(DJF$Year_2), FUN=mean)

MAM<-aggregate(MAM$BEST, by=list(MAM$Year), FUN=mean)

JJA<-aggregate(JJA$BEST, by=list(JJA$Year), FUN=mean)

SON<-aggregate(SON$BEST, by=list(SON$Year), FUN=mean)

#Check for autoregression

library(forecast) #for the auto.arima function

auto.arima(resid(lm(x~Group.1, data=DJF)), ic=c("bic"))

auto.arima(resid(lm(x~Group.1, data=MAM)), ic=c("bic"))

auto.arima(resid(lm(x~Group.1, data=JJA)), ic=c("bic"))

auto.arima(resid(lm(x~Group.1, data=SON)), ic=c("bic"))

#Construct the plot

plot(x~Group.1, data=DJF, type="l", col="blue", xlab="Year", ylab="Temperature anomaly (ºC)", main="Seasonal Climate Trends", ylim=c(-0.1, 0.8))

points(x~Group.1, data=MAM, type="l", col="green")

points(x~Group.1, data=JJA, type="l", col="red")

points(x~Group.1, data=SON, type="l", col="orange")

#Add the trend lines

abline(lm(x~Group.1, data=DJF), col="blue", lwd=2)

abline(lm(x~Group.1, data=MAM), col="green", lwd=2)

abline(lm(x~Group.1, data=JJA), col="red", lwd=2)

abline(lm(x~Group.1, data=SON), col="orange", lwd=2)

legend(1979, 0.8, c("DJF", "MAM", "JJA", "SON"), col=c("blue", "green", "red", "orange"), lwd=2)

#Get the slopesThat's all there was to it. I just repeated this code, modifying only the period of the initial subset to create the second graph (Monthly.2002<-subset(monthly, Time>=2001.92 & Time<2013.91) and the related seasonal subsets. To the person who requested my code: Hope this helps.

summary(lm(x~Group.1, data=DJF)

summary(lm(x~Group.1, data=MAM)

summary(lm(x~Group.1, data=JJA)

summary(lm(x~Group.1, data=SON)

## Comments

## Post a Comment