Tag Archives: R scripts

Arctic Amplification – November, 2013: Updated 12/20/13

Updated (bold italic) 12/20/13 to reflect comments from David. 

NASA’s GISS temperature anomaly map (link) for November, 2013 is reproduced  below. It uses a 2×2 degree grid cell for the globe. 

GISS_anom_map_11_13

The November, 2013 GISS temperature anomaly shows the critical global pattern that is important to recognize because it is fundamental to understanding why global warming is so dangerous.

First, the overall global anomaly for November, 2013 was 0.77 degrees C.  The 2 degree latitude zone mean anomaly varied from a low of 0.068 to a high of 2.33. So the mean global anomaly does not tell the full story, we need to look at the geographical distribution to really understand the global warming pattern.

As we examine the geographical distribution of the November, 2013 anomalies, we see that they tend to increase as we move from the equator toward the poles.  This pattern,  called polar amplification, means that the polar regions, particularly the Arctic region, warms much more rapidly than the overall global mean.

I developed this chart in R to display the mean zonal anomalies by 2 degree latitude zones to help me visualize the November, 2013 anomaly patterns.

art_amp

Here is the R Script that I used to produce the chart.
#### GISS Temperature Anomaly - Zonal mean by 2 degree latitude
##K O'Day, Dec. 18, 2013
##############################################################################################################################
 link <- c("http://data.giss.nasa.gov/tmp/gistemp/NMAPS/tmp_GHCN_GISS_ERSST_1200km_Anom11_2013_2013_1951_1980/nmaps_zonal.txt")
 mon <- "November, 2013"
 title <- paste("Mean Temperature Anomaly by 2 Degree Latitude Zones\n", mon, sep="" )
 note_1 <- "GISS Temperature Anomaly\n (1951-1980 base period)"
 df <- read.table(link, skip=4)
 par(las=1, oma=c(3,1,1,1), mar=c(5,5,3,1), ps=11)
 names(df)<- c("Zone", "Anom")
 #png(file="C://R_Home//Charts & Graphs Blog//RClimateTools//a_Revised_Blog//art_amp.png", bg="white")

 plot(df$Anom, df$Zone, xlim=c(0,3), ylim=c(-90,90), type="l", axes=F, xlab="Mean Anomaly for Zone - C",
      ylab = "Latitude",  xaxs="i", yaxs = "i", main=title)
   axis(1, at=NULL)
   axis(2, at=c(-90,-60,-30,0,30,60,90))
   abline(h=40, col="green")
   abline(h=0, col="darkgrey")
   abline(v=0.77,col = "black" )
   abline(h=64, col="blue")
   text(2.5, 43.5, "Philadelphia, Pa.", cex=0.7)
   text(2.7, 67, "Reykjavík, Iceland", cex=0.7)
   text(2.25, -20, note_1, cex=0.75, adj=0)
   rect(0.6,-65,0.85 , -50, col = "white", border = "white")
   text(0.77, -60, "Global Mean @ \n0.77 C", cex=0.7)
 mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,1, adj = 0, cex = 0.8, outer=T)
 mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.8, outer=T)
# dev.off()

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 - http://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")
}

Comparison of UAH and RSS Time Series with Common Baseline

In this post I set both UAH 5.4 and RSS 3.3 global temperature anomaly series to a common baseline period (1981-2010)  to compare them. Since both the UAH 5.4 and RSS 3.3 series are satellite based , they exhibit striking similarities.

Common Baseline

In this previous post, I showed how to convert temperature anomaly time series from one baseline period to another period.  I then used this technique in this post to directly compare UAH 5.4 (baseline 1981-2010) and GISS.

In this post, I compare the satellite based UAH 5.4 (baseline 1981-2010) and RSS 3.3 (baseline 1979-1998) series.

The offsets are as follows:

  • UAH:  -0.000978
  • RSS:      0.098772

Since the UAH TLT 5.4 series is based on a 1981-2010 baseline, the offset is nearly zero (-0.00098 versus 0.0). The RSS offset changes the baseline from 1979-1998 to 1981-2010.

Users can reproduce my analysis on their own by downloading my CTS.csv file and applying the offsets to the UAH and RSS series.

Comparison of 1981-2010 Baseline Series

Here is a plot of UAH and RSS 12 month moving averages for 1979 to current: Click to Enlarge

Continue reading

Comparison of UAH and GISS Time Series with Common Baseline

In this post I set both UAH and GISS global temperature anomaly series to a common baseline period (1981-2010)  and compare them. Even though the UAH series is satellite based and GISS series is station based, the series exhibit striking similarities.

Common Baseline

In this previous post, I showed how to convert temperature anomaly time series from one baseline period to another period.  I use this technique in this post to directly compare UAH (baseline 1981-2010) and GISS (baseline 1951-1980) series.

The offsets are as follows:

  • UAH:  -0.000978
  • GISS: 0.34958

Since the UAH TLT 5.4 series is based on a 1981-2010 baseline, the offset is nearly zero (-0.00098 versus 0.0).

Users can reproduce my analysis on their own by downloading my CTS.csv file and applying the offsets to the UAH and GISS series.

Comparison of 1981-2010 Baseline Series

Here is a plot of UAH and GISS 12 month moving averages for 1979 to current: Click to Enlarge Continue reading

