Category Archives: RClimate Script

Climate Examples Using RClimate

Andy Long at Northern Kentucky University  has created a RClimate web page that provides numerous examples of climate data analysis using RWeb,  RClimate scripts and publicly available climate  data.

RClimate

NOAA’s Annual Global Temperature Anomaly Trends

Here is an alternative way using dygraphs (link)

NOAA has released their December, 2014 global anomaly data, allowing us to examine the 1880 – 2014 global temperature anomaly trend.

This R script trend charts shows the annual and average decadal temperature anomalies for the NOAA data series.

 

NOAA’s data, like NASA GISS’s and the Japan Meteorological Agency, show that the global mean annual temperature continues to rise.  2014 was the warmest year in the warmest decade for all 3 global temperature data series.

Here is my R script for those who would like to reproduce my chart:

 


############## RClimate Script: NOAA Annual Temperature Anomaly Trend     ##################
##   1/16/15     Retrieve daa from NOAA's and plot annual - decadal trends                ##
#############################################################################################
 library(plyr); library(reshape); library(RCurl)
 NOAA_ann_link <- "http://www.ncdc.noaa.gov/cag/time-series/global/globe/land_ocean/ytd/12/1880-2014.csv"
  n <- read.table(NOAA_ann_link, skip = 2, sep = ",", dec=".",
                 row.names = NULL, header = T,  as.is = T, colClasses = rep("numeric",2),
                 col.names = c("yr", "anom") )
## Find last report year and last anomaly value
  num_rows <- nrow(n)
  NOAA_last_yr <- as.integer(n[num_rows,1])
  NOAA_last_anom <- signif(n[nrow(n),2],2)
# Decade calculations
 dec_mean<- as.numeric(14)
 dec_st <- as.numeric(14)
 dec_end <- as.numeric(14)
 base_yr <- 1880
 n$dec_n <-  (as.numeric((n$yr - base_yr) %/% 10) * 10) + base_yr
# df <- data.frame(df, dec_n)
 for (i in 1:13) {dec_st[i] = base_yr+ i*10
                 dec_sub <- subset(n, dec_n == dec_st[i], na.rm=T)
                 dec_mean[i] <- mean(dec_sub$anom)
         }
 dec_st[14] <- 2020              # Need to have for last step line across decade
 dec_mean[14] <- dec_mean[13]
 dec<- data.frame(dec_st, dec_mean)
# Trend chart function
  plot_func<- function() {
   par(las=1); par(ps=12); par(oma=c(2.5,1,1,1)); par(mar=c(2.5,4,2,1))
   p_xmin <- 1880;   p_xmax <- n[num_rows,1]
   title <- paste("NOAA Land and Sea Temperature Annual Anomaly Trend \n", p_xmin, " to ", NOAA_last_yr, sep="")
   plot(n$yr, n$anom, type = "l", col = "grey",
     xlim = c(p_xmin, p_xmax), ylab = "Temperature Anomaly - \u00B0C (1951-1980 Baseline)",
     xlab="", main = title,cex.main = 0.85)
   points(NOAA_last_yr, NOAA_last_anom, col = "red", pch=19)
  last_pt <- paste( NOAA_last_yr, " @ ", NOAA_last_anom, " \u00B0C",sep="")
  points(dec$dec_st, dec$dec_mean, type="s", col="blue")
  ## add legend
  legend(1882,0.6, c("Decadal Avg Anomaly" ,"Annual anomaly", last_pt), col = c("blue", "grey", "red"),
       text.col = "black", lty = c(1,1,0),pch=c(0,0,16),pt.cex=c(0,0,1),
       merge = F, bg = "white", bty="o", cex = .75, box.col="white")
   out <- paste("NOAA Annual Temperature Anomaly \nData updated through: " , NOAA_last_yr, sep="")
   t_pos <- p_xmin + 0.5*(p_xmax-p_xmin)
  text(t_pos, -0.55, out, cex = 0.7, adj = 0)
  data_source <- paste("Data Source: ", NOAA_ann_link, sep="")
# Plot Annotation
 mtext(data_source,1,0, cex = 0.75, adj=0.5, outer=T)
 mtext("D Kelly O'Day - https://chartsgraphs.wordpress.com", 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)
  }
##########################################################################################
plot_func()

NASA GISS’s Annual Global Temperature Anomaly Trends

Update: Here is an alternative way to chart the same data (link).

NASA’s Goddard Institute for Space Studies (GISS)  has released their December, 2014 anomaly data, showing that 2014 was the warmest year in the warmest decade in the 1880 – 2014 instrumental temperature record period.

NASA’s results are consistent with the Japanese Meteorological Agency report (here)

 

Here is my R script for those who would like to reproduce my chart.

 

############## RClimate Script: GISS Annual Temperature Anomaly ###########################
##                   http:chartsgraphs.wordpress.com    1/16/15                            ##
############################################################################################
  library(plyr); library(reshape)
 ## File Download and File
  url <- c("http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt")
  file <- c("GLB.Ts+dSST.txt")
  download.file(url, file)

