Working with bcmaps layers with point location data

We often want to be able to add point data to a map. This vignette will run through a simple example of converting a data.frame of latitude/longitude pairs into a SpatialPointsDataFrame, and reprojecting it so that the points can be overlaid on a layer from the bcmaps package.

For this example, say we have done surveys for a species across B.C., and we want to be able to add the presences and absences on to a map of British Columbia.

library(sf)
library(bcmaps)

We will create a mock dataframe of locations of species presence/absences (in real life these would probably be in a csv or Excel file that we would import):

set.seed(42)
spp <- data.frame(site_num = LETTERS[1:10], spp_present = sample(c("yes", "no"), 10, replace = TRUE),
                 lat = runif(10, 49, 60), long = runif(10, -128, -120), 
                 stringsAsFactors = FALSE)
head(spp)

Next we convert this to a SpatialPointsDataFrame using the sp package. The syntax for doing this is a little obscure - you specify the columns that hold the x and y data in a character vector on the right hand side of the <-:

spp <- st_as_sf(spp, coords = c("long", "lat"))
summary(spp)

plot(spp["spp_present"])

We know that our occurrence data are in decimal degrees in NAD83, so we will assign the corresponding proj4string. You can specify the projection using a full proj4 string ("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs") or the EPSG code (4269). We will use the EPSG code because it’s shorter and less likely to make a typo with.

spp <- st_set_crs(spp, 4326)

All of the maps in the bcmaps package are in BC Albers projection (EPSG:3005), which is the B.C. government standard. It is an ‘equal area’ projection, meaning that the size of areas is not distorted, and thus is suitable for analyses on large areas.

If we look at the proj4string for bc_bound and our spp_df, we see that they are different:

bc_bound <- get_layer("bc_bound")
st_crs(bc_bound)
## Coordinate Reference System:
##   User input: EPSG:3005 
##   wkt:
## PROJCS["NAD83 / BC Albers",
##     GEOGCS["NAD83",
##         DATUM["North_American_Datum_1983",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6269"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4269"]],
##     PROJECTION["Albers_Conic_Equal_Area"],
##     PARAMETER["standard_parallel_1",50],
##     PARAMETER["standard_parallel_2",58.5],
##     PARAMETER["latitude_of_center",45],
##     PARAMETER["longitude_of_center",-126],
##     PARAMETER["false_easting",1000000],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["Easting",EAST],
##     AXIS["Northing",NORTH],
##     AUTHORITY["EPSG","3005"]]
st_crs(spp)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

So let’s transform the dataframe of species presence/absences into the same CRS as bc_bound:

spp <- transform_bc_albers(spp)

plot(spp["spp_present"], expandBB = rep(0.2, 4), graticule = TRUE)
plot(st_geometry(bc_bound), add = TRUE)

Now we want to know what ecoregion of the province each of these observations was in. We can use the ecoregions data from bcmaps, and the st_join function from the sf package to extract ecoregions from the point data and add that information:

Working with bcmaps layers with point location data

2020-04-28