UAH Temperature Anomalies Following Predictable Pattern

In this post I show one simple  and 2 multiple regression models to assess the role of time, El Nino – La Nina SSTA and volcanic activity (SATO) on UAH global temperature anomaly trends. The 3rd model provides a reasonable  approximation of the actual UAH oscillations over the 1979 – Feb, 2011 period.

Click Image to Enlarge

This analysis is similar to previous temperature anomaly regressions (here, here, here) that I have done.

The simple trend line regression shows the overall trend is upward, however, there are several oscillations that the linear trend misses.  The yr_frac and Nino34 regression improves on the linear model, however, it undershoots in the early 1980s,  overshoots in the 1992-1994 period, periods following significant volcanic activity.

The yr_frac, Nino34 and SATO model improves the fit in the early 1980s and 1992-1994 period and is slightly better in the 1998 and 2010 El Nino periods.

The 3rd model matches the observed 2010 El Nino – La Nina oscillation pretty well so far, indicating that the 2010 – 2011 UAH anomalies are following a predictable pattern.

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.

 

September 2011 Arctic Sea Ice Extent Forecast

In this post, I use a quadratic regression model to forecast the  September, 2011  Arctic Sea Ice Extent. The model was developed with  1980 – 2010 data. Links to the R script, source data and  how-to article on polynomial regression are provided.

Arctic Sea Ice Extent Forecast for September, 2011

First, here is my forecast: (Click image to enlarge)

ASIE_forecast_2011

Based on the 1980 – 2010 downward Arctic Sea Ice trend,  my forecast is that September, 2011 SIE will decline  0.36 below 2010 levels, to 4.54 million km^2, with a confidence band of +- 0.59.

How Did I Develop My Forecast?

I have written a number of posts on Arctic Sea Ice Extent (here, here, here). In this post, I used the NSDIC‘s monthly data file (link)  to construct a quadratic regression model of September sea ice extent for the 1980 – 2010 period. I then used this model to predict the September, 2011  Arctic Sea Ice Extent.

I have 2 main learning curve sources for this model:

  • Tamino‘s post on Arctic Sea Ice decline provided the basic idea of using a quadratic model to fit Arctic SIE decline.
  • John Quick’s tutorial on polynomial regression provided the how-to instructions I needed to implement Tamino’s approach in R.

RClimate Script and Links

Here is the link to my RClimate script.

Climate Time Series In a Single CSV File: Update 1

I am pleased to announce my CTS.csv file which includes 18 climate monthly time series in one easy to access csv file. This is part of  my goal of having a user friendly way for do-it-yourself citizen climate scientists to get up-to-date agency climate time series in a painless way.

Update 1: Reader Scott asked if I could provide meta data for the columns in my CTS.csv. This page lists the source agency and data links for the climate data series.

Here’s a snap shot of the first 6 rows of my  CTS.csv file. The data extends from 1880 until the most recent month.  Click image to enlarge

My hope is to make the CTS.csv the go-to file for citizen climate scientists who may want to:

  • Check temperature anomalies trends by series (GISS, HAD, NOAA, RSS, UAH)
  • Assess climate oscillations(AMO, AO, MEI, Nino34,  PDO)  trends
  • Evaluate  CO2 versus temperature anomaly relationships
  • Evaluate relationship between Sunspot numbers and anomaly temperature anomaly trends
  • Compare atmospheric transmission, SATO index  and volcanic activity
  • Assess impact of volcanoes on temperature anomaly trends
  • Compare MEI versus Nino ENSO 34 indicators
  • Assess lower stratospheric trends using RSS’s TLS series

By having these climate time series in a single csv file, R and Excel users can work with up to date data in a convenient form. The file will be automatically updated monthly as the climate agencies release their latest data.

How can CTS.csv Help Do-It-Yourself Citizen Climate Scientists?

Interested climate observers who want to compare global SSTA versus Nino34 trends, for example, have to follow a multiphase process:

  1. Find data file – even with Google this can take time
  2. Download files
  3. Merge 2 or more files to get data  into a usable format – source files all have different formats
  4. Perform analysis

Steps 1-3 can be very time consuming, so many users don’t bother checking out their ideas. Rather, they may rely on climate blog  comments. With CTS.csv and some R or Excel analysis, they can find the facts themselves rather than just having opinions.  They can submit their analysis and charts to blog posts, hopefully increasing the rigor of blog discussions.

Climate bloggers can request that their readers submit charts to back up their climate trend claims.

Data & RClimate Scripts Are All Open Book

All of the RClimate script that I use to produce the CTS.csv is available on-line at this link. Source data links are included in the function for each series.

Using RClimate To Retrieve Climate Series Data

This post shows how to use RClimate.txt to retrieve a climate time series and write a csv file in 5 lines of R script.

One of my readers, Robert, wants to be able to download climate time series data and write it to a csv file.  The R script below shows how to  download the MEI data series and write a csv file.  For this example I will use the RClimate function (func_MEI) to retrieve the data. I then simply specify the path and file name link for the output file (note quotes around the output file name and then write  a csv file.

source("http://processtrends.com/files/RClimate.txt"
m <- func_MEI()
head(m)
output_link <- "C://R_Home/mei.csv"
write.csv(m, output_link, quote=FALSE, row.names = F)

Continue reading