Background

The raster data format is commonly used for environmental datasets such as elevation, climate, soil, and land cover. We commonly need to extract the data from raster objects using simple features (vector objects). For example if you had a set of points you collected using a GPS and wanted to know the mean annual temperature at each point, you might extract the data from each location in a raster-based map of temperature.

You could also be interested in some summary of the raster data across multiple pixels. For example, you might be interested in the mean elevation in a polygon.

In this case study I work with a timeseries of temperature data from WorldClim (http://worldclim.org). These are near-global rasters of various climatic variables available at several resolutions. For convenience, I work with the very coarse data (10 minutes, which is about 10km), but much finer data are available (~1km).

Objective

Identify the hottest country on each continent (not counting Antarctica) by intersecting a set of polygons with a raster image and calculating the maximum raster value in each polygon.

Method

Libraries

The following packages must be loaded

library(raster)
library(sf)
library(sp)
library(spData)
library(tidyverse)

knitr::opts_chunk$set(cache=TRUE) # cache the results for quick compiling

Data

Retrieve ‘world’ data for country polygons from the spData package

data(world)

Download the WorldClim maximum temperature dataset at the lowest resolution (10 degrees)

tmax_monthly = getData(name = "worldclim", var = "tmax", res = 10)

Prepare country polygon data

Remove “Antarctica” with ‘filter()’ because WorldClim does not have data there

world = filter(world, continent != "Antarctica", region_un != "Antarctica")

Convert ‘world’ object to sp format because raster does not accept sf objects

as(world,"Spatial")

Prepare climate data

Inspect tmax_monthly

tmax_monthly
## class      : RasterStack 
## dimensions : 900, 2160, 1944000, 12  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 
## names      : tmax1, tmax2, tmax3, tmax4, tmax5, tmax6, tmax7, tmax8, tmax9, tmax10, tmax11, tmax12
plot(tmax_monthly)

Convert to degrees C; WorldClim temperatures are recorded in 10*degrees C

tmax_monthly = tmax_monthly*0.1

Calculations

Calculate the annual maximum temperature around the world

Create ‘tmax_annual’ as the annual max temperature of each pixel in the raster stack

tmax_annual = max(tmax_monthly)

Rename the layer in ‘tmax_annual’ to ‘tmax’

names(tmax_annual) = "tmax"

Calculate the maximum temperature observed in each country

Identify the maximum temperature observed in each country, handling missing data along coastlines and accounting for small countries with less than 0.5 degree pixel

world_tmax = raster::extract(tmax_annual,world,small=T,fun=max,na.rm=T,sp=T)

Convert ‘world_tmax’ to sf format

world_tmax = st_as_sf(world_tmax)

‘tmax’ has now been added as a column to ‘world_tmax’ containing country polygons

Result

Plot the maximum temperature in each country polygon

Find hottest country in each continent

table = select(world_tmax, name_long, continent, tmax) %>% 
  st_set_geometry(NULL) %>% 
  group_by(continent) %>% 
  slice_max(tmax,n=1)
Country Continent Annual Maximum Temperature (C)
Algeria Africa 48.9
Iran Asia 46.7
Spain Europe 36.1
United States North America 44.8
Australia Oceania 41.8
French Southern and Antarctic Lands Seven seas (open ocean) 11.8
Argentina South America 36.5

Conclusion

The hottest country in each continent has been identified in the table above.