Species range maps
RStudio project
Open the RStudio project that we created in the first session. I recommend to use this RStudio project for the entire course and within the RStudio project create separate R scripts for each session.
- Create a new empty R script by going to the tab “File”, select “New File” and then “R script”
- In the new R script, type
# Session a4: Species range mapsand save the file in your folder “scripts” within your project folder, e.g. as “a4_RangeMaps.R”
We are living in an age of big data where biodiversity are becoming increasingly available in digital format and at a global scale (Wüest et al. 2020). Many different types of biodiversity data exist, e.g. from standardised monitoring schemes, citizen science platforms, or expert knowledge. Each of these data types comes with own challenges. Within this module, I want to provide some first impressions on how you can obtain and process typical types of species data. Specifically, in this session, we will work with range maps of terrestrial animals and plants, while occurrence data for terrestrial species in GBIF (the Global Biodiversity Information facility) and for marine species in OBIS (the Ocean Biodiversity Information System) will be covered later in the species distribution modelling part.
1 Obtaining range maps
We rarely have detailed biodiversity data available over large geographic extents. At broad (continental to global) extents, expert-drawn range maps (also called extent-of-occurrence maps) are often the primary data source on species distributions.
1.1 IUCN range maps
The IUCN (the International Union for the Conservation of Nature) provides expert range maps for a large number of species including mammals, birds (through BirdLife International), amphibians, reptiles, and freshwater fishes. There are also some range maps on plants and marine species, but these are very limited taxonomically. Have a look for which taxa range maps are available: https://www.iucnredlist.org/resources/spatial-data-download. You can download them for free, but you should provide some information on your work to obtain the data.
Most of the IUCN data are provided in the form of shapefiles. In
practical 2, we have already loaded the range map of the Alpine
shrew and learned how to use the package terra for
reading in the shapefiles. You may remember that the shapefile is
recognized as SpatVector object.
library(terra)
# Load the shapefile
(shrew <- terra::vect('data/IUCN_Sorex_alpinus.shp'))## class : SpatVector
## geometry : polygons
## dimensions : 1, 27 (geometries, attributes)
## extent : 5.733728, 26.67935, 42.20601, 51.89984 (xmin, xmax, ymin, ymax)
## source : IUCN_Sorex_alpinus.shp
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : id_no binomial presence origin seasonal compiler yrcompiled
## type : <chr> <chr> <int> <int> <int> <chr> <int>
## values : 29660 Sorex alpinus 1 1 1 IUCN 2008
## citation source dist_comm (and 17 more)
## <chr> <chr> <chr>
## IUCN (Internat~ NA NA# Plot the Central Europe
library(maps)
map('world',xlim=c(5,30), ylim=c(40,55))
# Overlay the range of the Alpine Shrew
plot(shrew, col='red', add=T)
Unfortunately, there is no API for the IUCN range maps. So, you need to register with IUCN and then you can download range maps for different taxonomic groups. The range maps for birds come as separate shape files per species, while the range maps of mammals are provided as single shapefile that contains all species. Course participants can download the mammal shapefile from the moodle course - of course, the license agreement by the IUCN apply!
# Read shapefile for all mammals using the raster package:
mammals <- terra::vect('data/MAMMTERR.shp')The shapefile contains the range polygons for all described mammal species. The attribute table contains information on species’ PRESENCE, ORIGIN, and SEASONALity. Please have a look at the metadata to understand the different values these attributes can be coded as.
mammals## class : SpatVector
## geometry : polygons
## dimensions : 43444, 14 (geometries, attributes)
## extent : -180.0001, 180.0025, -56.66595, 90 (xmin, xmax, ymin, ymax)
## source : MAMMTERR.shp
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : BINOMIAL PRESENCE ORIGIN SEASONAL COMPILER YEAR
## type : <chr> <int> <int> <int> <chr> <int>
## values : Amblysomus cor~ 1 1 1 IUCN 2008
## Amblysomus hot~ 1 1 1 IUCN 2008
## Amblysomus mar~ 1 1 1 IUCN 2008
## CITATION SOURCE DIST_COMM ISLAND (and 4 more)
## <chr> <chr> <chr> <chr>
## IUCN (Internat~ NA NA NA
## IUCN (Internat~ NA NA NA
## IUCN (Internat~ NA NA NA# Inspect attribute table
terra::head(mammals)## BINOMIAL PRESENCE ORIGIN SEASONAL COMPILER YEAR
## 1 Amblysomus corriae 1 1 1 IUCN 2008
## 2 Amblysomus hottentotus 1 1 1 IUCN 2008
## 3 Amblysomus marleyi 1 1 1 IUCN 2008
## 4 Amblysomus robustus 1 1 1 IUCN 2008
## 5 Amblysomus septentrionalis 1 1 1 IUCN 2008
## 6 Amblysomus septentrionalis 1 1 1 IUCN 2008
## CITATION SOURCE DIST_COMM ISLAND
## 1 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## 2 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## 3 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## 4 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## 5 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## 6 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## SUBSPECIES SUBPOP SHAPE_area SHAPE_len
## 1 <NA> <NA> 4.93503476 15.575123
## 2 <NA> <NA> 19.94471634 28.392110
## 3 <NA> <NA> 0.12779360 1.902574
## 4 <NA> <NA> 0.09936845 1.655898
## 5 <NA> <NA> 0.25763915 2.078154
## 6 <NA> <NA> 2.70333866 7.080348We can search for specific species or species groups in the attribute table in different ways:
# Range map for the species 'Lynx lynx'
terra::subset(mammals, mammals$BINOMIAL == "Lynx lynx")
# Show all entries for the species with the word 'Lynx' in their name
grep('Lynx',mammals$BINOMIAL, value=T)
# Range map for all species with the word 'Lynx' in their name
mammals[grep('Lynx', mammals$BINOMIAL),]We can use the attribute table subsets to select specific polygons that we want to look at.
# Assign range map of the Eurasian lynx to separate object
lynx_lynx <- mammals[mammals$BINOMIAL=='Lynx lynx',]
# Map the range
map('world')
plot(lynx_lynx, col='red', add=T)
Test it yourself
- Select another mammal species and plot the range map.
1.2 BIEN range maps
The BIEN
database (Botanical Information and Ecology Network) contains many range
maps on plants, but unfortunately only for the Americas. These can be
accessed using the BIEN package. As illustrative example,
we load the range map for the monkey puzzle tree (or Chilean pine -
Araucaria araucana).
library(BIEN)
# Load the range map for the monkey puzzle
(monkey_puzzle_raster <- BIEN_ranges_load_species('Araucaria_araucana'))## class : SpatialPolygonsDataFrame
## features : 1
## extent : -109.9819, -55.46083, -55.73974, 62.67197 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 1
## names : species
## value : Araucaria_araucanaCurrently, the BIENpackage still relies on the spatial R
package raster that was the predeccesor of
terra and is still being used by many packages. Here, we
simply convert the SpatialPolygonsDataFrame object into the
SpatVectorobject that we already know from the
terra package.
# Change into SpatVector object (terra package)
(monkey_puzzle_terra <- terra::vect(monkey_puzzle_raster, crs='+proj=longlat +datum=WGS84'))## class : SpatVector
## geometry : polygons
## dimensions : 1, 1 (geometries, attributes)
## extent : -109.9819, -55.46083, -55.73974, 62.67197 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : species
## type : <chr>
## values : Araucaria_araucana# Map
map('world',xlim = c(-180,-20),ylim = c(-80,80))
plot(monkey_puzzle_terra,col='red',add=T)
The native range of the mokey puzzle is in the Chilean Andes, so clearly the range maps also show areas where the species naturalized.
2 Working with range maps
2.1 Range size and range centre
The terra package allows us to easily calculate the area
of the polygons, meaning the range size of our species. The function
expanse() outputs the area in square meters, kilometers or
hectars.
# Range area of alpine shrew in square kilometers:
terra::expanse(shrew, unit="km")## [1] 490543.3# Range area of monkey puzzle in square kilometers:
terra::expanse(monkey_puzzle_terra, unit="km")## [1] 5783633We can also very easily calculate the centre of gravity or range
centroid from the SpatVector object.
# Range centroid:
terra::centroids(shrew)## class : SpatVector
## geometry : points
## dimensions : 1, 27 (geometries, attributes)
## extent : 16.02645, 16.02645, 46.90336, 46.90336 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : id_no binomial presence origin seasonal compiler yrcompiled
## type : <chr> <chr> <int> <int> <int> <chr> <int>
## values : 29660 Sorex alpinus 1 1 1 IUCN 2008
## citation source dist_comm (and 17 more)
## <chr> <chr> <chr>
## IUCN (Internat~ NA NA# Map the species range and add the centroid to the map
map('world',xlim=c(5,30), ylim=c(40,55))
plot(shrew, col='red', add=T)
points(terra::centroids(shrew), col='blue',cex=3,pch="+")
We need to be careful how to interpret these centroids. They represent the centre of gravity, i.e. the mean coordinates of the distribution (weighted by cell size) but obviously, if we have several patches, the centroid might not even fall within an occupied patch.
In case of the Monkey Puzzle, we have suspected that the range maps also contain non-native areas. We can clip the range maps to a desired spatial extent and then only calculate the range centroid for this (presumed) native range.
# Let's crop the range polygons to (rough coordinates of) South America
monkey_puzzle_SAm <- terra::intersect(monkey_puzzle_terra, ext(-85, -30, -55, 5))
# Map the range and range centroid
map('world',xlim = c(-100,-10),ylim = c(-60,15))
plot(monkey_puzzle_SAm,col='red',add=T)
points(terra::centroids(monkey_puzzle_SAm), col='blue',pch='+', cex=3)
Test it yourself
- Use the range map of the mammal species you selected earlier, and calculate the range size and range centroid. Add the point location of your centroid to the range map that you have plotted previously.
2.2 Rasterising range maps
For many applications in macroecology, we need to rasterise the polygons. The problem is that it is unclear at which spatial resolution the range maps accurately represent species occurrences. Hurlbert and Jetz (2007) and Jetz, McPherson, and Guralnick (2012) define the minimum spatial resolution as 100-200km (1-2°), although also resolutions of 50km (0.5°) and finer have been used (Krosby et al. 2015; Zurell et al. 2018).
2.2.1 Rasterising range
maps with terra
Rasterising polygon data is made very easy in the terra
package. We first have to define a SpatRaster object of the
desired resolution, and then transfer the polygon data to the raster
cells.
# By default, terra() will create a 1° resolution map in the *WGS 84* coordinate system (lon/lat).
(r_1deg <- terra::rast())## class : SpatRaster
## dimensions : 180, 360, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84# Now, rasterise the shrew polygon data to the 1° raster grid
(shrew_1deg <- terra::rasterize(shrew, r_1deg))## class : SpatRaster
## dimensions : 180, 360, 1 (nrow, ncol, nlyr)
## resolution : 1, 1 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source : memory
## name : layer
## min value : 1
## max value : 1map('world',xlim=c(5,30), ylim=c(40,55))
plot(shrew, col='red', add=T)
plot(shrew_1deg, add=T, alpha=0.6, legend=F)
Obviously, the margins of the range polgyon and the raster map differ at several places.
Test it yourself
Check out the help page ?rasterize and find out what the
argument touches is doing. Rasterise the shrew range map
again with setting a different touches argument and map the
result. What is the difference?
We look at a second example, the lynx:
# Define an empty SpatRaster of the world at 2° spatial resolution
(r_2deg <- terra::rast(res=2))## class : SpatRaster
## dimensions : 90, 180, 1 (nrow, ncol, nlyr)
## resolution : 2, 2 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84# Rasterize the eurasian lynx data
lynx_lynx_2deg <- terra::rasterize(lynx_lynx, r_2deg)
# Map the occupied grid cells
map('world')
plot(lynx_lynx_2deg, add = TRUE, legend = FALSE)
2.2.2 Rasterising range
maps with letsR
There are also specific macroecological packages in R that facilitate
working with range maps and rasterising them, for example the function
lets.presab() in the letsR package. At the
moment, this package still relies on the older spatial R package
raster.
library(letsR)## Lade nötiges Paket: raster## Lade nötiges Paket: splibrary(raster)
# The lets.presab() function expects "binomial" as one specific column name in the Polygons data frame
values(shrew)## id_no binomial presence origin seasonal compiler yrcompiled
## 1 29660 Sorex alpinus 1 1 1 IUCN 2008
## citation source dist_comm island
## 1 IUCN (International Union for Conservation of Nature) <NA> <NA> <NA>
## subspecies subpop legend tax_comm kingdom phylum class
## 1 <NA> <NA> Extant (resident) <NA> ANIMALIA CHORDATA MAMMALIA
## order_ family genus category marine terrestial freshwater SHAPE_Leng
## 1 EULIPOTYPHLA SORICIDAE Sorex NT False True False 88.67521
## SHAPE_Area
## 1 57.94637values(monkey_puzzle_terra)## species
## 1 Araucaria_araucanamonkey_puzzle_raster@data## species
## 1 Araucaria_araucana# As the BIEN data do not contain "binomial" as column name, we have to add it:
colnames(monkey_puzzle_raster@data) <- "binomial"
# We set the resolution to 1 degree (the default) and restrict the spatial extent to South America
r_monkey_puzzle <- lets.presab(monkey_puzzle_raster, resol=1, xmn = -100, xmx = -10, ymn = -57, ymx = 15)
# This time, we receive a raster object (instead of SpatRaster)
r_monkey_puzzle##
## Class: PresenceAbsence
## Number of species: 1
## Number of cells: 243
## Resolution: 1, 1 (x, y)# Map the range and range centroid
map('world',xlim = c(-100,-10),ylim = c(-60,15))
plot(monkey_puzzle_SAm,col='blue',add=T)
plot(r_monkey_puzzle, add=T, alpha=0.6, legend=F)
Test it yourself
- Rasterise the range map of the mammal species you selected above.
Try both approaches using the
terrapackage and theletsRpackage.
2.2.3 Bulk-rasterising
multiple species range maps with letsR
The letsR package also allows to bulk-download multiple
species and rasterise them to form a richness map. As first example, we
look at the Pinus genus in the Americas
# Extract the available Pinus species names
(pinus_names <- BIEN_ranges_genus("Pinus",match_names_only = T)[,1])## [1] "Pinus_albicaulis" "Pinus_aristata" "Pinus_arizonica"
## [4] "Pinus_armandii" "Pinus_attenuata" "Pinus_ayacahuite"
## [7] "Pinus_balfouriana" "Pinus_banksiana" "Pinus_brutia"
## [10] "Pinus_bungeana" "Pinus_canariensis" "Pinus_caribaea"
## [13] "Pinus_cembroides" "Pinus_clausa" "Pinus_contorta"
## [16] "Pinus_coulteri" "Pinus_culminicola" "Pinus_densiflora"
## [19] "Pinus_devoniana" "Pinus_douglasiana" "Pinus_durangensis"
## [22] "Pinus_echinata" "Pinus_edulis" "Pinus_elliottii"
## [25] "Pinus_engelmannii" "Pinus_flexilis" "Pinus_georginae"
## [28] "Pinus_glabra" "Pinus_greggii" "Pinus_halepensis"
## [31] "Pinus_hartwegii" "Pinus_herrerae" "Pinus_jeffreyi"
## [34] "Pinus_koraiensis" "Pinus_lambertiana" "Pinus_lawsonii"
## [37] "Pinus_leiophylla" "Pinus_longaeva" "Pinus_luchuensis"
## [40] "Pinus_lumholtzii" "Pinus_luzmariae" "Pinus_maximartinezii"
## [43] "Pinus_maximinoi" "Pinus_monophylla" "Pinus_montezumae"
## [46] "Pinus_monticola" "Pinus_mugo" "Pinus_muricata"
## [49] "Pinus_nelsonii" "Pinus_nigra" "Pinus_oocarpa"
## [52] "Pinus_palustris" "Pinus_parviflora" "Pinus_patula"
## [55] "Pinus_peuce" "Pinus_pinaster" "Pinus_pinceana"
## [58] "Pinus_pinea" "Pinus_ponderosa" "Pinus_praetermissa"
## [61] "Pinus_pringlei" "Pinus_pseudostrobus" "Pinus_pumila"
## [64] "Pinus_pungens" "Pinus_quadrifolia" "Pinus_radiata"
## [67] "Pinus_remota" "Pinus_resinosa" "Pinus_rigida"
## [70] "Pinus_roxburghii" "Pinus_rzedowskii" "Pinus_sabiniana"
## [73] "Pinus_serotina" "Pinus_strobiformis" "Pinus_strobus"
## [76] "Pinus_sylvestris" "Pinus_tabuliformis" "Pinus_taeda"
## [79] "Pinus_tecunumanii" "Pinus_teocote" "Pinus_thunbergii"
## [82] "Pinus_torreyana" "Pinus_virginiana" "Pinus_wallichiana"# Download the range maps for all Pinus species (as SpatialPolygonsDataFrame using the raster standard instead of terra)
(pinus <- BIEN_ranges_load_species(pinus_names))## class : SpatialPolygonsDataFrame
## features : 84
## extent : -179.1256, 179.8453, -55.95014, 63.61258 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 1
## names : species
## min values : Pinus_albicaulis
## max values : Pinus_wallichiana# As the BIEN data do not contain "binomial" as column name, we have to re-format the column names again before rasterising
colnames(pinus@data) <- "binomial"
r_pinus <- lets.presab(pinus, resol=1, xmn = -170, xmx = -10, ymn = -57, ymx = 60)
# Plot species richness
plot(r_pinus)
In the lets.presab() function, we can also specify which
PRESENCE, ORIGIN and SEASONAL information should be used and which not
(see IUCN metadata on moodle). Let’s look at Neotropical fruit bats
(Artibeus) as example. Here, we set
presence=1 meaning that we only consider extant species,
origin=1 meaning only native species, and
seasonal=1 meaning only resident species.
As the letsR package does not work with
SpatRaster objects yet, we have to use the older
raster package to prepare the range map data.
library(raster)
# Subset the SpatVector
artibeus_spp <- mammals[grep('Artibeus',mammals$BINOMIAL),]
# Convert to SpatialPolygonsDataFrame from raster package
artibeus_spp <- as(artibeus_spp, "Spatial")
# Rasterize the ranges using the letsR package
r_artibeus_spp <- lets.presab(artibeus_spp, resol=2,
presence = 1, origin = 1, seasonal = 1)
# Map the species richness
plot(r_artibeus_spp)
# Map single species - here, just the first two
par(mfrow=c(1,2))
plot(r_artibeus_spp, name = "Artibeus amplus")
plot(r_artibeus_spp, name = "Artibeus anderseni")
# Look at structure of the object and at the presence-absence matrix
str(r_artibeus_spp, 1)## List of 3
## $ Presence_and_Absence_Matrix: num [1:339, 1:22] -109 -109 -107 -105 -99 -83 -81 -105 -103 -99 ...
## ..- attr(*, "dimnames")=List of 2
## $ Richness_Raster :Formal class 'RasterLayer' [package "raster"] with 12 slots
## $ Species_name : chr [1:20] "Artibeus amplus" "Artibeus anderseni" "Artibeus aztecus" "Artibeus cinereus" ...
## - attr(*, "class")= chr "PresenceAbsence"head(r_artibeus_spp$Presence_and_Absence_Matrix)## Longitude(x) Latitude(y) Artibeus amplus Artibeus anderseni
## [1,] -109 27 0 0
## [2,] -109 25 0 0
## [3,] -107 25 0 0
## [4,] -105 23 0 0
## [5,] -99 23 0 0
## [6,] -83 23 0 0
## Artibeus aztecus Artibeus cinereus Artibeus concolor Artibeus fimbriatus
## [1,] 0 0 0 0
## [2,] 0 0 0 0
## [3,] 0 0 0 0
## [4,] 1 0 0 0
## [5,] 1 0 0 0
## [6,] 0 0 0 0
## Artibeus fraterculus Artibeus glaucus Artibeus gnomus Artibeus hirsutus
## [1,] 0 0 0 1
## [2,] 0 0 0 1
## [3,] 0 0 0 1
## [4,] 0 0 0 0
## [5,] 0 0 0 0
## [6,] 0 0 0 0
## Artibeus incomitatus Artibeus inopinatus Artibeus jamaicensis
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 1
## [5,] 0 0 1
## [6,] 0 0 1
## Artibeus lituratus Artibeus obscurus Artibeus phaeotis
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 1 0 0
## [4,] 1 0 1
## [5,] 0 0 0
## [6,] 0 0 0
## Artibeus planirostris Artibeus rosenbergii Artibeus toltecus
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 1
## [4,] 0 0 1
## [5,] 0 0 1
## [6,] 0 0 0
## Artibeus watsoni
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0
## [6,] 0Exercise:
- Select another genus of plants or mammals, and follow the workflow
to rasterise range maps (either using
terraorletsR) and map species richness of the species within that genus. - Plot the latitudinal species richness gradient for this genus. Interpret in light of previously discussed patterns.