Extracting raster data

R Skill Level: Intermediate - this activity assumes you have a working knowledge of R

This activity incorporates the skills learned in the Points, Polygons, Raster & Projections activities.

Objectives & Goals

Upon completion of this activity, you will:
  • know how to generate isoptope assignments
  • know how to perform basic raster calculations

Required packages

To complete the following activity you will need the following packages installed: raster sp rgeos

Installing the packages

If installing these packages for the first time consider adding dependencies=TRUE install.packages("raster") install.packages("velox") install.packages("sp") install.packages("rgeos")

Download R script Last modified: NA


Gathering spatial data

First we need a raster with the predicted stable-isotope of Hydrogen in precipitation. We can get that from waterisotopes.org. Below we use R to download the file and unzip it for us.

# Download the file #
utils::download.file(url = "http://wateriso.utah.edu/waterisotopes/media/ArcGrids/AnnualD.zip",
                     destfile = "../Spatial_Layers/AnnualD.zip",
                     quiet = TRUE,
                     extra = getOption("download.file.extra"))

# unzip the downloaded file #
utils::unzip(zipfile = "../Spatial_Layers/AnnualD.zip",
            files = NULL, list = FALSE, overwrite = TRUE,
            junkpaths = FALSE, exdir = "../Spatial_Layers", unzip = "internal",
            setTimes = FALSE)

# remove the zipped file to by tidy
## [1] TRUE

Read in the predicted surface for Hydrogen isotopes in precipitation. Generally, in North America Hydrogen isotope values are more depleted (more negative) in the north and more enriched (less negative). The plot below shows that pattern.

# read in Mean Annual D #
m_a_d <- raster::raster("../Spatial_Layers/AnnualD/mad/w001001.adf")

plot of chunk unnamed-chunk-4

Now that we have the base layer we need to determine the discrimination factor between the isotope values precipitation and the measured isotope values in the tissue of interest. Here, we’ll simulate some data using techniques we learned in other lessons.

back to top

Simulate ‘known’ values from capture location

To simulate these data we’ll first need capture locations that fall within the region of interest. For this exercise we’ll simulate 10 captures at 10 sites within the United States. Once we have the capture locations, we will extract the isotope value in precipitation. We’ll then simulate ‘measured’ isotope values - we’ll add a little noise to be more realistic.

# Read in a shapefile of the United States #
States <- raster::getData("GADM", country = "United States", level = 1)

# remove Alaska and Hawaii for this simulation (sorry) #

States <- States[!(States$NAME_1 %in% c("Alaska","Hawaii")),]

# dissolve all states so only have the border
USborder <- rgeos::gUnaryUnion(States, id = States$ISO)

Now that we have a boundary for the contiguous United States we can generate 10 randomly generated capture locations.

# Set the random seed generator

# Generate capture locations for calibration
CalibCapLocs <- sp::spsample(x = USborder, n = 10, type = "random")

plot of chunk unnamed-chunk-7

Extract precipitation values

Now that we have capture locations we’ll extract the predicted isotope values in precipitation.

# extract estimated H in precipitation
estHp <- raster::extract(m_a_d, CalibCapLocs)

Here are the predicted Hydrogen values in precipitation.

##  [1]  -32.96  -34.04 -117.70 -109.95  -51.22  -31.44  -75.62  -44.69
##  [9]  -24.39  -80.26

Simulate ‘measured’ tissue values

If the tissue values matched the precipitation values exactly making assignments would be pretty straightforward. Here we’ll add a little realism by adding some noise to the estimated value plus account for isotopic discrimination as they move through the food web. We’ll simulate 10 individuals at every capture location.

# set seed for random sampling

# Make a data frame to store the data
Hvalues <- data.frame(Hp = rep(estHp,each =10), # H in precip
                      Ht = NA)                  # H in tissue

# known slope 
knownBeta <- 0.93
Beta.sd <- 0.1

# known intercept
knownAlpha <- -15
Alpha.sd <- 0.2

