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.
# Session a4: Species range maps
and 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.
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.
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.080348
We 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
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_araucana
Currently, the BIEN
package 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
SpatVector
object 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.
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] 5783633
We 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
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).
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 : 1
map('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)
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: sp
library(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.94637
values(monkey_puzzle_terra)
## species
## 1 Araucaria_araucana
monkey_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
terra
package and the
letsR
package.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,] 0
Exercise:
terra
or
letsR
) and map species richness of the species within that
genus.