# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = FALSE,
# Disable warnings printed by R code chunks
warning = FALSE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
fig_count <<- fig_count + 1
paste0("Figure ", fig_count, ": ", caption)
}library(sf)
library(tmap)
library(mapview)
# set mapview options
mapviewOptions(basemaps = c("OpenStreetMap",
"Esri.WorldImagery",
"OpenTopoMap"))
source("jobbers/custom_bounding_box.R")Let us first show the area of interest. I show it at the four available levels: country, province, canton, and parish. The data is from the GADM project, which provides high-resolution maps of administrative areas worldwide.
Here are the specific steps:
gadm41_ECU.gpkgst_layers("gadm41_ECU.gpkg"), you will see
the available layers, which are ADM_ADM_0,
ADM_ADM_1, ADM_ADM_2, and
ADM_ADM_3. These correspond to country, province, canton,
and parish, respectively.Here are the steps I follow to obtain the files below (see this script for more details):
gadm41_ECU.gpkg fileI did this because at the moment I wanted to do the analysis on the entire mainland of Ecuador.Note: The area of interest is delimited by the coordinates defined in
this file.
The corresponding sf object is stored here.
However, a more adequate area for future analysis can be found here.
This area is not a box but a polygon that follows the administrative
boundaries of the country.
st_intersection() function from the sf
package.# read the files
gadm41_ECU_country <- readRDS("clean_data/gadm41_ECU_country_custom.RDS")
gadm41_ECU_provinces <- readRDS("clean_data/gadm41_ECU_province_custom.RDS")
gadm41_ECU_cantons <- readRDS("clean_data/gadm41_ECU_canton_custom.RDS")
gadm41_ECU_parishes <- readRDS("clean_data/gadm41_ECU_parish_custom.RDS")I choose the following polygon to delimit the area of interest. The polygon is defined by the coordinates of the four corners: Northeast, Southeast, Southwest, and Northwest.
coords <- rbind(NE, SE, SW, NW, NE) # from source("jobbers/custom_bounding_box.R")
crop_polygon <- st_polygon(list(coords)) %>%
st_sfc(crs = st_crs(gadm41_ECU_country)) # match CRS of your data
crop_polygon_sf <- st_sf(geom = crop_polygon)
mapview::mapview(crop_polygon_sf,
alpha.regions = 0.4, # fill transparency
color = "black", # border color
alpha = 1, # border transparency
legend = FALSE)