Ecosystem impact of invasive species

Biological invasions not only affect biodiversity and ecosystem structure they also impact ecosystem functioning via changes in biochemical pools and nutrient cycling, threatening the provision of ecosystem services (Ehrenfeld et al. 2010). Observational and experimental studies have already shown the effect of plant invasions (Vitousek et al. 1987, Liao et al. 2008) and the challenge now is to understand how invasions affect the ecosystems from a large scale perspective and how this knowledge can be used to inform management and conservation decisions.

Now, using plant functional traits as key indicators of ecosystem functioning (Chapin 2003) in combination with newly developed remote sensing tools to measure and predict such traits (Homolova et al. 2013), we can finally evaluate the impact of biological invasions, scale it from individuals to landscape and, generate spatially explicit assessments of their risks and consequences.

In this tutorial we will explore the use of hyperspectral data in combination with ground data (describing a selected set of biochemical leaf traits) of native tree species to assess the impact of an invasive shrub (Prunus serotina) on Phosphorus and Nitrogen foliar concentrations of the native vegetation.

The Dataset
The method and dataset presented here are thoroughly explained in Ewald et al. (in prep), a publication of the DIARS project. The dataset consists of an artificial hyperspectral image (click HERE to learn more about the artificial image) composed of a series of 150 field plots and 460 randomly chosen background plots accompanied by a table with P and N values for each of the field plots and the percentage cover of the invasive shrub Prunus serotina. You can download the data HERE, and remember to save it in the same directory you will use for your analysis.

The method

Lets start by installing and loading the packages with the commands for image processing and (raster and rgeos) and modeling (autopls) and, set the working directory:

# Run this command to install packages
# install.packages('raster','rgeos','autopls')

library(raster)
library(rgeos)
library(autopls)

# Set working directory
setwd("/home/garzonc/Desktop/DIARS/")

Once the working directory is set we can import the hyperspectral image and apply a brightness normalization which reduces the effects of shadows in the image. After that process is finished, we must aggregate the resulting image from its original size 3x3 meters per pixel to the new 12x12 meters per pixel resolution.

# import hyperspectral image
r.hs <- stack("mosaicCo_datasetComp.tiff")
# apply brightness normalisation
r.hs.bn <- prepro(r.hs, method = "bn")
# aggregate raster (using a factor of 4) for the prediction maps later.
r.hs.bn.agg <- aggregate(r.hs.bn, fact = 4)

Let’s import the field data. The first table contains the invasive species cover, recorded for each forest layer (herb, shurb and tree) and the geographic location of the plots (grid_allData1.csv). The second table contains the biochemical foliar measurements which consist of community weighted mean values of leaf nitrogen and phosphorous content for each plot (leafNP.csv).

# read field data
plotdata <- read.csv("grid_allData1.csv")
# plot band 100 of the hyperspectral image
plot(r.hs, (100), legend = FALSE)
# add the plot locations to the figure
points(plotdata$xgrid, plotdata$ygrid, pch = 3, cex = 0.5)

# read field data
np <- read.csv("leafNP.csv")

The np table we just imported contains the Phosphorus and Nitrogen concentrations per plot. We must now extract the mean wavelength value per plot from the hyperspectral images. To do so, create a spatial points object using the plot location data, generate a buffer of 24 meters around each point (to match the pixel size of the aggregated image), using that buffer extract the mean hyperspectral value for each band, and finally join both datasets in one table:

# create a SpatialPoints file
plots <- SpatialPoints(cbind(plotdata[1:50, "xgrid"], plotdata[1:50, "ygrid"]))
# create a buffer area
plots <- gBuffer(plots, byid = T, width = 12, capStyle = "SQUARE")

# extract the mean hyperspectral value for each band
extr.val <- extract(r.hs, plots, fun = mean, na.rm = TRUE)

# join the datasets
plotdata1 <- merge(np, plotdata, all.x = TRUE)

# Reorder the rows according to the grid number
plotdata1 <- plotdata1[order(plotdata1$gridN), ]

Now have 248 variables (the 248 wavelength bands of the hyperspectral image) that could potentially explain the concentrations of Nitrogen and Phosphorus on the vegetation and relate the patterns in P and N content to the Prunus serotina cover. A versatile method to analyze such multivariate data is known as partial least squares regression (PLSR). This method addresses the problem of making predictions in the presence of multivariate data. The package autopls contains a function using PLSR with a backward selection of predictors approach. To applied this method in the analysis of the data we have to follow these steps:

# calculate PLSR model for P and N selecting cross validation (CV) as the
# validation method, this method allows running a 10-fold cross validation
# calculate PLSR model for Nitrogen
plsrN <- autopls(plotdata1$N[1:50] ~ extr.val, val = "CV")
## 1   Pred: 248  LV: 2   R2v: -0.047 RMSEv: 2.355  
## 2   Pred: 58   LV: 2   R2v: -0.171 RMSEv: 2.49   Criterion: A1
## 3   Pred: 52   LV: 2   R2v: -0.1   RMSEv: 2.414  Criterion: A4
## 
## Predictors: 248   Observations: 50   Latent vectors: 2   Run: 1 
## RMSE(CAL): 2.16   RMSE(CV): 2.35   
## R2(CAL): 0.123    R2(CV): -0.0471
# calculate PLSR model for Phosphorus
plsrP <- autopls(plotdata1$P[1:50] ~ extr.val, val = "CV")
## 1   Pred: 248  LV: 2   R2v: -0.255 RMSEv: 0.28   
## 2   Pred: 223  LV: 2   R2v: -0.239 RMSEv: 0.279  Criterion: A4
## 3   Pred: 200  LV: 3   R2v: -0.173 RMSEv: 0.271  Criterion: A4
## 
## Predictors: 200   Observations: 50   Latent vectors: 3   Run: 3 
## RMSE(CAL): 0.241   RMSE(CV): 0.271   
## R2(CAL): 0.0725    R2(CV): -0.173
# create prediction map for P
Pmap <- predict(plsrP, r.hs.bn.agg)

# plot the map
plot(Pmap)

Is important to clarify that the green/yellow area in the predicted map based on the PLSR correspond to an area of a small lake and although is visible in the plot it is not included in the analysis. Now we can extract the validation plot values to validate the model outcome and examine if the concentration of P and N (as presented in the Pmap) relates to the Prunus serotina presence:

# extract values for validation points
val.plots <- plotdata[grep("validation", plotdata$Plot1), ]
val.plots$Ppred <- extract(Pmap, cbind(val.plots$xgrid, val.plots$ygrid))

# plot P values in ares with and without Prunus presence create new
# category: Prunus present in shrub or tree layer
val.plots$Prunus_pres <- factor(ifelse(val.plots$Prunus_tre + val.plots$Prunus_shr > 
    0, 1, 0))
plot(Ppred ~ Prunus_pres, data = val.plots, ylab = "Phosphorus values predicted", 
    xlab = "Prunus serotina presence")

The plot above shows the phosphorus (P) values predicted in the case of presence (1) or absence (0) of the invasive shrub Prunus serotina