Mapping Temperature Anomalies with R: UPDATE 12/17/13

In this post I show how to map NASA GISS’s  2×2 degree temperature anomaly data using R mapping tools. Rather than rely on a single value to reflect monthly global temperature anomaly, this map shows the anomalies  in each of the 16,200 cells in a  2 degree lon/lat grid. This lets us see the details that make up the global mean, we can see which areas are warmer and which are cooler.  I provide a link to my RClimate script and data file so that interested R users can make their own maps.

Here’s my R Climate map of NASA’s July 2010 2×2 degree data set. (Click map to check out the enlargement )

Users can use NASA’s interactive tool to make a comparable map at this link.  Here’s what NASA’s map version looks like.

NASA’s map is very good, so why bother making my version in R? I had several reasons:

  • I am not happy with NASA’s color scheme, I wanted to see if I could improve it
  • Northern Hemisphere polar amplification is an important issue that I think gets masked with NASA’s color scheme because so a large area shows up “reddish”, Does better color separation help to show NH polar amplification?
  • I don’t like NASA’s color bar scale, notice how the intervals vary by location along the scale axis.  Interval span range from a low of 0.1 to a high of 1.8 .  I’m sure there’s a good reason, but I wanted to see how other scales look.
  • Size, I wanted to make a larger image so that I could magnify the zoom farther to see specific areas in more detail.
  • I wanted to learn how to map climate data myself.

RClimate Map Features

Here are the features I was trying to incorporate into my map:

  • Larger size (1200 * 1000) – feel free to zoom in
  • Increased color separation to distinguish  really cool/warm areas
  • Uniform 1 degree C scale bar intervals
  • Missing data shown in white
  • Can show entire globe or map only a portion of the globe by adjusting the “bounding box”

RClimate Script Details

Here’s the R Script:

###################################################################################################
##                         NASA_global_2x2_grid_lota_map/R                                       ##
## Read NASA GISS global temp anom file for month, anom data in 2 degree grid cells              ##
## D Kelly O'Day       - http://chartsgraphs.wordpress.com                                       ##
## Uses R's powerfull geospatial analysis tools                                                  ##
## R script help provided by Andy South of UK's NERC                                             ##
## Orig script: 9/2/10; Updated 12/17/13 to reflect NASA's revised data file naming conventions  ##                                                                   ##
###################################################################################################
# Load map related libraries
  library(fields)    # needed for image.plot()
  library(sp)        # needed for coordinates, SpatialPointsDataFrame
  library(maptools)  # neeed for wrld_simpl map

# Get world shape for map background
  data(wrld_simpl)            ## from maptools package
  shp <- wrld_simpl

# Read source data file
 # link to NASA gistemps map page: http://data.giss.nasa.gov/gistemp/maps/
 # txt file of 2x2 degree anomalies available
 which_mo <- "09_2013"
 link_1 <- "http://data.giss.nasa.gov/tmp/gistemp/NMAPS/tmp_GHCN_GISS_ERSST_1200km_Anom"
 link_2 <- "_2013_1951_1980/nmaps.txt"
 link <- paste(link_1, which_mo, link_2, sep="")
 rdf <- read.table(link, skip = 1, sep = "", header=T)
 names(rdf) <- c("i", "j", "lon", "lat", "anom")

###############################
# Convert all anom data with 9999.0000 to NA
  rdf$anom[rdf$anom==9999.0000] <- NA         # convert all 9999.0000 to NA

## Promote to SpatialPointsDataFrame
  points_df <- rdf          # make copy of original file
  coordinates(points_df) = c("lon", "lat")     # convert to sp file with lon/lat cordinates
## Promote to SpatialPixelsDataFrame
  pixel_df <- points_df
  gridded(pixel_df) <- TRUE
## Promote to SpatialGridDataFrame
  rdf_sp = as(pixel_df, "SpatialGridDataFrame")
 main_title <- paste("GISS Temp Anomaly ", which_mo, sep="")
## image_func()  Plot Function
  image_func <- function() {
    par(mar=c(2.5,2,1,1)) ; par(oma=c(0,0,0,0))  # set plot par(mar=)
    par(pty="s") #square
    leg_title <- expression(paste("Anomaly - ",degree*C, " (Baseline: 1951-1980)", sep=""))
    g_plot <- as.image.SpatialGridDataFrame(rdf_sp["anom"])
    image.plot(g_plot, main=main_title, nlevel=12,las=1,
      axes=T, horizontal=T,  ylim=c(-90,90), xlim=c(-181,180), zlim = c(-6,6),
      legend.lab = leg_title, legend.mar=4,legend.shrink=0.5)
   plot(shp, border="black", add=T)
   grid(col = "lightgrey", lty=1)    # optional grid on/off?
  }

## Plot image on printer
  image_func()

About these ads

9 responses to “Mapping Temperature Anomalies with R: UPDATE 12/17/13

  1. When I click on Rclimate script link, I get this: The page cannot be found

  2. Sorry Kelly, where I can find the script? I am new in R

  3. I’d like to suggest that in addition to a white region around zero, that you pick lightly saturated tones for the next few bands out. Save the strongly saturated colors for the extreme values.

    Also, this map projection hugely overstates both the importance of and knowledge of the poles. Pick a projection that doesn’t have such a large area error. This is a hard problem, and there are no perfect answers.

    An underlying problem is the use of a 2×2 grid cell in the various data series. That is probably acceptable in the tropics (but even so the boxes are not equal in area). But it is just wrong for covering an entire sphere. There are much better meshes that could be used for surface data sets, where the individual mesh cells are constrained to equal (or nearly equal) area.

  4. One thing I think NASA does do right is the white area. This represents near normal temperatures. When you have small fluctuations from normal different colours for negative and positive values, large areas slightly above or below normal by an insignificant amount can give a very misleading visual impression.

    • swiftright

      I agree with your point about having a white area for near zero areas. I’m working on a revised color scheme for August, 2010 anomalies.

      link

  5. Hello,

    thanks, very nice!
    I found two little errors in your R script:
    * the url to the data is wrong
    * the separator in the read.table command should be ‘,’

    • Stefan

      Thank you for catching these and letting me know so that others won’t have the same problem you did.

      It should be working properly now. Please let me know if you still have any problems with my script.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s