# y = mx+b with some wiggle 
Hvalues$Ht <- Hvalues$Hp*rnorm(knownBeta,n = 100, sd = Beta.sd)+rnorm(knownAlpha,n=100, sd = Alpha.sd)

plot of chunk unnamed-chunk-11

back to top

Discrimination factor

Now that we’re finished simulating data we can begin the analysis. First we’ll use the simulated data to generate a discimination factor. In other words we need to determine how the tissue isotope signatures are related to the precipitation signatures. We’ll do this using a linear model so that we can make geographic assignments for individuals from unknown origins.

# discrimination factor 
fit <- lm(Hvalues$Ht~Hvalues$Hp)
## Call:
## lm(formula = Hvalues$Ht ~ Hvalues$Hp)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.0870  -3.9547  -0.9672   4.7238  21.8448 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -14.57726    1.70911  -8.529 1.86e-13 ***
## Hvalues$Hp    0.96013    0.02505  38.329  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 8.031 on 98 degrees of freedom
## Multiple R-squared:  0.9375,	Adjusted R-squared:  0.9368 
## F-statistic:  1469 on 1 and 98 DF,  p-value: < 2.2e-16

below we extract the information we need to make isotope assignments. First, we’ll need the intercept and slope from the linear equation above.

# extract the slope from lm oject
slope <- fit$coefficients[2]

# extract intercept from lm object
int <- fit$coefficients[1]

back to top

create a tissue-isoscape

Now that we know the relationship between our tissue of interest (say a tail feather) and the precipitation isoscape we can create a ‘tissue-isoscape’ to make predictions.

# make the tissue-isocape using y = mx+b
Tissue_isoscape <- m_a_d*slope+int

Great - now we’ll use this discrimination factor to make predictions for 100 individuals captured in the non-breeding season. First we need to simulate those values. In this simulation we’ll assume weak migratory connectivity where a single capture location in the non-breeding season is home to individuals from across the breeding distribution - similar to the Prothonotary warbler (see Tonra et al. 2019).

# simulate 100 tail feather values 
unk_Ht <- runif(n = 100, min = min(Hvalues$Ht), max = max(Hvalues$Ht)) 

plot of chunk unnamed-chunk-17

back to top

Making geographic assignments

Now that we 1) know the discrimination factor and 2) have a tissue-isoscape we are almost ready to make geographic assignments. First we need to understand what we are about to do. Below we make probabilistic geographic assignments based on a normal probability density function. The mean of the normal is the ‘measured’ tissue value. The location uncertainty is determined by the standard deviation. We derive the standard deviation using the standard deviation of the residuals from our discrimination function above.

Let’s extract the standard deviation of the residuals first.

(StdDev <- sd(fit$residuals))
## [1] 7.990335

Below is a figure of the general concept for how we make assignments.

plot of chunk unnamed-chunk-19

Now let’s apply the normal probability density function to a raster so that we can make spatially explicit assignments. The spatially explicit assignments are made where the likelihood that each ‘unknown’ tissue value, $y^*$, originates from a given location is

where $\mu_b$ is the specific cell in a given tissue isoscape and $\sigma_b$ is the standard deviation of the residuals from the discrimination equation.

Let’s make a function that we’ll call assign that makes the assignments for us.

# spatially explicit assignment function
assign <- function(isotope,isoscape,std) {
          ((1/(sqrt(2 * 3.14 * std))) * exp((-1/(2 * std^2)) * ((isotope) - isoscape)^2))

Let’s make our first geographic assignments

# Let's get an approximate idea where the 
# first individual molted 

MoltEst <- assign(isotope = unk_Ht[1],
                   isoscape = Tissue_isoscape,
                   std = StdDev)

plot of chunk unnamed-chunk-22

Great - now we have our first isotope assignment but we make predictions for the entire planet. Let’s confine our estimates to include only North America since our ‘study’ species isn’t found anywhere else during the breeding season.

plot of chunk unnamed-chunk-23

Coming soon

Assignments using the MigConnectivity package

back to top