Perhaps the one question I’m asked the most is how to extract values from a raster using a series of points. Today, I answer that question. I’ll cover the basics of how to do it in R.

First, lets get you some data. Since we’re extracting from rasters, users typically have more than one raster that they want to extract from. I’ll only be using two, but you should be able to figure out how to extrapolate this out to n rasters. While I know the example data sets can be called within R, I’ll be pointing you to the rasters as a download. This is in the event that you’ve generated your own data rasters.

It doesn’t really matter what variables you pick. For simplicity’s sake, I’m choosing bioclim layers for average temperature and precipitation. Rasters are being downloaded from Worldclim. They are in the process of updating the website, so its possible the links will lead to a 404 error depending on when you’re reading this.

So here’s where it gets a little tricky. You need to know what libraries you’re going to be working with and what the dependencies are. If the library, like adehabitathr, requires sp, then you need to stick with library sp. However, the newest standard is sf, short for simple features. It definitely should be your preferred library, if possible, because its designed to work with the tidyverse. If you’ve never heard of the tidyverse, do yourself a favor and check out Hadley Wickham’s R for Data Science. Now, I will proceed with library sf as it is the successor to library sp.

library(sf)
library(dplyr)
library(raster)

With these libraries loaded, what you want to do next is to create a raster stack. Think of this as layering your rasters on top of each other, much like a stack of blankets. The stack function comes courtesy of library raster. But first, you’ll need to load the rasters. Since I don’t know which directory you downloaded your rasters to, I’ll use the file.choose() function to open up a window to allow you to select the rasters. An important note about using the raster::stack() function. The extent of all the rasters must be the same. If they are not the same, you will either have to mask the rasters so that they have the same extent or skip raster::stack altogether.

avg_temp <- file.choose()
precip <- file.choose()
all_rasters <- stack(raster(avg_temp), raster(precip))

#Or if you're not using raster::stack

raster_temp <- raster(avg_temp)
raster_precip <- raster(precip)

Now you’ll need to locate the point files that you want to use. For our purposes, we’ll create a quick dataframe. You, I suspect, will be loading data from a csv into a dataframe using read.csv or some other similar function. I’m just going to pick 2 cities (Santa Cruz, CA and Oxford, OH). In case you didn’t know, these were the cities where I went to undergrad and grad school respectively. A quick google search later tells me that Santa Cruz is at 36.974117, -122.030792, and Oxford, OH is at 39.506996, -84.745232. So let’s create this dataframe. Important note before we begin. Know what projection your points are in. Since these are latitude and longitude coordinates, they are unprojected and will occur in epsg 4326, hence why the argument crs in st_as_sf is set to 4326. If you do not know what projection your data is in, you will not get the correct data returned.

df <- data.frame(city = c("Santa Cruz", "Oxford"), lat = c(36.974117, 39.506996), long = c(-122.030792, -84.745232))
df_sf <- st_as_sf(x = df,
                  coords = c("long", "lat"),
                  crs = 4326)

Lets go over what we just did. I created a dataframe in R with 3 columns (city, lat, long). However, I needed it as a simple features point, so I used the function st_as_sf to convert it. Additional documentation on st_as_sf here. While you’ll have multiple options on how to do this, if you don’t know what to do, separate your latitude and longitude to different columns.

Before we begin our extraction, there is some data cleaning we need to do. Now, for the above example, the data arrives clean. However, for your data set, that is not likely. So what we’ll want to do is to remove all points that have empty geometries and remove all points that are classified as multipoint. The extract function only works with geometry type point. Additional information on all simple feature geometry types can be found here.

df_sf <- df_sf[!st_is_empty(df_sf),,drop = FALSE]
df_sf <- st_cast(df_sf, 'POINT')

The first line of the above block of code removes all empty geometries. The second line makes sure that every record is geometry type point and not geometry type multipoint. If you prefer piping like I do, then you can use the code block below. If you don’t know what piping is, read Hadley’s ebook referenced above and come back after. For now, ignore the following block of code if you don’t know what piping is. In short, the immediately above block of code and the immediately below block of code do the same thing but in different ways.

df_sf <- df_sf %>%
  filter(!st_is_empty(df_sf))

df_sf <- st_cast(df_sf, 'POINT')

Now its time to make the magic happen and extract our points. Its important to note that I want to be very specific and call the extract function from library raster, not library dplyr. So pay attention to the syntax.

values <- raster::extract(all_rasters, df_sf)

Now to finish it off, we want to know what the values actually are, so we combine them with the simple features collection we created above.

complete <- cbind(df_sf, values)

And thats it. If you made it this far, the complete code is:

library(sf)
library(dplyr)
library(raster)

avg_temp <- file.choose()
precip <- file.choose()
all_rasters <- stack(raster(avg_temp), raster(precip))

df <- data.frame(city = c("Santa Cruz", "Oxford"), lat = c(36.974117, 39.506996), long = c(-122.030792, -84.745232))
df_sf <- st_as_sf(x = df,
                  coords = c("long", "lat"),
                  crs = 4326)

#Pick one
#Either this one, with no pipes

df_sf <- df_sf[!st_is_empty(df_sf),,drop = FALSE] #drops all points that have empty geometries
df_sf <- st_cast(df_sf, 'POINT')

#Or this one with piping

df_sf <- df_sf %>%
  filter(!st_is_empty(df_sf))

df_sf <- st_cast(df_sf, 'POINT')

values <- raster::extract(all_rasters, df_sf)
complete <- cbind(df_sf, values)

Featured image attribution: Photo by Z on Unsplash