An Introduction to species distribution models

Load the packages we will need (if a package is still not in installed, use install.packages(“PackageName”) to install it).

library(rworldmap);library(rgdal);library(raster);library(plotfunctions);library(bRacatus);library(rgbif)

## Loading required package: sp

## ### Welcome to rworldmap ###

## For a short introduction type :   vignette('rworldmap')

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson

Download species data from GBIF

koala <- occ_search(scientificName="Phascolarctos cinereus",
                      hasCoordinate=T)

koala_pts <- koala[[3]]

koala_pts2 <- giveOcc(koala_pts, 
                         longitude = "decimalLongitude",
                         latitude = "decimalLatitude")

coordinates(koala_pts) <- ~decimalLongitude+decimalLatitude  #inform coordinates
  
#visualise
plotOcc(koala_pts2)

## Warning in wkt(obj): CRS object has no comment

download bioclimatic variables (raster stack)

bioclim <- getData("worldclim",var="bio",res=10)

## Warning in getData("worldclim", var = "bio", res = 10): getData will be removed in a future version of raster
## . Please use the geodata package instead

class(bioclim)

## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

nlayers(bioclim)

## [1] 19

bioclim[[1]]

## class      : RasterLayer 
## dimensions : 900, 2160, 1944000  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio1.bil 
## names      : bio1 
## values     : -269, 314  (min, max)

anual_mean_temp <- bioclim[[1]]

class(anual_mean_temp)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

plot(anual_mean_temp)

Crop the raster to the area o study

extent(anual_mean_temp)

## class      : Extent 
## xmin       : -180 
## xmax       : 180 
## ymin       : -60 
## ymax       : 90

extent(koala_pts)

## class      : Extent 
## xmin       : 135.5122 
## xmax       : 153.5639 
## ymin       : -38.8 
## ymax       : -17.41922

study_area <- crop(anual_mean_temp,extent(koala_pts))

plot(study_area)
plot(koala_pts,add=T,pch=19,cex=.5)

Visualising another layer of the raster stack, and crop the whole area of Australia

anual_prec <- bioclim[[12]]

plot(anual_prec)

australia_prec <- crop(anual_prec,extent(109,156,-44,-10))
australia_temp <- crop(anual_mean_temp,extent(109,156,-44,-10))

plot(australia_prec)
plot(koala_pts,add=T,pch=19,cex=.5)

Extracting values from rasters at point locations

australia_var <- stack(australia_prec,australia_temp)
australia_var

## class      : RasterStack 
## dimensions : 204, 282, 57528, 2  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : 109, 156, -44, -10  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : bio12, bio1 
## min values :   129,   51 
## max values :  4234,  293

values <- extract(australia_var,koala_pts)
head(values)

##      bio12 bio1
## [1,]  1120  109
## [2,]   527  164
## [3,]  1501  201
## [4,]   862  143
## [5,]   698  154
## [6,]  1284  199

class(values)

## [1] "matrix" "array"

values <- as.data.frame(values)
class(values)

## [1] "data.frame"

min_temp <- min(values$bio1)
max_temp <- max(values$bio1)

min_prec <- min(values$bio12)
max_prec <- max(values$bio12)

Produce a very simple SDM for the Koala

climatic_suitability <- australia_prec

for(i in 1:length(climatic_suitability[]))
{
  climatic_suitability[i] <- ifelse(australia_temp[i] > min_temp &
                            australia_temp[i] < max_temp &
                            australia_prec[i] > min_prec &
                            australia_prec[i] < max_prec,
                            1,0)
  
}

plot(climatic_suitability)