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:

#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 slopes

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)

 That'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.