ZULE SpatialR Workshop

Isabella C Richmond

Spatial Analysis Software

Our lab’s expertise lies in MANY softwares (which is extremely cool), they include:

  • ArcGIS (proprietary/click & point)
  • ENVI (proprietary)
  • QGIS (click & point)
  • Google Earth Engine
  • R

Very Briefly - Some Advantages of R

 

  • Reproducibility/documentation
  • Continuity
  • Visualization
  • Resources (see the end of this presentation)
  • All the other usual stuff

Package network

LOTS of packages for spatial analysis in R - can get really specific for certain analysis types/goals. The packages below will allow you do 95% of the spatial analyses and visualizations that you want to do in R:

  1. Vector data -> sf
  2. Raster data -> stars or terra
  3. Interactive visualization -> mapview
  4. Basemaps -> osmdata, basemapR, ggmaps
  5. Static visualization -> ggplot2

Spatial Data Types

Projections

  • Most common projection is WGS 84 (EPSG: 4326)
  • Often a projection that uses meters as units is useful, for us NAD83/MTM Zone 8 can be useful (EPSG: 32188)
  • Projection codes can be found at epsg.io

 

General Workflow

  1. Install & load packages

  2. Identify type of spatial data (vector vs raster)

  3. Import data (csv, shapefile, geopackage, tif) (+ visualize)

  4. Assign projection and transform data to match (+ visualize)

  5. Do any necessary spatial calculations (+ visualize)

    • e.g., buffers, intersections, etc.

  6. Prepare for mapping

    • bounding boxes, basemaps, etc.

  7. Map

Step 1: Install and load packages

install.packages(c('sf', 'stars', 'mapview', 'ggplot2', 
                   'dplyr', 'osmdata', 'basemaps'))

source('scripts/0-packages.R')

Step 2/3: Identify type of data & Import

# let's load in the Ruelles Vertes data 
# data from https://donnees.montreal.ca/ville-de-montreal/ruelles-vertes
# shapefile containing polygons for 935 ruelles vertes across Montreal

rv <- read_sf("input/ruelles-vertes.shp")
# okay, now let's load in some canopy cover data 
# tif file containing pixels of land cover type across Montreal
cc <- read_stars("input/canopy-cover.tif")

Step 4: Projections

# let's check the projections of both of our datasets 
rv

cc 

# unfortunately, they don't match (they rarely do!)
# the projection of canopy is in metres, so we will project the ruelles to match it 
rv_t <- st_transform(rv, st_crs(cc))

# since the crs is NAD83 MTM Zone 8, we could also use the ESPG code to reproject
rv_t_code <- st_transform(rv, 32188)

Step 5: Spatial Calculations

# merge polygons with the same IDs
rv_m <- rv_t |>
  group_by(RUELLE_ID) |>
  summarise(geometry = st_union(geometry), 
                   nhood = first(PROPRIETAI),
                   codenhood = first(CODE_ARR),
                   date = first(DATE_AMENA))

# calculate area of our sampling sites 
rv_a <- rv_m |> 
  mutate(area = st_area(geometry))
# when we inspect what kind of data this is, it is a vector of class "units"
class(rv_a$area)
# if we want to use this info in future operations - we need to convert from units to double
rv_a <- rv_a |> 
  mutate(area = as.double(area))

# produce buffers surrounding sites of interest
buffers <- st_buffer(rv_a, 50) #NOTE: units of projection are in meters

# sample 3 random sampling points within each ruelle 
samp <- st_sample(rv_a, c(3,3), type = "random", by_polygon = F) |>
  st_as_sf()
# let's intersect our land use/canopy cover with the buffers we created 
cc_int <- cc[buffers]


# for each ruelle's buffer, let's calculate % canopy cover
# calculate canopy cover by dividing number of pixels == 4 (canopy) by total number of pixels within a buffer
buff_can <- aggregate(cc, buffers, FUN = function(x) sum(x == 4)/length(x)) |>
    st_as_sf() |>
    rename(percan_2021 = `canopy-cover.tif`) |>
    mutate(percan_2021 = round(percan_2021, 2))

Step 6: Mapping Prep

