Using Domain to model Butterfly Distribution in Pakistan

1. Introduction

In this tutorial basic features of the algorithm Domain will be presented by modeling a prediction of the specie´s distribution of butterflies for the country Pakistan. The tutorial is part of the project seminar “Species Distribution Modeling” at Philipps-University Marburg.

2. Theory of Domain

Domain is an profile method that was developed by Carpenter G., A.N. Gillison and J. Winter in 1993. It is an algorithm that computes the dissimilarity between a test location and a location of known occurrence of a certain species (Carpenter et al. 1993). Therefor it only needs presence records and some biophysical attributes (ebd.). Same can be said about Bioclim, a profile method that was published in 1986 already. In the eyes of Carpenter et al. Bioclim had some deficiencies. They criticized that the model treats each climatic axis independently. This result can lead to problematic predictions, as in the following plot:

Figure 1. Bioclim environmental envelope (Carpenter et al. 1993).

While point A, which is a know occurrence site, is outside of core bioclimate, point B is inside. However Point B is not very similar to the actual occurrence sites. The following is the adaptation of Carpenter et al. known as Domain.

Figure 2. DOMAIN environmental envelope (Carpenter et al. 1993).

Note that the Points A and C are in the core bioclimate, while the Point B is not. In general the prediction is more precise.

In order to test for the similarity of locations, various abiotic factors are taken into account. These factors however can be qualitative as well as quantitative. In order to work with mixed data types the algorithm uses the so-called Gower Distance.

What is the Gower Distance?

The Gower distance does not describe a spatial distance but the similarity of two locations to each other. The data types are combined as follows:

“Briefly, to compute the Gower distance between two items you compare each element and compute a term. If the element is numeric, the term is the absolute value of the difference divided by the range. If the element is non-numeric the term is 1 if the elements are different or the term is 0 if the elements are the same. The Gower distance is the average of the terms.” (McCaffrey 2020)


We have a dataset of a group of 5 people described by 5 different elements: age, height, gender, eye color and do they have pets. This data set could look like this:

##   age height gender eyecolor pets
## 1  12    161   male    green  YES
## 2  31    180 female     blue   NO
## 3  61    163 female    brown  YES
## 4  49    189   male     blue   NO
## 5  25    172 female    brown   NO

If we were interested in the similarity of person #2 and #5 we would get one value between 0 and 1 for each element. As McCaffrey states, the average of these values is the Gower Distance of the people.

##   age height gender eyecolor pets
## 2  31    180 female     blue   NO
## 5  25    172 female    brown   NO
diff(range(group$age))          # range of age 
## [1] 49
(31-25)/49                      # age value
## [1] 0.122449
diff(range(group$height))       # range of height
## [1] 28
(180-172)/28                    # height value
## [1] 0.2857143
#0                              # gender value
#1                              # eye color value
#0                              # having pets value
mean(0.122449,0.2857143,0,1,0)  # Gower Distance between person #2 and #5
## [1] 0.122449

3. Preparation and loading data packages

## Loading required package: sp
## Checking rgeos availability: TRUE
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/super/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/super/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/super/Documents/R/win-library/4.0/rgdal/proj
## Loading required package: lattice
## Loading required package: ggplot2
pak <- getData("GADM", country = "Pakistan", level = 0)

4. Data input

In order to fly efficiently butterflies need to maintain the temperature of their flight muscles well above the ambient temperature (Turner et al. 1987). An increased thermal and radiant energy significantly also increases their reproductive rate (Kingsolver 1983). Furthermore high levels of humidity increases the species richness of butterflies. In an multi-species analysis the observed species richness in a wet forest was two times higher than in a comparable dry one (Checa et al. 2019).

This input allows us to incorporate various climatic variables. For this tutorial we will choose the following:

  1. precipitation (mm)

  2. solar radiation (kJ m-2 day-1)

  3. average temperature (°C)

  4. maximum temperature (°C)

  5. minimum temperature (°C)

  6. water vapor pressure (kPa)

