This tutorial/workshop assumes that you have existing vector spatial data that you want to analyze. You can create polygons and such directly in R, but we will save that for another time…

1 Reading shapefiles

Here are two different ways to read in spatial vector data, each with advantages and disadvantages:

1.0.1 reading shapefiles using rgdal::readOGR()

dir_spatial   <- 'spatial_data'
layer_bc <- 'ohibc_rgn_simple'

poly_bc_rgn <- readOGR(dsn = dir_spatial, layer = layer_bc, stringsAsFactors = FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "spatial_data", layer: "ohibc_rgn_simple"
## with 8 features
## It has 3 fields
### you can easily plot polygons, though complex polys take a while
plot(poly_bc_rgn, col = 'light blue', border = 'blue')

  • dsn is the data source name. It can be a .gdb or a directory with shapefiles in it. It can be relative or absolute file path.
    • dsn is not a fan of ~ as a shortcut for home directory, e.g. dsn = '~/github/ohiprep.
    • To get around this, one option is to use path.expand(), e.g. dsn = path.expand('~/github/ohiprep')
  • layer is the layer name.
    • If you’re reading from a .gdb, you can tell it which layer to pull out of that .gdb
      • ogrListLayers() can help identify the layers within the .gdb if you’re not sure.
    • If you’re reading a .shp (and the associated .dbf etc), then layer should be the base name of the shapefile (without the .shp extension). E.g. if you want rgn_all_mol_med_res.shp, then layer = 'rgn_all_mol_low_res'.
  • p4s allows you to input a proj.4 string to indicate projection information (CRS - coordinate reference system). If there’s already a .prj file associated with this layer, readOGR() will automatically read it in.

Quick aside on projections and coordinate reference systems:

  • WGS 84: lat and long in degrees. A basic projection, very common to see esp for global data. Not great for rasters in spatial analysis, since each cell will contain a different amount of area depending on your latitude.
    • EPSG:4326 WGS 84
    • proj.4 string: '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'
  • Mollweide: equal area projection, units in meters. Works well for global datasets; looks like a map though the edges get distorted. Equal area projections work well in conjunction with rasters - each cell will contain the same amount of area.
    • EPSG:54009
    • proj.4 string: '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
  • BC Albers: equal area projection (so it’s good for rasters), units in meters. Albers projections are often used for smaller regions, each with its own version. British Columbia’s official projection is BC Albers, which is what we will be using today (mostly).
    • EPSG:3005 - NAD83 / BC Albers
    • proj.4 string: '+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'

1.0.2 reading shapefiles using maptools::readShapePoly()

shp_bc <- file.path(dir_spatial, layer_bc)
p4s_bc <- CRS('+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0')

poly_bc_rgn <- readShapePoly(fn = shp_bc, proj4string = p4s_bc)

plot(poly_bc_rgn, col = 'light blue', border = 'blue')