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 (such as the buffered points above, a transect, or within a polygon). For example, you might be interested in the mean elevation within the entire polygon in the above figure.
In this case study we’ll work with a HadCRUT temperature data from the Climatic Research Unit at the University of East Anglia. These are near-global rasters of surface temperature on a five degree grid.
Identify the hottest country on each continent by intersecting a set of polygons with a raster image and calculating the maximum raster value in each polygon.
world
datasetDownload starter R script (if desired). Save this directly to your course folder (repository) so you don’t lose track of it!
The details below describe one possible approach.
You will need to load the following packages
library(terra)
library(spData)
library(tidyverse)
library(sf)
Loading the spData()
package may return a warning:
To access larger datasets...install spDataLarge...
. This is
not required - you can use the standard lower resolution files and
safely ignore this message.
Download the mean annual temperatures over the reference period 1961-1990 Climatic Research Unit data (CRU) here. Absolute temperatures for the base period 1961-90 on a 5° by 5° grid. Download these data in netcdf format using the code below:
library(ncdf4)
download.file("https://crudata.uea.ac.uk/cru/data/temperature/absolute.nc","crudata.nc")
# read in the data using the rast() function from the terra package
tmean=rast("crudata.nc")
Note: If the above code returns an error about
nc_open()
, try adding method="curl"
at the end
of the download.file()
command.
tmean=rast("crudata.nc")
).tmean
object (you can start by just
typing it’s name tmean
, then perhaps making a
plot()
). How many layers does it have? What do these
represent? You can read more about the data heremax()
to calculate the maximum value observed in
each pixel across all months. You may want to plot()
the
output of this and compare the plot of the original (full) raster. How
many layers does it have now?terra::extract()
to identify the maximum
temperature observed in each country (fun=max
). Also set
na.rm=T, small=T
to 1) handle missing data along coastlines
and 2) account for small countries that may not have a pixel centroid in
them.bind_cols()
to bind the original world
dataset with the new summary of the temperature data to create a new
object called world_clim
with the outputs from the
extract()
function above.ggplot()
and geom_sf()
to plot the
maximum temperature in each country polygon (not the original raster
layer). To recreate the image below, you will also need
+scale_fill_viridis_c(name="Maximum\nTemperature (C)")
.
Note the use of \n
to insert a line break in the text
string. You can move the legend around with
+theme(legend.position = 'bottom')
.dplyr
tools to find the hottest country in each
continent. You may need group_by()
and
slice_max
. To create a nice looking table, you may also
want to use select()
to keep only the desired columns,
arrange()
to sort them, st_set_geometry(NULL)
to drop the geometry column (if desired). Save this table as
hottest_continents
.Your final result should look something like this:
And the summary table will look like this:
name_long | continent | max |
---|---|---|
Mali | Africa | 35.7 |
Algeria | Africa | 35.7 |
Kuwait | Asia | 34.8 |
Saudi Arabia | Asia | 34.8 |
Australia | Oceania | 32.0 |
United States | North America | 29.5 |
Brazil | South America | 28.2 |
Albania | Europe | 25.0 |
French Southern and Antarctic Lands | Seven seas (open ocean) | 7.5 |
Antarctica | Antarctica | -2.0 |
Note that these data are based on 0.5 degree resolution data and thus will miss small hot places and cannot be directly compared with station-based data.
Build a leaflet map of the same dataset.