## 1st 8 rows and the last 12 rows contain instructions
## Find out the number of rows in the file, and exclude the last 12
    rows <- length(readLines(file)) - 12
## Read file as  char vector, one line per row, Exclude first 8 rows
    lines <- readLines(file, n=rows)[8:rows]
## Data Manipulation, R vector
## Use regexp to replace all the occurences of **** with NA
    lines2 <- gsub("\\*{3,5}", " NA", lines, perl=TRUE)
## Convert the character vector to a dataframe
    df <- read.table(
      textConnection(lines2), header=TRUE, colClasses = "character")
    closeAllConnections()
## Select monthly data in first 13 columns
    df <- df[,c(1,14)]
## Convert all variables (columns) to numeric format
    df <- colwise(as.numeric) (df)
    df[,2] <- df[,2]/100
    names(df) <- c("Year", "anom")
## Remove rows where Year=NA from the dataframe
    df <- df [!is.na(df$Year),]
## Find last report month and last value
    GISS_last <- nrow(df)
    GISS_last_yr <- df$Year[GISS_last]
    GISS_last_temp <- df$anom[GISS_last]
## Calc decade averages
  dec_mean<- as.numeric(14)
  dec_st <- as.numeric(14)
  dec_end <- as.numeric(14)
 # yr_n <- as.integer(df$Year)
  base_yr <- 1870
  df$dec_n <-  (as.numeric((df$Year - base_yr) %/% 10) * 10) + base_yr
 # df <- data.frame(df, dec_n)
  for (i in 1:13) {dec_st[i] = base_yr+ i*10
     dec_sub <- subset(df, dec_n == dec_st[i], na.rm=T)
     dec_mean[i] <- mean(dec_sub$anom)
     }
 dec_st[14] <- 2020              # Need to have for last step line across decade
 dec_mean[14] <- dec_mean[13]
 dec<- data.frame(dec_st, dec_mean)
#### Plot function
  plot_func <- function() {
  par(las=1); par(ps=12)
  par(oma=c(2.5,1,1,1)); par(mar=c(2.5,4,2,1))
# specify plot yr min & max
  p_xmin <- 1880;   p_xmax <- GISS_last_yr+10
  title <- paste("GISS Land and Sea Temperature Annual Anomaly Trend \n", p_xmin, " to ",
   GISS_last_yr, sep="")
  plot(df[,1], df[,2], type = "l", col = "grey",
     xlim = c(p_xmin, p_xmax), ylab = "Temperature Anomaly - C (1951-1980 Baseline)",
     xlab="", main = title,cex.main = 1,cex.lab=0.8,cex.axis=0.85)
  points(GISS_last_yr, GISS_last_temp, col = "red", pch=19)
  last_pt <- paste( GISS_last_yr, ", ", GISS_last_yr, " @ ", GISS_last_temp, "C",sep="")
  points(dec$dec_st, dec$dec_mean, type="s", col="blue")
## add legend
  legend(1880,0.6, c("Decade Mean Anomaly", "Annual Anomaly" ,GISS_last_temp), col = c("blue", "grey", "red"),
       text.col = "black", lty = c(1,1,0),pch=c(0,0,16),pt.cex=c(0,0,1),
       merge = T, bg = "white", bty="o", cex = .75, box.col="white")

Animated Images of Arctic Sea Ice Extent Decline

This post shows how to download and animate a series of Arctic Sea Ice Extent images using R and the animation package.

Continue reading

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

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")
}

Visualizing the Arctic Sea Ice Extent Decline

Understanding what is happening to Arctic sea ice is critical to recognizing the serious consequences of global warming. So I want to help people visualize the 30+ year trend in Arctic sea ice extent.

The source data file is here:  http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv

Comparison of GISS LOTAs During 5 El Nino – La Nina Cycles

In this post I compare GISS LOTAs during 5 El Nino – La Nina cycles (2010, 1998, 19883, 1973 and 1970).

El Nino – La Nina Cycles

In a previous post I showed the Nino 34 SSTA cycles for 2010, 1998, 1983, 1973 and 1970 here. In this post, I want to see how GISS Land Ocean Temperature Anomalies (LOTA) vary over El Nino – La Nina cycles.  Here is my RClimate chart showing GISS anomalies for 6 months prior to cycle year,the cycle year and the 12 months after cycle year (30 month period).

Click chart to enlarge

GISS_NINO34_cycles

While the 2010 cycle is only partially complete, there are a number of interesting aspects in this chart. The average temperatures during the cycles have clearly risen with the latest cycle showing the highest maximum anomaly. All 5 cycles all have similar patterns, with a buildup in 6 months prior to cycle year. The maximum – minimum range for the 5 cycles are comparable, ranging from 0.45 (2010) to 0.60 (1998).

Here is a data summary of the 5 cycles.