Arctic Sea Ice Extent Trends With R

This post shows how to retrieve on-line Arctic Sea Ice Extent data from the National Snow and Ice Data Center (NSIDC), consolidate the data files, generate a csv file, summarize and plot the data and post it as a Google Docs so that interested readers can download and analyze this data series themselves.

Links to the NSIDC source data files, R script and my Google Docs csv file are provided.

Sea Ice Extent Data and Definition

The Japan Aerospace Exploration Agency (JAXA) provides daily updates to the Arctic Sea Ice Extent and maintains a downloadable csv  file of daily values from June, 2002 at this link.

JAXA and other  climate science data sites define  sea ice extent as the area of the  ocean where sea-ice concentration is 15% of more.

NSIDC maintains a series of 12 ASCII files (one for each month) that include the monthly average sea ice extent values for the period 1979 to present. I am using the NSIDC data set because it covers a longer period (31 years) than the JAXA data series.

Here’s the 1979 – 2009 trend chart of monthly Arctic sea ice extent measurement data based on NSIDC files. (Click image to enlarge)

This chart shows the monthly sea ice extent measurements from January, 1979 to December, 2009, a 31 year period.  The annual minimum and maximum monthly values are highlighted in red and blue respectively.  The annual average value is also included as well as the linear trend lines for the minimum and maximums.

Visual inspection, the trend in the minimums and maximums all show a downward trend in sea ice extent over the 31 year period.

I have posted a monthly csv file on the 1979-2009 NSIDC Arctic sea ice extent and area data on Google Docs so that you can view the data here or download a csv file here.

R Script to Consolidate NSIDC Files

Here’s  my R script to download the 12 monthly files, process the data and generate the csv file that I have posted on Google Docs. Note that I have a 12 line file csv file (view here) for the links to the Monthly NSIDC files.

#### NSIDC Arctic Sea Ice Extent Chart     ###############################################
##                       NSIDC_trend_plot.R                                             ##
##                                                                                      ##
## Reads & plots monthly NSIDC data file developed by DK O'Day                          ##
## and stored on GoogleDocs as spreadsheet to simplify user access                      ##
##  D Kelly O'Day         http://chartsgraphs.wordpress.com     12/30/09                ##
## Data link: http://spreadsheets.google.com/pub?key=tNZHAYnHUvPbyGjCfj2SJ9g&output=csv ##
##########################################################################################
  par(las=1); par(ps=10); par(oma=c(3.5,2,1,1)); par(mar=c(2,4,2,1))
  par(bg="white")
  library(zoo)

########### NSIDC Monthly Artcic Sea Ice Data: 1979 - Present       ######################
link_n <- "http://spreadsheets.google.com/pub?key=tNZHAYnHUvPbyGjCfj2SJ9g&output=csv"
 sie <- read.table(link_n, skip = 1, sep = ",", dec=".",
             row.names = NULL, header = FALSE,
             as.is = T,
             colClasses = rep("numeric", 5),
             comment.char = "#", na.strings = "NA",
             col.names = c("y_fr", "yr", "mo", "ext", "area")   )
  sie_data <- sie[,c(1,3,4)]
  sie_data <- subset(sie_data, y_fr >= 1979) # get 1st full year of data
  max_yr <- (trunc(max(sie_data$y_fr, na.rm=T)))
  min_yr <- trunc(min(sie_data$y_fr, na.rm=T))
  n_yrs <- as.integer(max_yr - (min_yr-1))
## Calc min & max values & month in each year
# Set up arrays
   min_ext <- numeric(n_yrs)
   min_yr_frac <- numeric(n_yrs)
   min_mo <- numeric(n_yrs)
   avg_ext <- numeric(n_yrs)
   avg_yr_frac <- numeric(n_yrs)
   max_ext <- numeric(n_yrs)
   max_yr_frac <- numeric(n_yrs)
   max_mo <- numeric(n_yrs)
