Map of study area

Load Libraries

In [1]:

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggspatial)
library(sf)

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

library(rnaturalearth)
library(patchwork)
library(OpenStreetMap)

Important coordinants

Regional site location

In [2]:

site <- st_sfc(st_point(c(-99.924176, 50.056015)), crs = 4326)
box = c(xmin = -99.99, xmax = -99.924176 + 0.12, ymax = 50.056015 +.065, ymin = 50.056015 - 0.065)

Study plot locations

In [3]:

# sites <- st_as_sf(data.frame(site = c(1,2,3,4), lat = c(50.052525, 50.052209, 50.059208, 50.060235), long = c(-99.924242, -99.918429, -99.912849, -99.931591)), coords = c("long", "lat"), crs = 4362)
                                                                                              
sites <- data.frame(site = c(1,2,3,4), lat = c(50.052525, 50.052209, 50.059208, 50.060235), long = c(-99.924242, -99.918429, -99.912849, -99.931591))

Canada

In [4]:

canada <- ne_states(country = "Canada", returnclass = "sf") %>%
  st_transform(crs = 3348) 

Maps

Canada

In [5]:

p1 <- ggplot() +
  theme_bw() +
  layer_spatial(canada, fill = "white") +
  layer_spatial(site, shape = 23, colour = "red", fill = "red", size = 2) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #panel.border = element_blank(),
        panel.background = element_blank())

In [6]:

sf_use_s2(FALSE)

Spherical geometry (s2) switched off

mb_river <- ne_download(scale = 10,  type = 'rivers_lake_centerlines', category = "physical", returnclass = "sf") %>%
  st_transform(crs = 4326) %>%
  st_crop(st_bbox(c(xmin = -101.5889, ymin = 48.99267, 
                    xmax = -95.15399, ymax = 51.399), 
                  crs = st_crs(.)))

Reading layer `ne_10m_rivers_lake_centerlines' from data source 
  `/tmp/RtmpRXaZVX/ne_10m_rivers_lake_centerlines.shp' using driver `ESRI Shapefile'
Simple feature collection with 1473 features and 38 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -164.9035 ymin: -52.15775 xmax: 177.5204 ymax: 75.79348
Geodetic CRS:  WGS 84

although coordinates are longitude/latitude, st_intersection assumes that they
are planar

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

mb_lakes <- ne_download(scale = 10,  type = 'lakes', category = "physical", returnclass = "sf") %>%
  st_transform(crs = 4326) %>%
  st_crop(st_bbox(c(xmin = -101.5889, ymin = 48.99267, 
                    xmax = -95.15399, ymax = 51.399), 
                  crs = st_crs(.)))

Reading layer `ne_10m_lakes' from data source `/tmp/RtmpRXaZVX/ne_10m_lakes.shp' using driver `ESRI Shapefile'
Simple feature collection with 1355 features and 41 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -165.9656 ymin: -50.66967 xmax: 177.1544 ymax: 81.95521
Geodetic CRS:  WGS 84

although coordinates are longitude/latitude, st_intersection assumes that they
are planar

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

mb <- ne_states(country = "Canada", returnclass = "sf") %>%
  filter(name == "Manitoba") %>%
  st_transform(crs = 4326) %>%
  st_crop(st_bbox(c(xmin = -101.5889, ymin = 48.99267, 
                    xmax = -95.15399, ymax = 51.399), 
                  crs = st_crs(.)))

although coordinates are longitude/latitude, st_intersection assumes that they
are planar

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

mb_cities <- maps::canada.cities %>%
  filter(country.etc == "MB") %>%
  mutate(name = str_remove(name, " [A-Z]{2}$")) %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>%
  st_transform(crs = 4326) %>%
  st_crop(st_bbox(c(xmin = -101.5889, ymin = 48.99267, 
                    xmax = -95.15399, ymax = 51.399), 
                  crs = st_crs(.))) %>%
  filter(pop >5000)

although coordinates are longitude/latitude, st_intersection assumes that they
are planar

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

p4 <- ggplot() +
  theme_bw() +
  geom_sf(data = mb, fill = "white") +
  geom_sf(data = mb_river, colour = "skyblue4") +
    geom_sf(data = mb_lakes, fill = "lightblue") +
  geom_sf(data = mb_cities) +
  ggrepel::geom_label_repel(
    data = mb_cities,
    aes(label = name, geometry = geometry),
    stat = "sf_coordinates",
    min.segment.length = 0,
    nudge_y = 0.02,
    size = 2
  ) +
  layer_spatial(site, shape = 23, colour = "red", fill = "red", size = 2) +
  labs(tag = "a)") +
  annotation_scale(location = "bl",
                   height = unit(0.05, "cm"), 
                   pad_y = unit(0.5, "cm"),
                   pad_x = unit(1.1, "cm")) +
  annotation_north_arrow(location = "br", 
                         height = unit(0.5, "cm"),
                         width = unit(0.5, "cm"), 
                   pad_y = unit(0.5, "cm"),
                   pad_x = unit(1.1, "cm")) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        plot.margin = unit(c(0,0,0,0), "mm"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())



map <- openmap(c(50.047, -99.936), c( 50.065, -99.909), zoom = NULL,
                  type = "esri-imagery", mergeTiles = TRUE) 

sa_map2 <- openproj(map, projection = "+proj=longlat +datum=WGS84 +no_defs +type=crs")

sites <- data.frame(sites = c("1", "2", "3", "4"), lat = c(50.052472, 50.052288, 50.059197, 50.060114), long = c(-99.924327, -99.918839, -99.912919, -99.931479)) 


p5 <- OpenStreetMap::autoplot.OpenStreetMap(sa_map2) +
  theme_void() +
  geom_point(data = sites,
           aes(x = long , y = lat ), 
           colour = "red", size =  2.5) +
  ggrepel::geom_label_repel(data = sites, 
            aes(long, lat, label = sites), 
            size = 3, colour = "black",
            min.segment.length = 0) + 
  geom_segment(aes(x = -99.934824, y = 50.048, xend = -99.927655, yend = 50.048), colour = "black", linewidth = 1) +
  geom_segment(aes(x = -99.931150, y = 50.048, xend = -99.927655, yend = 50.048), colour = "white", linewidth = 1) +
  annotate("label", x = -99.926, y = 50.048, label = "500m", fill = "lightgray") +
  labs(tag = "b)")

In [7]:

p6 <- p4 + inset_element(p1, left = 0.6, bottom = 0.6, right = 1, top = 1) + p5
p6

Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data

#ggsave(plot = p6, filename = here::here("./Figures/Figure 1.tiff"), width = 174, height = 100, units = "mm", dpi = 400) 

Showing a) the location of the study site in southern Manitoba with an inset map of Canada; and b) the locations of the four riparian areas included in this study