# set up bounding box - order: xmin, ymin, xmax, ymax
bb <- c(xmin = -74.0788,
        ymin = 45.3414,
        xmax = -73.3894,
        ymax = 45.7224)

# or use layer to create bounding box 
# this creates a bounding box around all the buffers we created
bb_layer <- st_bbox(buffers)

# view all available maps 
get_maptypes()

# you can add the basemap as a layer
bbox <- st_bbox(bb)
st_crs(bbox) = 4326

ggplot() + 
  basemap_gglayer(bbox, map_service = "carto", map_type = "light") + 
  scale_fill_identity() + 
  coord_sf()

# you can also make it as a ggplot
basemap_ggplot(bbox, map_service = "carto", map_type = "dark")

# and there are different spatial classes you can use, such as stars where the basemap is returned as a stars object
basemap_stars(bbox, map_service = "osm", map_type = "streets")
# set up bounding box - order: xmin, ymin, xmax, ymax
bb <- c(xmin = -74.0788,
        ymin = 45.3414,
        xmax = -73.3894,
        ymax = 45.7224)

# Island boundary using osmdata
mtl <- opq(bb) |> # this make a call for any data found within our coordinates
  add_osm_feature(key = 'place', value = 'island') |> # we select any features that are places + islands
  osmdata_sf() # transform it into sf object
# we are returned a list of different types of geometries that match our request
# we will select the ones we are interested in and combine
multipolys <- st_make_valid(mtl$osm_multipolygons) # grab multipolygons (large islands)
polys <- st_make_valid(mtl$osm_polygons) # grab polygons (small islands)
polys <- st_cast(polys, "MULTIPOLYGON")
allpolys <- st_as_sf(st_union(polys, multipolys)) # combine geometries and cast as sf
st_crs(allpolys) = 4326 # set CRS as EPSG 4326 

# Water
# going to do the same thing as above but we want water features within our coordinates
water <- opq(bb) |>
  add_osm_feature(key = 'natural', value = 'water') |>
  osmdata_sf()
wmpols <- water$osm_multipolygons
wmpols <- st_cast(wmpols, "MULTIPOLYGON")
wmpols <- st_as_sf(st_make_valid(wmpols))
st_crs(wmpols) = 4326


# osmdata is super flexible -- you can query just about anything!

Step 7: Map!

# Theme -----------------------------------------------------
# Palette set-up in separate script for ease if we want same colours across multiple maps 
source('scripts/0-palette.R')

# define theme for ggplot - can do this in the ggplot script as well if desired
th <- theme(panel.border = element_rect(linewidth = 1, fill = NA),
            panel.background = element_rect(fill = canadacol),
            panel.grid = element_line(color = gridcol2, linewidth = 0.2),
            axis.text = element_text(size = 11, color = 'black'),
            axis.title = element_blank())

# Plot -------------------------------------------------------
# transform all layers so they have the same CRS 
# use EPSG 3347 - same projection that Statistics Canada uses
mtlcrs <- 3347
rv_m <- st_transform(rv, mtlcrs)
allpolys_m <- st_transform(allpolys, mtlcrs)
wmpols_m <- st_transform(wmpols, mtlcrs)

# want the bounds of our map to be slightly smaller than the entire island
bbi <- st_bbox(st_buffer(allpolys_m, 0.75))


# plot all layers together with theme we set above 
# if we were plotting the raster we could use the geom_stars function
plot <- ggplot() +
  geom_sf(fill = montrealcol, data = allpolys_m) + 
  geom_sf(aes(fill = percan), data = rv_m) + 
  geom_sf(fill = watercol, col = "#5b5b5b", data = wmpols_m) + 
  scale_fill_viridis_c(direction = -1) + 
  coord_sf(xlim = c(bbi['xmin'], bbi['xmax']),
           ylim = c(bbi['ymin'], bbi['ymax'])) +
  th


# Save -------------------------------------------------------
ggsave(
  'graphics/ruelles-vertes.png',
  plot,
  width = 10,
  height = 10,
  dpi = 320
)

Some Troubleshooting Tips

  • Do data exploration (always) but especially spatial data every step of the way (mapview!!)
  • Always make sure you have appropriate, matching projections
  • Inspect your attributes - including geometry type
  • If you use GitHub, pay attention to file size

Resources