## loop through all years
 for (i in 1:n_yrs) {
   temp <- subset(sie_data,trunc(sie_data$y_fr) == i+1978)
     min_ext[i] <- min(temp$ext, na.rm=T)
     row_min <- which.min(temp$ext)
     min_yr_frac[i] <- temp[row_min,1]
     min_mo[i] <- temp[row_min,2]
     avg_ext[i] <- mean(temp$ext, rm.na=T)
     avg_yr_frac[i] <- i + 1978
     max_ext[i] <- max(temp$ext, na.rm=T)
     row_max <- which.max(temp$ext)
     max_yr_frac[i] <- temp[row_max,1]
     max_mo[i] <- temp[row_max,2]
    }
  mins <- data.frame(min_yr_frac,min_mo, min_ext)
  maxs <- data.frame(max_yr_frac, max_mo, max_ext)
  avgs <- data.frame(avg_yr_frac, avg_ext)
  overall_avg <- mean(avgs$avg_ext, na.rm=T)
## Trend Rates
  lm_min <- lm(mins$min_ext ~ mins$min_yr_frac)
  min_sum <- summary(lm_min)
  min_r2 <-  min_sum$adj.r.squared
   a_min <- coef(lm_min)[1]
   b_min <- coef(lm_min)[2]
  lm_max <- lm(maxs$max_ext ~ maxs$max_yr_frac)
  max_sum <- summary(lm_max)
  max_r2 <-  max_sum$adj.r.squared
   a_max <- coef(lm_max)[1]
   b_max <- coef(lm_max)[2]
 cat(c("Max r2= ", signif(max_r2,3), "Min r2= ", signif(min_r2,3)), fill = 20)
## Generate Plot
  title = "Arctic Sea Ice Extent Trend\nBased on NSIDC Monthly Data: 1979-2009"
  y_title <- expression(paste("Monthly Arctic Sea Ice Extent - Millions k",m^2))
  plot(sie_data$y_fr, sie_data$ext, col = "grey", type = "l",
   xlim = c(1979, 2010), xlab = "", ylim=c(0, 16.5), ylab = y_title, main = title,
   cex.main=0.9)
  abline(h=overall_avg, col = "darkgrey", lty="solid")
  points(avgs$avg_yr_frac, avgs$avg_ext, type="S", col = "black")
  points(mins$min_yr_frac, mins$min_ext, type = "p",pch=19, col = "red", cex=0.5)
  points(maxs$max_yr_frac, maxs$max_ext, type = "p", pch=19,col = "blue", cex=0.5)

  abline(a=a_min, b=b_min, col = "red")
  abline(a=a_max, b=b_max, col = "blue")
  tr_min <- paste("Ann Min Regr Trend = ", as.character(signif(b_min,2)))
  tr_max <- paste("Ann Max Regr Trend = ", as.character(signif(b_max,2)))
  text(1992,5 ,tr_min, adj=0, cex=0.8, col = "red")
  text(1992,16.5, tr_max, adj=0, cex=0.8, col = "blue", bg="red")
  overall_note <- paste("Overall avg - ", as.character(signif(overall_avg,3)))
  rect(1978,10.7, 1983,11.3, col="white", border="white")
  text(1978, 11, overall_note, adj=0,cex=0.7, col = "black")
  note<- "See post for link to on-line file\nof NSIDC monthly time series available\nfor download as CSV file from\nGoogle Docs"
  text( 1995, 1.5, note, cex = 0.7, col ="darkgrey", adj = 0)
  legend(1978.5,4, c("Monthly SIE(1979 to present)", "Ann Max", "Ann Max Regr", "Ann Mean Trend", "Ann Min", "Ann Min Regr"),
               col = c("grey","blue","blue", "black", "red","red"),
               text.col = "black", lty = c(1,0,1,1,0,1), pch=c(1,19,1,1,19,1),pt.cex=c(0,1,0,0,1,0),
               merge = F, bg = NULL, bty="o", box.col="white",cex = .65)

## Plot Annotations

  source_1 <- "NSIDC Data Source files: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/"
  source_2 <- "Climate Charts & Graphs time series file:"
  source_3 <- link_n
  mtext(source_1, 1,0,  outer=T, adj=0.5, cex=0.7)
  mtext(source_2,1, .75, outer=T, adj=0.5, cex = 0.7)
  mtext(source_3, 1,1.5, outer=T, adj=0.5,cex = 0.7)
  mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,2.25, adj = 0, cex = 0.8, outer=TRUE)
  mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 2.25, adj = 1, cex = 0.8, outer=TRUE)

About these ads