Heat Wave

The graphic shows the temperature anomaly in the US for the month of July (actually July 1 – 27). Parts of Kansas/Oklahoma experienced temperatures a full 4°C higher than normal over the month, while there was a cool anomaly in the Pacific Northwest.




The map uses an equal area projection (Albers) with horizontal and vertical scales in meters.

R aficionados will recognize use of the ggplot2 R package. We are experimenting with ggplot2 for weather map-making at Biospherica™ .


Data from CFSv2 Analysis.

ggplot2 : Elegant Graphics for Data Analysis, Hadley Wickham


Daily 2m temperature maps are available from operational CFSv2 here (tmp2m.l.gdas.201107.grib2). These correspond to initial conditions used by the CFSv2. To find a mean temperature for July, an average was taken over  4xdaily 6h forecast records from this grib2 file. There are probably many good ways to do this.  One way is to first interpolate from gaussian to latlon grid, then use the wgrib2 -netcdf to produce files which can be read directly by R’s raster package. To get a temperature anomaly the monthly 2m temperature climatology needs to be subtracted. This is available here (flxf06.cfsr.fclm.m07.1982.2010.grb2).

The temperature anomaly raster map can be projected from latlon to another PROJ4 map projection using raster::projectRaster. The resulting July global temperature anomaly raster object is called tmp2m.ras. USA country maps (level 0 & 1) are from gadm. These are reprojected and called USA.poly & USA.poly1. Finally, with tmp2m.ras, USA.poly and USA.poly1 objects in R memory,  ggplot2 was used to produce the map as follows:


#crop tmp2m.ras
x.min <- bbox(USA.poly)[1,1]
x.max <- bbox(USA.poly)[1,2]
y.min <- bbox(USA.poly)[2,1]
y.max <- bbox(USA.poly)[2,2]
tmp2m.ras <- crop(tmp2m.ras, extent(x.min,x.max,y.min,y.max))
#mask tmp2m.grd using country outline
USA.ras <- rasterize(USA.poly, tmp2m.grd)
tmp2m.ras <- mask(tmp2m.ras,USA.ras)
#convert data to dataframe for ggplot
tmp2m.ps <- rasterToPoints(tmp2m.ras)
df <- data.frame(tmp2m.ps)
colnames(df) <- c("easting","northing","tmp2m")
USA.poly <- fortify(USA.poly,region="ISO")
USA.poly1 <- fortify(USA.poly1,region="ID_1")
#now the pretty stuff
g <- ggplot(df,aes(x=easting,y=northing)) + scale_x_continuous(limits=c(x.min,x.max)) + scale_y_continuous(limits=c(y.min,y.max))
g <- g+geom_tile(aes(fill=tmp2m),alpha=0.6)+scale_fill_gradientn(name=expression(paste(degree,"C",sep="")),colours=rev(brewer.pal(9,"RdYlBu")))
g <- g+coord_equal()
g <- g+stat_contour(aes(z=tmp2m,colour=..level..),size=0.5, bins=4)+scale_colour_gradient(name=expression(paste(degree,"C",sep="")))
#country outline
g <- g + geom_path(data=US.poly, aes(x=long,y=lat,group=group), colour="grey40", alpha=1.0,size=0.5)
#state outline
g <- g + geom_path(data=USA.poly1, aes(x=long,y=lat,group=group), colour="grey80", alpha=1.0,size=0.3)
g <- g + opts(title = "USA Temperature July 2011")



  1. that map looks great! Can you make the code available? I’d love to see how you used both the CFS data and ggplot2…

  2. Hi Joe,
    Thank you for sharing this illustrative example.
    I am trying to make similar analysis over climate data.
    I am wondering about how to show easting and northing as Latitude and longitude values keeping the projection as it is.
    In other words, in the display map i want to show latitude and longitude values and not eastings and northings. Can you suggest me how can i do that ?
    Thank you,

  3. Hi Joe,
    I need help with something similar but i’ve been trying to use 2d- Density plots (ggplot2)…

    I’m trying to plot ecological distribution of some species of organisms i’m studying over the Arabian/Persian Gulf. Here is a sample of code i’ve tried,

    Backround layer

    nc <- get_map("Persian Gulf", zoom = 6, maptype = 'terrain', language = "English")
    ncmap <- ggmap(nc, extent = "device")

    i have the shape file for Arabian Gulf as well.

    Other layers

    stat_density2d(data=sample.data3, aes(x=long, y=lat, fill=..level.., alpha=..level..),geom="polygon")+
    geom_point(data=sample.data3, aes(x=long, y=lat))+
    geom_point(aes(x =50.626444, y = 26.044472), color="red", size = 4)+
    scale_fill_gradient(low = "green", high = "red") + scale_alpha(range = c(0.00, 0.25), guide = FALSE)

    but , i will like to use the stat_density2d to show the distributions of hundreds of species (which are recorded in columns e.g SP1….SPn) over the water body rather than just displaying latitude and longitude.

    Also, is it possible to restrict my heat map to just the water body ? I'll appreciate any help and recommendations i can get on this please

    Here is a link to stackoverflow where I initially asked the question with an image of what the problem look like.




Leave a Comment

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.