From point occurrences to range maps
Biodiversity Informatics Workshop
Madagascar
dismo package (Hijmans & Elith 2011) and from our friend Vijay Barve (https://vijaybarve.wordpress.com)So, we’ve already seen where do the primary biodiversity data comes from (e.g. point occurrences from inventories, collections, museums, etc.). Some of this data is already available in the internet!!
Install and load the following packages
library(rgdal)
library(rgeos)
library(sp)
library(raster)
library(maptools)
library(maps)
library(rgbif)
library(dismo)
library(ggmap)
library(rvertnet)Get some occurrence data for a particular species from GBIF, directly within R. This may take some time, given the number of occurrences for the selected species. Here, we’ll work with the fossa! NOTA: we need to have an internet connection
fossa <- gbif("Cryptoprocta","ferox")Take a look at the object dimensions (i.e. how many lines and columns)
dim(fossa)Check the column names, so we can look for the coordinates later
names(fossa)Create another object keeping only the coordinates (longitude/latitude)
fossa_points <- fossa[,c("lon","lat")] NOTE: the variable names’ may be different (.e.g “LATITUDE”, “Latidude”, “lat”, etc.). Hence, it is always good to check the column names (as above)
Plot the points (x,y)
plot(fossa_points[,1],fossa_points[,2],col="blue",pch=19)Add a world map, for reference. We need to load a dataobject from the maptools package
data(wrld_simpl)
plot(wrld_simpl,add=T)Are there any “weird” points?
Do the opposite, first plot the worldmap and then the points
plot(wrld_simpl,col="light yellow")Add the points
points(fossa_points[,1],fossa_points[,2],col="blue",cex=0.2)NOTE: the arguments pch above and cex here are used simply to define the size and tipe of the points. You can play with this!
Depending on the species, the map/points can be too large or be outside our region of interest. We can focus on a particular region/country. Let’s do this for Madagascar
#Firts, get the "position" of Madagascar in the "wrld_simpl" object
which(wrld_simpl$NAME=="Madagascar")
#Then, use these indices (numbers) to subset the original object
mdg.nmb_map <- wrld_simpl[108,]Now you can repeat the steps above (plotting) for the selected region (Madagascar) only.
Check if there are any duplicated points
fossa_dups <- duplicated(fossa_points)
### NOTE: the function "duplicated" returns the results of a logical test (e.g. TRUE or FALSE)
# How many are duplicates?
length(which(fossa_dups==TRUE))
# How many are NOT duplicates?
length(which(fossa_dups==FALSE))
# Keep only those lines that are not duplicates
fossa_dups_row <- which(fossa_dups==TRUE)
# What's the size? That is, how many points are duplicates
length(fossa_dups_row)
# Create another object withoyt the duplicate records
fossa_nodups <- fossa_points[-fossa_dups_row,]
# What are the dimensions of the new object?
dim(fossa_nodups)
# Take a look at the first rows of data
head(fossa_nodups)Check if there are data with longitude/latitude of zero
fossa_lonzero = subset(fossa_nodups, lon==0)
fossa_latzero = subset(fossa_nodups, lat==0)NOTE: Remember that variable names can be different (Latitude, latitude, lat, etc…)
Check if there are data with NAs
fossa_nas <- which(is.na(fossa_nodups[,1])==TRUE)
fossa_points_nonas <- fossa_nodups[-fossa_nas,]
#how many NAs?
dim(fossa_nodups)
dim(fossa_points_nonas)Plot with other maps, using other data (e.g. download other species using another package: “rgbif”).
civet <- occ_data(scientificName = 'Fossa fossana', hasCoordinate = TRUE)NOTA: You can choose any other species (try!)
Get a world map
world = map_data("world")Plot the points using the map as a polygon and defining the points’ attributes (colour, size, etc)
ggplot(world, aes(long, lat)) +
geom_polygon(aes(group = group), fill = "white",
color = "gray40", size = .2) +
geom_jitter(data = civet$data,
aes(civet$data$decimalLongitude, civet$data$decimalLatitude), alpha=0.6,
size = 4, color = "red") +
labs(title = "Malagasy civet")Is it good? What if we use another species?
Now with GoogleMaps
ayeaye = occ_data(scientificName = 'Daubentonia madagascariensis', hasCoordinate = TRUE)
# Define the "extent", the map's "window"
e = extent(43, 51, -26, -12)
# Take that window from GoogleMaps
r = gmap(e)
# plot
plot(r, interpolate=TRUE, main="Map")
# join with the lat/long values
xy1=cbind(ayeaye$data$decimalLongitude,ayeaye$data$decimalLatitude)
# Add the points, using the same projection as the map (i.e. "Mercator")
points(Mercator(xy1) , col='blue', pch=20)
# Write the species name
text(4500000,-2400000,"\n\n\nDaubentonia \nmadagascariensis", adj = c(0,1), col="red", cex=0.75)This part can be used to create “simple” range maps based on geometry (e.g. minimum convex polygons, etc.), without considering environmental variables (no ENMs or SDMs).
This model draws a convex hull around all ’presence’ points.
#create a polygon around the species' records
fossa_hull <- convHull(fossa_points_nonas)
#Plot the created polygon
plot(mdg.nmb_map)
plot(fossa_hull,col="green",add=T)
#Create a raster to be used as a mask
r <- raster()
values(r) <- 0.5
ext <- c(43, 51, -26, -12)
r1 <- crop(r,ext)
#Use the raster and the polygon to create a new "map"
fossa_hullmodel <- predict(fossa_hull, r1)
plot(fossa_hullmodel, main='Convex Hull')
plot(mdg.nmb_map, border='dark grey',add=T)
points(fossa_points_nonas, pch='+')This model draws circles around all ’presence’ points.
fossa_circ <- circles(fossa_points_nonas, lonlat=TRUE)
fossa_circles <- predict(r1, fossa_circ, mask=TRUE)
plot(fossa_circles, main='Circles')
points(fossa_points_nonas, pch=19, cex=0.5, col="red")
#Play with the distance to create the "circles"
fossa_circ2 <- circles(fossa_points_nonas, d=100000, lonlat = TRUE)
fossa_circles2 <- predict(r1, fossa_circ2, mask=TRUE)
plot(fossa_circles2, main='Circles')
points(fossa_points_nonas, pch=19, cex=0.5, col="red")NOTA: Save the R “workspace”. That is, everything that was done in these exercises must be saved in an “.Rdata” file. These objects may be used later on. DO NOT CLOSE R IF WITHOUT SAVING THE FILE!!