This post shows an R script to automatically generate a trend chart from web based global temperature data. The R script allows me to update my plot each month as soon as the source file is updated. The plot is self documenting so that the chart reader can see the data period, the overall trend line, slope for trend line; the last data point is highlighted and value given. Date stamp and name are included in the margin.
There are numerous climate data sites that provide continuously updated data files for all aspects of climatology. While many sites provide useful charts of their data, I find that I want to combine data from several files and/or data sites, so I need the ability to download data files and generate my own climate charts. When I started my climate change learning expedition a few years ago, I thought that I’d be able to do everything in Excel. As I progressed in my efforts, I learned a few things:
- There’s an incredible amount of climate data available on the web
- Excel’s data and charting tools are limited for multi-variate analysis and charting
- It’s not realistic to store climate data on my PC, it’s more effective to use the web based climate data
- I need tools that let me automatically access the latest data from several data sites and produce my own charts
Along the way I realized that I needed to add R to my data analysis toolkit if I wanted to continue my study of climate trends. So I’ve been working my way up the R learning curve. While I have a great deal more to learn about R, I’m now at the point where I can resume my climate change data analysis and charting interests using R.
Many visitors to my ProcessTrends site download my trend analysis and global warming Excel workbooks. My goal is to provide comparable R tools so that visitors can download R tools. This post starts my effort to provide R tools for users to access and analyze web based climate data similar to my Excel tools.
Remote Sensing System (RSS) Temperature Anomaly Trends
RSS provides monthly satellite based temperature measurement and anomaly data files by various zones. Here’s a sample from the February 4, 2009 land and ocean temperature anomaly data file.
This file snippet shows monthly temperature anomaly data for 9 zones across the globe. The file is updated monthly. The data series, started in January, 1979, now has 30 years of measurements.
My goal is to be able to download these types of file each month and generate trend charts like the one shown in Figure 2. This chart includes the data series, a trend line based on a simple linear model, the slope of the trend line, and plot annotation to fully document the chart.
RSS’s 1979-2009 temperature anomaly trends show a positive trend of 0.0157 C per year, with significant up and down swings. In a previous post I compared El Nino – La Nina event and volcanic eruption time-lines to the GISS temperature anomalies. With R’s lattice plotting capabilities, I can show multiple time series in a single plot.
Now that I can automatically download and analyze web based data, I want to be able to compare climate change trends by zone, land versus ocean and RSS versus GISS. I’ll then move on to regression models to see how much effect El Nino – La Nina has on temperature trends.
Since I plan to continue exploring climate data, I have tried to develop my R script so that it is reusable from project to project.
R Script For Web File Download & Charting
Here’s the R script listing.
# RSS MSU Satellite Data - TLT - Temp Anomaly Trend - Lower Troposhpere - Land & Ocean C ## R script by D Kelly O'Day, 2/5/09 ## Get Monthly anomaly data link1<- "http://www.remss.com/data/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_TLT_Anomalies_Land_and_Ocean_v03_2.txt" mo_RSS_in <- read.table(link1, skip = 3, sep = "", dec=".", row.names = NULL, header = FALSE, as.is = T, colClasses = c(rep("numeric",3), rep("NULL", 8)), comment.char = "#", na.strings = c("*", "-",-99.9, -999.9), col.names = c("yr", "mo", "RSS_lo", rep("",8))) yr_frac <- mo_RSS_in$yr + (mo_RSS_in$mo-1)/12 # yr_frac simplifies calcs # Create new data frame & attach mo_RSS_lo <- data.frame(mo_RSS_in, yr_frac) attach(mo_RSS_lo) # Get last mo, yr and value c <- nrow(mo_RSS_lo) # Find number of data rows lmo <- mo_RSS_lo$mo[c] # Fimd last month of data lyr <- mo_RSS_lo$yr[c] # Find last year of data lRSS <- RSS_lo[c] # Find last reading # Create plot # Set plot parameters par(las=1); par(ps = c(10)) ; par(bty="l") par(oma=c(0,0,0,0)) # no outer margin plot(RSS_lo ~ yr_frac, type="l", col ="grey60", xlab = "", ylab = expression(paste("RSS Temperature Anomaly - ",degree,"C")), xlim=c(1979, 2010), ylim=c(-0.6, 0.80), xaxs="i", yaxs="i", cex.axis = 0.95, cex.lab = 0.95) abline(h=0, col = "grey") # Calculate Linear Model for rate of change lm_fit <- lm(RSS_lo ~ yr_frac) a <- coef(lm_fit) b <- coef(lm_fit) ## Determine min & max yr_frac yr1 <- min(yr_frac) yr2 <- max(yr_frac) y1 <- a + b* yr1 y2 <- a + b *yr2 x_val <- c(yr1, yr2) y_val <- c(y1, y2) lines(x_val,y_val, type = "l", col = "red") b <- signif(b, 3) # Plot annotation n_text <- 0.9 # use to control note font size # Establish plot titles title_3 <- expression(paste("RSS Land & Ocean Temperature Anomalies - ",degree, "C")) title_2 <- paste("Period: 1/1979 to ",lmo,"/",lyr) mtext(title_3, side = 3, line = 3, adj = 0.5) mtext(title_2, side = 3, line = 2, adj = 0.5, cex = n_text) mtext( "RSS Version 3.2 - MSU/AMSU TLT", side = 3,line = 1, adj=.5, cex = n_text) # Establish plot notes text(1995,-0.4, "Overall slope = ", pos = 1, col = "red", cex = n_text) text(1998.7, -0.4, b, pos = 1, col ="red", cex = n_text) text(2002,-0.4,expression(paste(degree,"C/yr")), pos = 1,cex = n_text, col = "red") text(1980, 0.7, " Area: 70S to 82.5N", col = "black", pos =4, cex = n_text) # Highlight last reading points( yr2, lRSS, pch=19, col = "black") points(1996, -0.5, pch=19) note <- paste(lmo, "/", lyr, " @ ", lRSS," C") text(1996, -0.5, note, pos = 4, col = "black", cex = n_text) # Add margin notes for name and plot date mtext("K O'Day, ChartsGraphs.Wordpress.Com",side=1,line = 3,adj = 0,cex = 0.7) mtext(format(Sys.time(), "%m/%d/ %Y"), side = 1, line = 3, adj = 1, cex = 0.7) detach(mo_RSS_lo)