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

[This article was first published on Climate Charts & Graphs I » RClimate Script, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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!!

library(animation)
  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 <- "http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv"
  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",
       main=title)
    text(20, 1.5, note_1, cex = 0.8, adj=0, col = "black")
    text(20,1,"Data Source: http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv", cex = 0.8, adj=0,col = "black")
    mtext("D Kelly O'Day - https://chartsgraphs.wordpress.com", 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") dev.off()
}
## 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")
}

Filed under: Arctic Sea Ice, Do-it-yourself Climate Science, Global Warming, Learn R, R Example and Scripts, RClimate Script, Time Series Charts Tagged: Arctic SeaIce, Climate Trends, R scripts

To leave a comment for the author, please follow the link and comment on their blog: Climate Charts & Graphs I » RClimate Script.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)