Category Archives: R Example and Scripts

R Script to Build Animation of Arctic Sea Ice Extent – Update 12/20/13

In my previous post I showed an animation of Arctic Sea Ice Extent from the 1980’s through August, 2012 (link).  In this post, I show how to build this Arctic Sea ice Extent  animated chart.

Source Data

The Arctic Ice Sea Monitor (link)   updates their daily csv file with the latest satellite based arctic sea ice measurements.  Here is the daily csv file link.

R script

To develop my animation of the daily Arctic Sea Ice extent, I decided to produce a plot for each year that showed the current year in red and the previous years in grey.  I go this idea from Tamino at Open Mind.

Here is my R script:
Be sure to set your working directory to appropriate location!!

  ani.options(convert=shQuote('C:\\Program Files (x86)\\ImageMagick-6.7.9-Q16\\convert.exe'))
## Use setwd() to specify directory where you want png images to be saved
  setwd("<strong>C:\\R_Home\\Charts & Graphs Blog\\RClimateTools\\Arctic_sea-ice_extent</strong><em>")
# use png_yn to toggle between plot output to png file or screen
  png_yn <- "y"
# Establish chart series patterns and colors to be able to distinguish current yr from previous years in plot
  pattern <- c(rep("dashed", 5), rep("solid", 12))
  ser_col <- c(rep("black",5),rep("grey",12))
# Establish chart annotations for date, chart title,
  what_date <- format(Sys.Date(), "%b %d, %Y")  # with month as a word
  title <- paste("IARC-JAXA Daily Arctic Sea Ice Extent*\n", what_date)
  note_1 <- "*Extent - Area of Ocean with at least 15% Sea Ice"
  par(oma=c(2,1,1,1)); par(mar=c(2,4,2,1))
#  Day of year axis setup
## Set up basic day of year vectors (mon_names, 1st day of mon)
  mon_names <- c("Jan", "Feb", "Mar", "April", "May", "June", "July", "Aug", "Sept", "Oct","Nov","Dec")
  mon_doy <- c(1,32,60,91,121,151,182,213,244,274,305,335,366)
  mon_pos <- c(16, 46, 75, 106,135, 165, 200, 228, 255, 289, 320, 355)
# Read JAXA Arctic Sea ice Extent csv file
# Data File: Month,Day,1980's Avg,1990's Avg,2000's Average,2002:2012
  link <- ""
  j_data <- read.csv(link, header = F, skip=1, na.strings = c(-9999))
 series_id <-  c("mo", "day", "1980s", "1990s", "2000s","2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009",
                "2010", "2011", "2012", "2013")
 colnames(j_data) <- series_id
# File has data for each day in 366 day year
# Establish Day of year
  for (i in 1:366)   j_data$yr_frac[i] <- i
    #convert ASIE to millions Km^2
   j_data[,c(3:17)] <- j_data[,c(3:17)]/1000000
