index.qmd

---
title: "ZULE SpatialR Workshop"
author: Isabella C Richmond
format: 
  revealjs:
    theme: simple
    slide-number: c
editor: visual
---

## 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

::: {layout="[[-1], [1], [-1]]"}
::: columns
::: {.column width="50%"}
-   Reproducibility/documentation
-   Continuity
-   Visualization
-   Resources (see the end of this presentation)
-   All the other usual stuff
:::

::: {.column width="50%"}
![](presentation_imgs/esri.jpeg)
:::
:::
:::

## 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`](https://r-spatial.github.io/sf/)
2.  Raster data -\> [`stars`](https://r-spatial.github.io/stars/) or [`terra`](https://rspatial.org/terra/pkg/1-introduction.html)
3.  Interactive visualization -\> [`mapview`](https://r-spatial.github.io/mapview/)
4.  Basemaps -\> [`osmdata`](https://docs.ropensci.org/osmdata/), [`basemapR`](https://github.com/Chrisjb/basemapR), [`ggmaps`](https://github.com/dkahle/ggmap)
5.  Static visualization -\> [`ggplot2`](https://ggplot2.tidyverse.org/)

## Spatial Data Types

![](presentation_imgs/rastervector.jpg)

## Projections

::: columns
::: {.column width="50%"}
-   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](https://epsg.io/)
:::

::: {.column width="50%"}
::: {layout="[[-1], [1], [-1]]"}
![](presentation_imgs/PROJ.jpg)
:::
:::
:::

## General Workflow {.scrollable}

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

```{r}
#| echo: true
#| eval: false
#| label: install

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

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

```

## Step 2/3: Identify type of data & Import

::: panel-tabset
## Vector: `sf`

```{r}
#| echo: true
#| eval: false
#| label: read-sf

# 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")

```

## Raster: `stars`

```{r}
#| echo: true
#| eval: false
#| label: read-stars

# 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

```{r}
#| echo: true
#| eval: false
#| label: proj

# 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

::: panel-tabset
## Vector: `sf`

```{r}
#| echo: true
#| eval: false
#| label: calc-sf

# 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()
```

## Raster: `stars`

```{r}
#| echo: true
#| eval: false
#| label: calc-stars

# 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

::: panel-tabset
## `basemaps`

```{r}
#| echo: true
#| eval: false
#| label: basemaps

# 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")

```

## `osmdata`

```{r}
#| echo: true
#| eval: false
#| label: osmdata

# 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!

::: panel-tabset
## Vector
```{r}
#| echo: true
#| eval: false
#| label: plot-vector

# 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
)

```



## Raster

## Vector + Raster
:::

## 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 {.scrollable}

-   [Coordinate systems and projections](https://pro.arcgis.com/en/pro-app/latest/help/mapping/properties/coordinate-systems-and-projections.htm)
-   [sf](https://r-spatial.github.io/sf/)
-   [stars](https://r-spatial.github.io/stars/articles/stars1.html)
-   [terra](https://rspatial.org/terra/)
-   [osmdata](https://docs.ropensci.org/osmdata/)
-   [ggplot2](https://ggplot2.tidyverse.org/)
-   [Intro to sf and stars](https://keen-swartz-3146c4.netlify.app/sf.html)
-   [terra vs stars](https://www.r-bloggers.com/2021/05/a-comparison-of-terra-and-stars-packages/)
-   [osmdata feature keys](https://wiki.openstreetmap.org/wiki/Map_features)
-   [terra methods](https://rspatial.github.io/terra/reference/terra-package.html)