The used Data can be found on the worldclim page.

Preparing the environmental data

#they need to be stacked together and cropped to the extent of Pakistan.
wprec <- list.files("C:/Users/super/Documents/sdm/Domain/Process/Data/envdata/wc2.1_2.5m_prec",pattern = ".tif", full.names = TRUE) # This creates a vector with all tif-files in this location
wprecstack <- raster::stack(wprec)
wsrad <- list.files("C:/Users/super/Documents/sdm/Domain/Process/Data/envdata/wc2.1_2.5m_srad",pattern = ".tif", full.names = TRUE)
wsradstack <- raster::stack(wsrad)
wtavg <- list.files("C:/Users/super/Documents/sdm/Domain/Process/Data/envdata/wc2.1_2.5m_tavg",pattern = ".tif", full.names = TRUE)
wtavgstack <- raster::stack(wtavg)
wtmax <- list.files("C:/Users/super/Documents/sdm/Domain/Process/Data/envdata/wc2.1_2.5m_tmax",pattern = ".tif", full.names = TRUE)
wtmaxstack <- raster::stack(wtmax)
wtmin <- list.files("C:/Users/super/Documents/sdm/Domain/Process/Data/envdata/wc2.1_2.5m_tmin",pattern = ".tif", full.names = TRUE)
wtminstack <- raster::stack(wtavg)
wvapr <- list.files("C:/Users/super/Documents/sdm/Domain/Process/Data/envdata/wc2.1_2.5m_vapr",pattern = ".tif", full.names = TRUE)
wvaprstack <- raster::stack(wvapr)
# In a step in between the data got cropped and clipped to the extent of pak
cropped_prec <- crop(wprecstack, extent(pak))
clipped_prec <- mask(cropped_prec, pak)
clipped_prec_sum <- sum(clipped_prec)
cropped_rads <- crop(wsradstack, extent(pak))
clipped_rads <- mask(cropped_rads, pak)
clipped_rads_mean <- mean(clipped_rads)
cropped_tavg <- crop(wtavgstack, extent(pak))
clipped_tavg <- mask(cropped_tavg, pak)
clipped_tavg_mean <- mean(clipped_tavg)
cropped_tmax <- crop(wtmaxstack, extent(pak))
clipped_tmax <- mask(cropped_tmax, pak)
clipped_tmax_mean <- mean(clipped_tmax)
cropped_tmin <- crop(wtminstack, extent(pak))
clipped_tmin <- mask(cropped_tmin, pak)
clipped_tmin_mean <- mean(clipped_tmin)
cropped_vapr <- crop(wvaprstack, extent(pak))
clipped_vapr <- mask(cropped_vapr, pak)
clipped_vapr_mean <- mean(clipped_vapr)
wvaprstack <- raster::stack(wvapr)
cropped_vapr <- crop(wvaprstack, extent(pak))
clipped_vapr <- mask(cropped_vapr, pak)
clipped_vapr_mean <- mean(clipped_vapr)

A write out of our now stacked raster data looks like this.

pakclim <- stack(clipped_prec_sum,clipped_rads_mean,clipped_tavg_mean,clipped_tmax_mean,clipped_tmin_mean,clipped_vapr_mean)
plot(pakclim,main=c("precipitation","solar radiation","average temperature","maximum temperature","minimum temperature","water vapor pressure"))

Now the butterfly data are read in and converted into spacial points.

butterflies <- read.csv("C:/Users/super/Documents/sdm/Domain/Process/Data/PakistanLadakh.csv", header=TRUE, sep=',')
butterflies[,1] <- as.factor(butterflies[,1])
list_butterflies <- levels(butterflies[,1]) 
coordinates(butterflies) <- c("x","y")
points <- crop(butterflies, pak)
plot(pakclim$layer.1, xlab = "Latitude", ylab = "Longitude", main= "Precipitation Sum", 1)
plot(pak, add=TRUE)
points(butterflies, pch=19, cex=0.1, col='black')