# Loop through years
   for (j in 3:17)
     png_name <- paste("asie",series_id[j],".png",sep="")
      if (png_yn =="y") png(filename=png_name)
      which_yr <- j
      no_yrs <- j
  # Calc min asie for year
    min_asie <- min(j_data[,j], na.rm = T)  # must remove na's to get valid answer
    lab_asie <- round(min_asie,3)
    min_r <- which(j_data[,j] == min_asie)
    min_d <- j_data[min_r,2]
    min_m <- j_data[min_r,1]
    min_date <- paste(min_m,"/",min_d,"/",series_id[j], sep="")
    plot(j_data[,17],  type="n", col = "grey",axes=F, xlab="",
       ylab="Arctic Sea Ice Extent - Millions Sq KM",
       ylim=c(0,15),xaxs="i", yaxs = "i",
    text(20, 1.5, note_1, cex = 0.8, adj=0, col = "black")
    text(20,1,"Data Source:", cex = 0.8, adj=0,col = "black")
    mtext("D Kelly O'Day -", 1,0.5, adj = 0, cex = 0.8, outer=T)
  # custom x & y axes
    axis(side = 1, at=mon_doy, labels=F, xaxs="i")
    axis(side=1, at= mon_pos, labels=mon_names, tick=F, line=F, xaxs="i")
    axis(side=2,  yaxs="i", las=1)
    points(70, min_asie, col = "red",pch=19, cex = 2)
  # Add each previous yr data series as light grey line
  for (n in 3:no_yrs)
    points(j_data[,18], j_data[,n], type="l",lwd=1,lty=pattern[j], col=ser_col[j])
    text(182,14,series_id[j], col = "red", cex = 1.1)
  points(j_data[,18], j_data[,j], col="red", type="l",lwd=2.5)
  text(182,14,series_id[j], col = "red", cex = 1.1)
  text(120,min_asie+0.5, min_date, col="red", cex=0.9)
  text(120,min_asie, lab_asie, col="red", cex=0.9)
  if(png_yn == "y")
## copy last png file 3 times to provide pause in animation
if(png_yn== "y")
  for (c in 1:2)
    file_name <- paste("asie2012",c, ".png",sep="")
    file.copy(from= "asie2012.png", to = file_name, overwrite=T)
  ani.options(outdir = getwd())    # direct gif output file to working dir
  ani.options(interval= 0.80)
  im.convert("asie*.png", "last_animation.gif")

Do It Yourself Climate Trend Analysis

This post describes my consolidated global temperature anomaly CSV file that users can easily download to Excel or R to do their own trend analysis.

Do It Yourself Global Temperature Anomaly  Trend Analysis

As I wrote in my July 10, 2009 post,

“There are many blogs and web sites (small sample: 1, 2, 3, 4, 5, 6) with multiple opinions on global climate trends. Some sites are data oriented and others are opinion oriented. What is a [data analyst] charter to think?

My advice, take a look at the data for yourself. As an Excel or R charter, why not analyze it yourself to get a better appreciation for what is going on.

To help you get started, I’ve developed a consolidated monthly CSV file that presents  the 5 major global land and ocean temperature anomaly data series: GISS, NOAA, HADCrut3, RSS and UAH.

Here’s the link to my consolidated temperature anomaly CSV file.

I update the consolidated file regularly by downloading the latest agency source files so that the consolidated file includes all source agency data revisions. This way you can get the most up-to-date temperature anomaly data without having to reformat/ manipulate the 5 individual files.

Tracking Long Term and Recent UAH Channel 5 Anomally Trends

In this post I present  combined long term and recent trend charts of UAH Channel 5 temperature anomalies using R’s figure inside figure capability.  This approach provides a better picture of global temperature trends than UAH’s day-of-year plots.


The University of Alabama Huntsville (UAH) provides a daily update to their Channel 5 global average temperature at 14,000 feet. I have previously posted about this data set here. The source data file link is here.

Continue reading

RClimate Script: Global Mean Sea Level Trend by Month

This RClimate Script lets users read a NetCDF file and plot the latest satellite altimetry based global mean sea level data  from 1993 to the latest completed month. The trend chart shows NOAA’s Laboratory for Satellite Altimetry altimeter global mean sea level with the seasonal signal removed and inverted barometer.

I will add this chart to my Climate Trend update sidebar and update it each month as NOAA releases updated mean sea level data.

Continue reading

RClimate Script: Sea Surface Temperature (SST) Anomaly Trends

This RClimate Script lets users retrieve and plot the National Climatic Data Center’s monthly Sea Surface Temperature (SST) Anomaly dat series .  Links to the NCDC data file and my RClimate script are included. Users can run my script with a simple R source() statement.

NCDC  SST Anomaly

I’ve discussed the Hadley SST anomaly trends in this earlier post.

Here’s NCDC’s global SST anomaly data series trend since 1880.

This chart shows the decadal means and highlights the latest monthly SST anomaly.

Here are the data and RClimate Script links:

RClimate Script: NINO 3.4 SST Anomaly Trends

This RClimate Script lets users retrieve and plot the weekly NOAA NINO 3.4 SST  anomaly data for 1990 to the most recent value.  Links to the NOAA data file and my RClimate script are included. Users can run my script with a simple R source() statement.

NINO 3.4 SST Anomaly

I’ve discussed ENSO and NINO 3.4 in this earlier post.

SST’s in 4 equatorial Pacific zones are closely monitored to assess the status of El Nino – Southern Oscillation (ENSO).

NINO 3.4 SST anomaly provides a quick and effective indication of ENSO conditions. NOAA updates the NINO 1,2,3,3.4 and 4 SST and SSTA series weekly.

Here’s theweekly NINO 3.4 trend from 1990 to the most recent weekly reading. Click image to enlarge.

Here is the Data Link.

Here is the R script that I used to generate this chart.

## set link & read data
 #link <- ""
link <- ""
 nino_34 <- read.fwf(link, widths=c(10,-31,4,4,-17),skip = 4, header=F)
 names(nino_34) <- c("cdt", "sst", "ssta")

## calc yr_frac
 dt <- as.Date(nino_34$cdt, format="%d%b%Y")
 yr <- as.numeric(format(dt, format="%Y"))
 mo <- as.numeric(format(dt, format="%m"))
 dy <- as.numeric(format(dt, format="%d"))
 yr_frac <- as.numeric(yr + (mo-1)/12 + (dy/30)/12)
 nino_34 <- data.frame(dt,yr_frac, nino_34[,2:3])

## Determine Month, Year for last reading
 c <- nrow(nino_34) # Find number of data rows
 ldt <- nino_34[c,1]
 ldt_c <- paste(mo[c],"/", dy[c],"/", yr[c],sep="")
 l_yr_frac <- nino_34[c,2]
 l_nino_34 <- nino_34[c,4]
 lp_note <- paste(ldt_c, " @ ",l_nino_34,"C",sep="")

## subset data to add color for LaNina & El Nino cnditions
 la_nina <- subset(nino_34, ssta < 0)
 el_nino <- subset(nino_34, ssta>= 0)

## Plot titles
 title = "NINO 3.4 SST Anomaly Trend (Baseline: 1950-1979) \n NOAA: Weekly Data Centered on Wed"
 y_lab <- expression(paste("SST Anomaly ",degree,"C (1950-1979)", sep=""))

p_fun <- function() {
## set plot pars
 par(ps=10); par(las=1); par(oma=c(3.5,1,0,1)); par(mar=c(2,4,2,0))

plot(nino_34$yr_frac, nino_34$ssta, type = "n", main=title, xlab="", ylab = y_lab,
 cex.main=0.95, cex.lab=0.95)
 points(la_nina$yr_frac, la_nina$ssta, col = "blue", type = "h")
 points(el_nino$yr_frac, el_nino$ssta, col = "red", type = "h")
 points(l_yr_frac, l_nino_34, pch=19, col = "black",cex=0.75)
 points(1991.5, -2, pch=19, col="black", cex=0.75)
 text (1992, -2, lp_note, adj=0, cex=0.8)

## Generate and add bottom footer KOD, source, System date notes
 source_note <- paste("Data Source: ", link)
 mtext("D Kelly O'Day -", 1,1, adj = 0, cex = 0.8, outer=TRUE)
 mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.8, outer=TRUE)
 mtext(source_note, 1,0,adj=0.5, cex=0.8, outer=T)

RClimate Script: CO2 Trends

This RClimate Script lets users retrieve the latest data file on monthly Mauna Loa  CO2 levels and generate a trend chart with the latest reading highlighted.

CO2 TrendsKeeling Curve

Here’s the Mauna Loa Observatory CO2 trend from 1958 to Dec., 2009.

Here are the data and RClimate Script links: