GBIF Occurrence Data

We’ll start by sourcing our data. For observations, let’s take a look at Purple Loosestrife, a wetland species that is native to Europe, and invasive in North America. For actual research work, I normally download the files directly from GBIF, and examine them carefully to check for errors or missing data. For this demo we’ll use the rgbif package to download the data directly into R, and we’ll assume there are no problems with the data that need to be corrected.

lsGBIF <- occ_search(scientificName = "Lythrum salicaria",
                    limit = 10000,
                    basisOfRecord = "Preserved_Specimen",
                    hasCoordinate = TRUE,
                    fields = c("decimalLatitude",
                               "decimalLongitude", "year",
                               "country", "countryCode"))

save(lsGBIF, file = "../data/2021-07-29-ls-gbif-recs.Rda")

This returned an object with 7969 records. I saved that locally, so that I’m not making GBIF search their database everytime I work on this demo.

load("../data/2021-07-29-ls-gbif-recs.Rda")
lsOccs <- lsGBIF$data

lsGBIF$data is the table with the actual records in it. That’s what we’ll be working with. The other components of lsGBIF are metadata related to the original GBIF search. That’s useful to have, but not needed for the rest of this example.

Next, we tell R which columns are the coordinates, which allows us to map the observations. This also converts our observation matrix to a SpatVector object.

The crs argument here tells R that our points are in lat/lon (unprojected) coordinates. If all of your data use lat/lon, you don’t need to specify this, but it’s important if you need to reproject your data, or combine data with different projections.

lsOccs <- vect(lsOccs, geom = c("decimalLongitude",
                                "decimalLatitude"),
               crs = "+proj=longlat +datum=WGS84")

Basemap

We’ll also need a world map to use in our plots. the world function from the geodata package will download one for us. The first time you call this function in a directory, it downloads the data from the internet, and saves it locally according to the path argument. Subsequent calls will load your local copy of the data, to speed things up.

In this case I’ve set the path argument to store the downloaded files in a location that is convenient for me. You can set this to anything you like, or leave it out to use the current R working directory.

wrld <- world(path = "../data/maps/")

This is a relatively low resolution map of the continents. For higher resolution maps, see the gadm function.

Now we can plot our data. Note that we set the margins (via mar) inside the plot call, not via the par function used in most base R plots.

plot(wrld, border = "gray80", mar = c(0, 0, 0, 0))
points(lsOccs, col = 2, cex = 0.3)

Climate Data

To get our climate data, we can use geodata’s worldclim_* functions. I’m using the coarsest resolution (10 minutes) to speed things up for this demonstration. The path argument works the same was as for world, storing a local copy of the files. In this instance we’ll use worldclim_global to get the bioclim variables for the world all at once:

wclim <- worldclim_global(var = "bio", res = 10,
                          path = "../data/") 

We can take a look at one layer:

plot(wclim$wc2.1_10m_bio_1, main = "bio1",
     mar = c(0.1, 0.1, 2, 4), legend = TRUE,
     axes = FALSE)
par(mar = c(0.1, 0.1, 2, 4))
box()

Sampling Bias

One of the challenges we deal with with herbarium data is that observations tend to be clustered together in non-random ways. Sites near universities, museums, or in well-used field stations will have more records than more remote locations. We can’t completely account for this bias, but we can reduce it by spatial thinning. This is the process of randomly selecting a small number (often just one) of records for each grid cell in our analysis. The result is that those highly-sampled locations will be represented by a single record, meaning they will have a less exaggerated influence on the results.

We will need the full set of records below, so I’ll make a copy of the un-thinned data named lsOccsAll.

lsOccsAll <- lsOccs ## keep this for later

## select one record per cell:
lsOccs <- spatSample(lsOccs, size = 1, strata = wclim)

Linking Climate Data to Species Observations

Next, we need to extract the environmental values from the climate rasters for each of our observation records:

lsOccs <- cbind(lsOccs, extract(wclim, lsOccs))

In the process of extracting wclim values for our observations, we usually end up with a few missing values. This is a consequence of mismatches between the observation coordinates and the climate rasters. In some cases, the observations are placed off the coast in the ocean, or in another area where there is no climate data available. We need to exclude these missing values from our analysis.

