R code for simple RSS graph

CombinedUAHRSSHadCrutGISSAMOPDOCET

Now available as a pdf file at the code page


#RLH Simple RSS plot - Mar 2014

#required libraries
require(stats)
require(signal)

#simple moving average FIR filter
MovingAverage <- function(x,n=5)
{
stats::filter(x,rep(1/n,n), sides=2)
}

#Cascaded Triple Running Mean 'Gaussian' FIR filter on data using 1.2067 inter stage multiplier from V. Pratt
CTRM <- function(data, period=12)
{
f1 = period
f2 = round(f1/1.2067)
f3 = round(f2/1.2067)
CTRM = MovingAverage(data,f1)
CTRM = MovingAverage(CTRM,f2)
CTRM = MovingAverage(CTRM,f3)
}

#"I ran a 5 pass-multipass with second order polynomials on 15 year data windows as per the Savitzky-Golay method." Nate Drake PhD
SavitzkyGolay <- function(data, period=12)
{
f1 = period * 2 + 1
SavitzkyGolay = signal::sgolayfilt(data,n=f1)
SavitzkyGolay = signal::sgolayfilt(SavitzkyGolay,n=f1)
SavitzkyGolay = signal::sgolayfilt(SavitzkyGolay,n=f1)
SavitzkyGolay = signal::sgolayfilt(SavitzkyGolay,n=f1)
SavitzkyGolay = signal::sgolayfilt(SavitzkyGolay,n=f1)
}

#RSS global data url from the 'net
datasource = "http://data.remss.com/msu/graphics/TLT/time_series/RSS_TS_channel_TLT_Global_Land_And_Sea_v03_3.txt&quot;

#get the data
RSS = read.table(datasource, skip=6, col.names=c("Year","Month","Anomalies"))

#remove the invalid monthly values
RSS = subset(RSS, Anomalies != -99.9)

#get/save the last month name as a string
month = month.abb[as.numeric(tail(RSS[2],1))]

#convert the date columns to a fractional value in col 1
RSS[1] = RSS[,1] + (RSS[,2]/12) - 1/24
RSS[2] = NULL

#make up the title text
main = sprintf("%s RSS Global with Annual CTRM low pass filter and trend", month)

#change to TRUE for plot to file
PlotToFile = FALSE

#change the directory, size, etc. as required
if(PlotToFile)
{
png(filename=sprintf("C:/R/Images/%s.png", main), width=1500, height=700)
}

plot(RSS, main=main) # main plot
abline(h=(seq(-1,1,0.1)), col="gray", lty=2) # add the grid
abline(v=(seq(1975,2020,1)), col="gray", lty=2)
lines(RSS[,1], SavitzkyGolay(RSS[,2], 15*12), col="red", lwd=1, lty=2) # S-G 15 year trend
lines(RSS[,1], CTRM(RSS[,2]), col="green", lwd=3) # Annnual summary
points(tail(RSS[,1], 1), tail(RSS[,2], 1), col="green", pch = 19) # highlight the current value
legend('bottomright', c("Annual LP","Trend"), lty=c(1,2), lwd=c(3,1), col=c("green","red"))
mtext(side=1, line=3, adj=1, format(Sys.Date(), format="RLH %b %Y"))

#close the file if required
if(PlotToFile)
{
dev.off()
}

Advertisements

7 thoughts on “R code for simple RSS graph

      • I hope so! The link below leads to the images from the code, one for the RSS data, the other for HadCrut4. They compare very closely to yours, although the phase of the CTRM data may be slightly off. This is due to the same reason as your Excel example from WUWT, which is the selection of the starting index for the running mean.

        Image album: http://imgur.com/a/aEOtq

  1. I wouldn’t worry about that small phase offset. It is almost always of no concern. I do always try and plot the original data as points though (as you have done) to ensure that the context of the output curves are preserved.

    Helps keep some relevance to the output. It does hide some, possibly relevant, information if you omit them.

  2. Richard, I played with your filter using the HADCRU hemispheric datasets.
    We see a lovely sine wave in the SST, with a background rise.
    We get the same with the NH land.
    If you do a simple fit, using log (CO2), made into months, and a sine-wave, you get nice enough fits.
    So far, so boring.
    However, one cannot but help but notice that the northern hemisphere land and oceans sine waves have the same periodicity and very similar amplitudes, but the land leads.
    The NH land changes before the oceans.
    Could you have a little look and see that I haven’t lost my remaining marble.
    If the AMO does follow a land temperature change, the stadium wave hypothesis suddenly become a big drive.

    • Hi Doc,

      Well there are lots of interesting things to look at. A full treatment of HadCRUT is on my list of things to do. If you look on the HadCRUT page now you will see the ‘cycle’ plot I did for WUWT. A more detailed look into HadCRUT will follow shortly.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s