Using remote sensing to map invasive species

Success in the strategies to manage biological invasions rely in our ability to detect and monitor species distributions. However, the field surveys designed to detect and map invasive species often take enormous logistic efforts, in time and space, and might lead to observational bias (e.g. decreased sampling effort at less accessible areas) which can negatively affect the outcome and its applicability.

Remote sensing provides standardized large-scale spatially and temporally explicit data to map and monitor changes in the environment and the distribution of native and invasive species. It helps field surveys overcome the logistical barrier and facilitates the detection and monitoring of species at multiple temporal and spatial scales. Here, as part of the DIARS project, we present a method to map invasive species using hyperspectral images.

The dataset

In this tutorial we will map the distribution of a invasive bryophyte species (Campylopus introflexus), based on the study presented in Skowronek et al. 2016 at Sylt island (Germany) as part of the DIARS project.

For the purpose of this tutorial and in order to reduce the computer processing time (while using real remote sensing data in the analyzes), we provide an already prepared artificial image (click HERE to read more information on the process to create this image).

Throughout the tutorial you will find this gray boxes. They contain general “best practices” tips that will help you improve the quality of the analysis when applying this framework to your own study case.

The framework

Figure 1. Framework used to map invasive species distribution using remote sensing data.

You can find more information on the software we will be using in this tutorial HERE.

Let’s start!

Setting the R environment

Let’s start by loading the libraries that we are going to be using in this tutorial. To load the packages, it is necessary to first install them using the command install.packages. After installation you can load the packages in R:

# Run this command if the packages are not already installed
# install.packages(c('raster','rgdal','hyperSpec','rgeos','rgrass7'))

# Load the packages
library(raster)
library(rgdal)
library(hyperSpec)
library(rgeos)
library(rgrass7)

Now create a folder where the datasets are going to be located and the outputs and results will be saved. To set the working directory in R to the folder created, execute the following command line, replacing the path to the folder you just created:

setwd("/home/garzonc/Desktop/DIARS")

Back to the top

Data preparation

Three datasets are required to map the species through remote sensing: i) two shapefiles with the x,y locations of the field plots that correspond to the calibration and validation points, ii) a shapefile containing a set of points randomly distributed within the study area that do not overlap with the field plots and will serve as background points and, iii) a hyperspectral image in tiff format. You can download a sample dataset HERE and save them in the working directory you just created.

Remember to write down the resolution, projection and geographic region of all your shape and raster files. Make sure they match!

Setting up the field and remote sensing based-data

The first step in the mapping process is to link the field data with the hyperspectral data, that is, to extract the hyperspectral values (reflectances at different wavelenghts) at each x,y location of the validation, calibration and background points. Let’s start!

To extract the spectral values for calibration, validation and background points we first need to load the hyperspectral image in R. Here we also plot the hyperspectral imagery and the field plots to check that all of the field plots are located within the hyperspectral image.

# Create a RasterStack object with the hyperspectral bands in tiff format
layers <- stack("mosaic_datasetSyltMix.tiff")
# Explore and visualize the bands in the RasterStack
layers

## class       : RasterStack 
## dimensions  : 180, 175, 31500, 248  (nrow, ncol, ncell, nlayers)
## resolution  : 1.8, 1.8  (x, y)
## extent      : 469443.6, 469758.6, 6075848, 6076172  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : mosaic_datasetSyltMix.1, mosaic_datasetSyltMix.2, mosaic_datasetSyltMix.3, mosaic_datasetSyltMix.4, mosaic_datasetSyltMix.5, mosaic_datasetSyltMix.6, mosaic_datasetSyltMix.7, mosaic_datasetSyltMix.8, mosaic_datasetSyltMix.9, mosaic_datasetSyltMix.10, mosaic_datasetSyltMix.11, mosaic_datasetSyltMix.12, mosaic_datasetSyltMix.13, mosaic_datasetSyltMix.14, mosaic_datasetSyltMix.15, ...

plot(layers[[150]], legend = FALSE)  #Plot band 150

To read the point shapefiles make sure they are located in the current working folder. Initiate the GRASS software from within R, create a directory that will contain the GRASS GIS database for the working session, set the gisBase argument to GRASS binaries and libraries directory path, this should be something like C:\Program Files\GRASS GIS 7.00beta4 if your are using Windows (check HERE for other operating systems) and, set the correct projection, geographic region and spatial resolution before opening any shapefile:

dir.create("grassdata")  #Create a grassdata directory

myGRASS <- "/home/garzonc/src/grass72/grass72_release/dist.x86_64-pc-linux-gnu"  #Set the location of the grass binaries and libraries
myGISDbase <- "/home/garzonc/Documents/grassData"  #Set the directory path to your GISDbase
myLocation <- "DIARS"  #Name the Location
myMapset <- "PERMANENT"  #Name the mapset

# Initiate the GRASS GIS session
initGRASS(myGRASS, home = tempdir(), gisDbase = myGISDbase, location = myLocation, 
    mapset = myMapset, override = TRUE)

## Error in initGRASS(myGRASS, home = tempdir(), gisDbase = myGISDbase, location = myLocation, : could not find function "initGRASS"

# Set the projection
execGRASS("g.proj", flags = c("c"), parameters = list(epsg = 32632))

## Error in execGRASS("g.proj", flags = c("c"), parameters = list(epsg = 32632)): could not find function "execGRASS"

# Set the region and resolution
execGRASS("g.region", parameters = list(n = "6080248.2", s = "6079942.2", e = "464307.5969", 
    w = "463974.5969", res = "1.8"))

## Error in execGRASS("g.region", parameters = list(n = "6080248.2", s = "6079942.2", : could not find function "execGRASS"

Now the R environment is set and you can import the calibration points shapefile using the v.in.ogr GRASS module:

execGRASS("v.in.ogr", flags = c("overwrite", "quiet"), parameters = list(input = "/home/garzonc/Desktop/DIARS/cal_pointsR.shp", 
    output = "cal_pointsR"))

pointdata_cal <- readVECT("cal_pointsR")  #Read the imported vector file

## Exporting 90 features...
##    5%  11%  17%  23%  30%  36%  42%  48%  54%  60%  66%  72%  78%  84%  90%  96% 100%
## v.out.ogr complete. 90 features (Point type) written to <cal_pointsR>
## (SQLite format).
## OGR data source with driver: SQLite 
## Source: "/home/garzonc/Documents/grassData/DIARS/PERMANENT/.tmp/localhost.localdomain/790.0", layer: "cal_pointsR"
## with 90 features
## It has 9 fields

layers <- stack("mosaic_datasetSyltMix.tiff")
plot(layers[[150]], legend = FALSE)  #Plot hyperspectral image band 150
points(pointdata_cal, pch = 3, cex = 0.5)  #Add the calibration points to the plot