lsOccs <- lsOccs[complete.cases(data.frame(lsOccs)), ]

PCA

The Niche Quantification analysis starts with a Principal Components Analysis of the environmental data. The actual ordination uses the global data, with the observation records and the native and invasive background environment treated as supplemental rows.

pca.clim <- dudi.pca(globalEnvM, center = TRUE,
                    scale = TRUE, scannf = FALSE, nf = 2)
global.scores <- pca.clim$li

nativeLS.scores <-
  suprow(pca.clim,
         data.frame(lsEA)[, colnames(globalEnvM)])$li   
invasiveLS.scores <-
  suprow(pca.clim,
         data.frame(lsNA)[, colnames(globalEnvM)])$li

nativeEnv.scores <- suprow(pca.clim, naEnvM)$li
invasiveEnv.scores <- suprow(pca.clim, eaEnvM)$li

Let’s break that down. dudi.pca does a PCA analysis on globalEnvM, which is a matrix of all the environmental variables over the entire study area. We use that to create a two-dimensional summary of the total environmental variability.

Next, we map our observation data (lsEA and lsNA) into that 2-dimensional ordination, using the suprow function. lsEA and lsNA are SpatialPointsDataFrame objects. Sometimes you can treat them as if they were data.frames, but other times you need to explicity convert them. This is one of those times, hence I’ve wrapped them in data.frame().

Recall that lsEA and lsNA have more columns than the environmental matrix: they also include year, countryCode, country. We only want to include the environmental variables when you project the observations into the ordination. To make sure that we use the same variables as in the original ordination of globalEnvM, in the same order, I select the columns explicitly to match that object:

data.frame(lsEA)[, colnames(globalEnvM)]

The output of dudi.pca and suprow includes a lot of information that we aren’t using here. We only need the li element, so I’ve selected that from each of the function outputs.

Occurrence Density Grids

Finally we’re ready to do the Niche Quantification/Comparisons. We’ll use the PCA scores for the global environment, the native and invasive environments, and the native and invasive occurrence records.

nativeGrid <- ecospat.grid.clim.dyn(global.scores,
                                   nativeEnv.scores,
                                   nativeLS.scores)

invasiveGrid <- ecospat.grid.clim.dyn(global.scores,
                                   invasiveEnv.scores, 
                                   invasiveLS.scores)

ecospat.plot.niche.dyn(nativeGrid, invasiveGrid,
                       quant = 0.05) 

The resulting plot shows us the environmental conditions present in Eurasia (inside the green line) and North America (inside the red line). The green area represents environments occupied by Lythrum salicaria in Eurasia, but not in North America, the red area shows environments occupied in North America and not Eurasia, and the blue area shows environments occupied in both ranges. We can also see that there are a few areas in Eurasia with environments not present in North America, and vice versa. However, for the most part, Lythrum salicara doesn’t occur in this environments.

To get the index values we use ecospat.niche.dyn.index:

indexVals <- ecospat.niche.dyn.index(nativeGrid,
                                     invasiveGrid) 
indexVals$dynamic.index.w

##   expansion   stability   unfilling 
## 0.005940337 0.994059663 0.022704396

Projecting Climate Space to Geographic Space

We can take these results and project them onto a geographic map. This will show us which areas of the native and invasive range of Lythrum salicaria fall into each of the three categories.

Note that the function ecospat.niche.dynIndexProjGeo doesn’t yet handle the SpatRaster objects created by the terra package, so we have to convert it to a raster’s stack when we call it. For convenience and consistency with the rest of the code, I convert the results back to a SpatRaster.

geoProj <-
  ecospat.niche.dynIndexProjGeo(nativeGrid,
                                invasiveGrid,
                                env = globalEnvR)

plot(geoProj, legend = FALSE, 
     col = c("grey", "green", "red", "blue"),
     ylim = c(20, 80), axes = FALSE, mar = c(0, 0, 0, 0))
points(lsEA, col = 'darkgreen',
       cex = 0.25, pch = 23, bg = "white")
points(lsNA, col = 'darkred',
       cex = 0.25, pch = 23, bg = "white")
par(mar = c(0,0,0,0))
box()