Code chunks for the #30DayMapChallenge 2025

Author

Miguel Alvarez

Published

November 1, 2025

Here are the pieces of code used to produce the maps.

Day 1

# Load required packages
library(rgbif)
library(ggplot2)
library(sf)
library(dplyr)
library(rnaturalearth)
library(magick)

# Get world map
worldmap <- ne_countries(scale = 'medium', type = 'map_units',
    returnclass = 'sf')

# Get occurrence data for Cyperus papyrus
cypepap <- occ_search(scientificName = "Cyperus papyrus", hasCoordinate = TRUE)

cypepap_map <- st_as_sf(cypepap$data, 
  coords = c("decimalLongitude", "decimalLatitude"),
  crs = "+proj=longlat +datum=WGS84") %>%
  mutate(origin = continent == "AFRICA")

# Create a global bounding box (in lon/lat)
N <- 200
bg_coords <- cbind(c(-180, rep(180, N), rep(-180, N)),
  c(90, seq(90, -90, length.out = N), seq(-90, 90, length.out = N)))

bg <- st_as_sf(
  st_sfc(st_polygon(list(
    bg_coords
  ))),
  crs = 4326
)

# Create the map
water <- "#ccfcff"
m <- ggplot() + 
  geom_sf(data = bg, fill = water, color = NA) +
  geom_sf(data = worldmap, fill = "#b4b4b4", color = water) +
  geom_sf(data = cypepap_map, col = "red", show.legend = FALSE,
    cex = 1.2) +
  geom_sf(data = cypepap_map, aes(color = origin), show.legend = FALSE,
    cex = 0.5) +
  coord_sf(crs = st_crs("ESRI:54030")) +
  scale_color_manual(
    values = c("TRUE" = "#f5b53e", "FALSE" = "#fdfd8c")
  ) +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    panel.grid.major = element_blank(),
    plot.background = element_rect(fill = "transparent", color = NA)
  )

# Save the map as PNG file
img_name <- "01-points.png"
ggsave(img_name, plot = m, width = 8, height = 6, dpi = 300,
  bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 2

# Load required packages
library(rnaturalearth)
library(sf)
library(ggplot2)
library(dplyr)
library(units)
library(magick)

# download the South America map and crop continental Chile
south_america <- ne_countries(continent = "South America", scale = "large",
  returnclass = "sf") %>%
  st_make_valid() %>%
  st_crop(y = st_bbox(c(xmin = -77, ymin = -57, xmax = -65, ymax = -17)))

# Extract Chile and convert to UTM
chile <- south_america %>%
  st_transform(crs = "EPSG:32718") %>%
  filter(admin == "Chile")

# Extract nodes
nodes <- chile %>%
  st_cast("MULTILINESTRING") %>%
  st_cast("POINT")

# Compute the pairwise distance matrix
dist_matrix <- st_distance(nodes)

# Find the indices of the two farthest points
max_idx <- which(dist_matrix == max(dist_matrix), arr.ind = TRUE)[1, ]

max_line <- st_sfc(
  st_linestring(
    rbind(
      st_coordinates(nodes[max_idx[1], ]),
      st_coordinates(nodes[max_idx[2], ]))
    ), crs = st_crs(nodes))

# Measure length of the line
line_length <- st_length(max_line) %>%
  set_units("km")

# Calculate centroid of line
line_centroid <- st_centroid(max_line)
line_centroid <- st_sf(line_centroid)
line_centroid$distance <- paste0(round(as.numeric(line_length), 1), " km")

# Draw the map
m <- ggplot() + 
  geom_sf(data = chile, fill = "gray", color = NA) +
  geom_sf(data = max_line, color = "red") +
  geom_sf_label(data = line_centroid, aes(label = distance), color = "red",
    angle = 90) +
  theme(
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90)
  )

# Save the map as PNG file
img_name <- "02-lines.png"
ggsave(img_name, plot = m, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 3

# Load required packages
library(rnaturalearth)
library(sf)
library(ggplot2)
library(dplyr)
library(magick)

# download the South America map and crop continental Chile
sam <- ne_countries(continent = "South America", scale = "large",
  returnclass = "sf") %>%
  st_make_valid() %>%
  st_crop(y = st_bbox(c(xmin = -77, ymin = -57, xmax = -65, ymax = -41)))

# Merge neighbouring polygons
touch_matrix <- st_relate(sam, pattern = "F***T****")

# Build groups of connected features
library(igraph)
g <- graph_from_adjacency_matrix(touch_matrix)
groups <- components(g)$membership

# Merge polygons within each group
merged <- sam %>%
  mutate(group = groups) %>%
  group_by(group) %>%
  summarise(geometry = st_union(geometry))

# Split into single polygons and transform
merged <- merged %>%
  st_transform(crs = "EPSG:32718") %>%
  st_cast(to = "POLYGON")

merged <- merged %>%
  mutate(area = st_area(merged))

# Filter islands and classify them by their area
islands <- merged %>%
  filter(area < max(area)) %>%
  mutate(quartile = as.factor(5 - ntile(area, 4)))

# Crop a piece of continent around the islands
bbox <- st_bbox(islands)
bbox["ymax"] <- bbox["ymax"] + 100000
bbox["xmax"] <- bbox["xmax"] + 100000
bbox["xmin"] <- bbox["xmin"] - 100000

continent <- merged %>%
  filter(area == max(area)) %>%
  st_crop(y = bbox)

# Draw the graphic
m <- ggplot() +
  geom_sf(data = continent, fill = "#78d7e4ff", color = NA, alpha = 0.5) +
  geom_sf(data = islands, aes(fill = quartile), color = NA) +
  scale_fill_manual(
    values = c(
      "1" = "#ffca57ff",
      "2" = "#d0352e",
      "3" = "#4ab85cff",
      "4" = "#2427b1ff"
    ),
    name = "Area\nQuartile"  # legend title
  ) +
  xlim(450000, 1200000.) +
  ylim(3775178.7, 5340000) +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.title = element_text(size = 9),
    panel.background = element_rect(fill = "#d7eef1ff")
  )
m

# Save the map as PNG file
img_name <- "03-polygons.png"
ggsave(img_name, plot = m, height = 9, width = 6, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 4

# Load packages
library(dplyr)
library(sf)
library(terra)
library(leaflet)
library(htmlwidgets)
library(webshot2)

# Import Data
track_points <- read.csv("track_points.csv") %>%
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>%
  st_transform(crs = "EPSG:21037")

# Get bounding box (expand M and round to 100 m)
My <- 100000 # Margin added to top and bottom
Mx <- 400000 # Margin added to left and right
Z <- 100 # Zell resolution in meters
bbox <- (round(st_bbox(track_points)/Z) * Z) + c(-Mx, -My, Mx, My)

# Get raster of 100 m resolution
r <- rast(
  nrows = (bbox["ymax"] - bbox["ymin"])/Z,
  ncols = (bbox["xmax"] - bbox["xmin"])/Z,
  xmin = bbox["xmin"], xmax = bbox["xmax"],
  ymin = bbox["ymin"], ymax = bbox["ymax"],
  crs = "EPSG:21037"
)
values(r) <- 0

# pierced rasters
v_buf <- buffer(vect(track_points), width = 1000)
r_masked <- mask(r, v_buf, inverse = TRUE)

# visualize with leaflet
m <- leaflet() %>%
  addTiles() %>%
  addRasterImage(r_masked, color = "black", opacity = 0.5) %>%
  setView(lng = 36.45, lat = -0.16, zoom = 8) %>%
  addScaleBar()

# Save map
saveWidget(m, "04-my-data.html", selfcontained = TRUE)
webshot2::webshot("04-my-data.html", file = "04-my-data.png")

Day 5

# Load packages
library(elevatr)
library(sf)
library(terra)
library(ggplot2)
library(tidyterra)
library(magick)

# create an sf polygon for the bbox (CRS WGS84)
bbox <- st_as_sfc(st_bbox(c(
    xmin = 7.18,
    ymin = 50.64,
    xmax = 7.3,
    ymax = 50.70
  ), crs = st_crs(4326)))
bbox <- st_sf(geometry = bbox)

# z = zoom level / resolution control (higher = finer). Typical values 9-14.
dem_raster <- get_elev_raster(locations = bbox, z = 14, clip = "locations")

# convert to terra SpatRaster
dem_terra <- rast(dem_raster)
dem_terra

# project to UTM
dem_utm <- project(dem_terra, "EPSG:25832")

# radius from m to zells
r <- 100 # radius in m
r <- round(r/mean(res(dem_utm))) # radius converted to raster cells

# Prepare windows

# circle (From https://github.com/kamapu/spatialist)
d <- 2*r + 1 # diameter
r_ind <- matrix(rep(1:d, times = d), nrow = d) # row index
c_ind <- matrix(rep(1:d, times = d), nrow = d, byrow = TRUE) # col index

Dist <- sqrt(((r + 1) - r_ind)^2 + ((r + 1) - c_ind)^2) # Distance to Center in Zells
Circle <- Dist <= r # Circle

North <- matrix(c(rep(TRUE, times = r*d), rep(FALSE, times = d + r*d)),
    ncol = d, byrow = TRUE) & Circle
South <- matrix(c(rep(FALSE, times = d + r*d), rep(TRUE, times = r*d)),
    ncol = d, byrow = TRUE) & Circle

North[North == TRUE] <- 1
North[North == 0] <- NA
South[South == TRUE] <- 1
South[South == 0] <- NA

# Extract Elevations towards North and South using Focal
# This takes some time

North_r <- focal(dem_utm, North, fun = mean)
South_r <- focal(dem_utm, South, fun = mean)

# Northness
Northness <- South_r - North_r

# Hillshades
slope <- terrain(dem_utm, "slope", unit = "radians")
aspect <- terrain(dem_utm, "aspect", unit = "radians")
hill <- shade(slope, aspect, 40, 270)

# Stack and secure
sieben_gebirge <- rast(list(elevation = dem_utm, northness = Northness, shades = hill))
plot(sieben_gebirge)

writeRaster(sieben_gebirge, "sieben_gebirge.tif")

# Produce map with ggplot
m <- ggplot() +
  geom_spatraster(data = sieben_gebirge, aes(fill = northness)) +
  scale_fill_viridis_c(name = "Northness") +
  geom_spatraster_contour(data = sieben_gebirge, aes(z = elevation),
      breaks = c(125, 150, 175, 225, 250, 275, 325, 350, 375)) +
  geom_spatraster_contour_text(data = sieben_gebirge, aes(z = elevation),
      breaks = c(100, 200, 300, 400), color = "black") +
  theme_minimal()

# Write the image
img_name <- "05-earth.png"
ggsave(img_name, plot = m, height = 6, width = 9, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 6

# Load packages
library(dplyr)
library(sf)
library(rnaturalearth)
library(ggplot2)
library(RColorBrewer)
library(magick)

# Download South American map (without Falkland Islands and French Guyana)
sam <- ne_countries(continent = "South America", scale = "small", type = "map_units",
    returnclass = "sf") %>%
  filter(!adm0_iso %in% c("B12", "GUF"))

# Collect population data (https://www.worldometers.info/world-population/south-america-population/)
pop <- c(
  ARG = 45696159,
  CHL = 19764771,
  URY = 3386588,
  BRA = 211998573,
  BOL = 12413315,
  PER = 34217848,
  COL = 52886363,
  VEN = 28405543,
  GUY = 831087,
  SUR = 634431,
  ECU = 18135478,
  PRY = 6929153
)
pop <- data.frame(iso = names(pop), population = pop)

# Calculate centroids and side aspect
bboxes <- as.data.frame(do.call(rbind, lapply(st_geometry(sam), st_bbox))) %>%
  mutate(centr_x = (xmax + xmin)/2, centr_y = (ymax + ymin)/2) %>%
  mutate(aspect = (xmax - xmin)/(ymax - ymin)) %>%
  mutate(name = sam$admin, iso = sam$adm0_iso) %>%
  inner_join(y = pop)

# Calculating new bboxes
c <- 0.0021 # Scaling factor
bboxes_sf <- bboxes %>%
  mutate(height = sqrt(population/aspect) * c) %>%
  mutate(width = aspect * height) %>%
  mutate(
    xmin = centr_x - width/2,
    xmax = centr_x + width/2,
    ymin = centr_y - height/2,
    ymax = centr_y + height/2 
  ) %>%
  select(name, iso, population, xmin, xmax, ymin, ymax) %>%
  rowwise() %>%
  mutate(geometry = st_as_sfc(st_bbox(c(
    xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax
  ), crs = "EPSG:4326"))) %>%
  st_as_sf()

# Mapping bounding boxes
bboxes_sf$color <- brewer.pal(12, "Set3")

m <- ggplot() +
  geom_sf(data = sam, fill = NA, color = "wheat3", linewidth = 0.8) +
  geom_sf(data = bboxes_sf, aes(fill = color), color = NA, alpha = 0.7) +
  geom_sf_text(data = bboxes_sf, aes(label = iso), color = "grey20", cex = 5) +
  scale_fill_identity() +
  theme_minimal() +
  theme(
    axis.title = element_blank()
  )

# Write the image
img_name <- "06-dimensions.png"
ggsave(img_name, plot = m, height = 10, width = 8, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 8

# Load packages
library(osmdata)
library(sf)
library(lwgeom)
library(dplyr)
library(ggplot2)
library(magick)

# Do query to osm
bbox <- st_bbox(c(xmin = 7.012, ymin = 50.761, xmax = 7.13, ymax = 50.861), crs = 4326)

roads <- opq(bbox = bbox) %>%
  add_osm_feature(key = "highway") %>%
  osmdata_sf()

roads_sel <- roads$osm_lines %>%
  select(osm_id, name, highway, lanes, surface)

# Get bounding box and apply coordinates transformation
bbox <- st_as_sfc(bbox) %>%
  st_as_sf() %>%
  st_transform(crs = "EPSG:25832")

roads_sel <- roads_sel %>%
  st_transform(crs = st_crs(bbox)) %>%
  st_crop(y = bbox)

# Split the background
blocks <- st_split(bbox, roads_sel) %>%
  st_collection_extract(type = "POLYGON")

# Produce sizes and size categories
blocks$size <- as.numeric(st_area(blocks))/1000000 # Size in square kilometers
blocks$size_class <- cut(blocks$size,
  breaks = c(0, 0.01, 0.025, 0.1, max(blocks$size)), labels = FALSE)

# Some labels for References
cities <- rbind(
  Niederkassel = c(50.815267, 7.037569),
  Ranzel = c(50.832243, 7.030632),
  Rheidt = c(50.789460, 7.049923),
  Mondorf = c(50.779176, 7.066840),
  Bergheim = c(50.776256, 7.095724),
  Eschmar = c(50.792068, 7.109497),
  Kriegsdorf = c(50.808077, 7.102283),
  Uckendorf = c(50.820749, 7.064948),
  Libur = c(50.838601, 7.070700),
  Hersel = c(50.769734, 7.045189),
  Uedorf = c(50.783207, 7.032862),
  Widdig = c(50.794806, 7.024796),
  Urfeld = c(50.804954, 7.023736)
)
cities <- data.frame(name = rownames(cities), lon = cities[,2], lat = cities[,1])

# Cites to space
cities <- cities %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(crs = st_crs(blocks))

# Produce map
m <- ggplot() +
  geom_sf(data = blocks, aes(fill = as.factor(size_class)), color = NA) +
  scale_fill_manual(
    values = c(
      "4" = "#44bd4aff",
      "3" = "#c9ee8cff",
      "2" = "#e68e0bff",
      "1" = "#c53408ff"
    ),
    name = "Size\nCategories"  # legend title
  ) +
  geom_sf_label(data = cities, aes(label = name), col = "black", alpha = 0.8) +
  theme_minimal() +
  theme(
    axis.title = element_blank()
  )

# Write the image
img_name <- "08-urban.png"
ggsave(img_name, plot = m, height = 10, width = 8, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 11

# Load packages
library(sf)
library(rnaturalearth)
library(dplyr)
library(rmapshaper)
library(ggplot2)
library(magick)

# Download first-level administrative boundaries (states) for Germany
germany <- ne_states(country = "Germany", returnclass = "sf") %>%
  select(adm1_code, name, woe_label, geometry)

# Calculate Size in square kilometers
germany$size <- as.numeric(st_area(germany))/1000000
germany$main <- "all nodes"

# Simplify polygons, keep 10% of points
germany_s10 <- ms_simplify(
  germany,
  keep = 0.1,
  keep_shapes = TRUE
)
germany_s10$main <- "10% of nodes"

# Simplify polygons, keep 1% of points
germany_s01 <- ms_simplify(
  germany,
  keep = 0.01,
  keep_shapes = TRUE
)
germany_s01$main <- "1% of nodes"

# Simplify polygons, keep 0.5% of points
germany_s005 <- ms_simplify(
  germany,
  keep = 0.005,
  keep_shapes = TRUE
)
germany_s005$main <- "0.5% of nodes"

# Combine all version in one object
germany_all <- bind_rows(germany, germany_s10, germany_s01, germany_s005)
germany_all$main <- factor(germany_all$main,
  levels = c("all nodes", "10% of nodes", "1% of nodes", "0.5% of nodes"))

# Do ggplot
m <- ggplot(germany_all) +
  geom_sf(aes(fill = size), color = "black") +
  facet_wrap(~ main, nrow = 1) +  # all plots in a row
  scale_fill_viridis_c(name = "State\nArea" ~ (km^2), direction = -1) + # single legend
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

# Write the image
img_name <- "11-minimal-map.png"
ggsave(img_name, plot = m, height = 8, width = 15, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 14

# Load packages
library(osmdata)
library(sf)
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(magick)

# Do query to osm
bbox <- st_bbox(c(
    xmin = 7.081,
    ymin = 50.771,
    xmax = 7.109,
    ymax = 50.785),
  crs = 4326)

roads <- opq(bbox = bbox) %>%
  add_osm_feature(key = "highway") %>%
  osmdata_sf()

# Extract lines and reclassify
roads_mod <- roads$osm_lines %>%
  select(osm_id, name, highway, lanes, surface) %>%
  mutate(highway = recode(highway,
    "living_street" = "residential",
    "service" = "residential",
    "footway" = "path",
    "track" = "path",
    "unclassified" = "path",
    "tertiary" = "secondary"
  )) %>%
  filter(!highway %in% c("platform", "proposed", "steps"))

# Visualize
ggplot() +
  geom_sf(data = roads_mod, aes(color = highway)) +
  facet_wrap(~highway)

# Get bounding box and apply coordinates transformation
bbox <- st_as_sfc(bbox) %>%
  st_as_sf() %>%
  st_transform(crs = "EPSG:25832")

roads_mod <- roads_mod %>%
  st_transform(crs = st_crs(bbox)) %>%
  st_crop(y = bbox)

# Create buffers
buffers <- roads_mod %>%
  mutate(buffer_width = case_when(
    highway == "path" ~ 4,
    highway == "residential" ~ 6,
    highway == "secondary" ~ 10
  )) %>%
  mutate(geometry = st_buffer(geometry, dist = buffer_width)) %>%
  st_union()

# Split background
bg <- st_difference(bbox, buffers) %>%
  st_cast("POLYGON")

# Select colors
palette <- brewer.pal(8, "Accent")

set.seed(42)
bg <- bg %>%
  mutate(rcols = sample(palette, nrow(bg), replace = TRUE))

m <- ggplot(data = bg) +
  geom_sf(aes(fill = rcols), color = NA, alpha = 0.7) +
  geom_sf(fill = NA, color = "#927a4fff", linewidth = 0.8) +
  theme_void() +
  theme(legend.position = "none")
m

# Write the image
img_name <- "14-openstreetmap.png"
ggsave(img_name, plot = m, height = 7, width = 9, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 16

# Load packages
library(elevatr)
library(sf)
library(terra)
library(isocubes)
library(grid)
library(magick)

# create an sf polygon for the bbox (CRS WGS84)
bbox <- st_as_sfc(st_bbox(c(
    xmin = 7.18,
    ymin = 50.64,
    xmax = 7.3,
    ymax = 50.70
  ), crs = st_crs(4326)))
bbox <- st_sf(geometry = bbox)

# z = zoom level / resolution control (higher = finer). Typical values 9-14.
dem_raster <- get_elev_raster(locations = bbox, z = 14, clip = "locations")

# convert to terra SpatRaster
dem_terra <- rast(dem_raster)
dem_terra

# project to UTM
dem_utm <- project(dem_terra, "EPSG:25832")

# resample to a grid
N <- 150   # Horizontal cell size in m

template <- rast(ext(dem_utm), 
                 resolution = N, 
                 crs = crs(dem_utm))

# Resample DEM to 100 m grid (bilinear or nearest)
dem100 <- resample(dem_utm, template, method = "bilinear")

# Reclssify into 100 m elevation classes
N <- 50   # Vertical cell size in m
max_elev <- ceiling(max(values(dem100), na.rm = TRUE) / N) * N

class_matrix <- cbind(
  from = seq(0, max_elev - N, N),
  to   = seq(N, max_elev, N),
  becomes = seq(0, max_elev, N) / N
)

dem_class <- classify(dem100, class_matrix)

# Convert raster to data.frame
df <- as.data.frame(dem_class, xy = TRUE)
colnames(df) <- c("xcoord", "ycoord", "class")

# Remove NA
df <- df[!is.na(df$class), ]

# Convert geographic coordinates to grid indices
# (isocubes uses integer grid positions)
df$col <- terra::colFromX(dem_class, df$xcoord)
df$row <- terra::rowFromY(dem_class, df$ycoord)

# Build cube coordinates
coords <- do.call(rbind, lapply(1:nrow(df), function(i) {
  if (df$class[i] == 0) return(NULL)
  z_vals <- 0:(df$class[i] - 1)
  cbind(
    x = df$col[i],
    y = df$row[i],
    z = z_vals
  )
}))

coords <- as.data.frame(coords)

# Color palette
max_z <- max(coords$z)
pal <- terrain.colors(max_z + 1)

coords$fill <- pal[coords$z + 1]

# visualize with isocubes
g <- isocubesGrob(
  coords = coords,
  fill = coords$fill,   # outlines only
  col  = "black",
  size = 2.5,
  y = 0.25,
  angle = 30
)

# Open PNG device
img_name <- "16-cell.png"

png(img_name,
    width = 3000,   # adjust
    height = 3000,  # adjust
    res = 300,
    bg = NA)      # high-res print quality

grid.newpage()
grid.draw(g)

dev.off()

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 19

# Load required packages
library(sf)
library(lwgeom)
library(cowplot)
library(dplyr)
library(ggplot2)
library(patchwork)
library(rnaturalearth)
library(magick)

# An overview of projections here:
# https://map-projections.net/imglist.php

# Get world map
worldmap <- ne_countries(scale = 'medium', type = 'map_units',
    returnclass = 'sf')

# Get places to connect
places <- rbind(
    "Puerto Guadal" = c(-46.84587, -72.70271),
    Valdivia = c(-39.82040, -73.23898),
    Llaver<c3><ad>a = c(-34.17339, -71.32848),
    Bonn = c(50.73616, 7.10202),
    "Addis Abeba" = c(8.9919, 38.7622),
    Marigat = c(0.47189, 35.98589),
    Korogwe = c(-5.15626, 38.45412)
)

places <- data.frame(name = rownames(places), lon = places[ , 2],
    lat = places[ , 1]) %>%
    st_as_sf(coords = c("lon", "lat"), crs = 4326)

# Produce a connector
# For segmentation, see this discussion:
# https://github.com/r-spatial/sf/issues/1078

D <- 1 # segment length in degrees

link <- places %>%
    st_coordinates() %>%
    st_linestring() %>%
    st_sfc(crs = st_crs(places)) %>%
    st_as_sf() %>%
    st_set_crs(NA) %>%
    st_segmentize(dfMaxLength = D) %>%
    st_set_crs(st_crs(places))

# Produce the background
bg <- list(cbind(x = c(rep(c(-180, 180), each = 2), -180),
    y = c(-90, rep(90, 2), rep(-90, 2)))) %>%
    st_polygon() %>%
    st_sfc(crs = st_crs(places)) %>%
    st_as_sf() %>%
    st_set_crs(NA) %>%
    st_segmentize(dfMaxLength = D) %>%
    st_set_crs(st_crs(places))

# Draw projections -------------------------------------------------------------

# plain carre
water <- "#ccfcff"

p0 <- ggplot() + 
  geom_sf(data = bg, fill = water, color = NA) +
  geom_sf(data = worldmap, fill = "#b4b4b4", color = water) +
  geom_sf(data = link, color = "darkgreen") +
  geom_sf(data = places, color = "darkgreen", cex = 1.8) +
  geom_sf(data = places, color = "yellow", cex = 1) +
  theme_map() +
  theme(
    plot.title = element_text(hjust = 0.5)  # Center the title
  )

p1 <- p0 +
    ggtitle("Plate carr<c3><a9>e")
p1

# Robinson
p2 <- p0 +
    coord_sf(crs = st_crs("ESRI:54030")) +
    ggtitle("Robinson")
p2

# Winkel tripel
# See https://wilkelab.org/practicalgg/articles/Winkel_tripel.html

wintri <- list(crs = "+proj=wintri +datum=WGS84 +no_defs +over")

for(i in c("worldmap", "bg", "link", "places"))
    wintri[[i]] <- st_transform_proj(get(i), crs = wintri$crs)

p3 <- ggplot() +
    geom_sf(data = wintri$bg, fill = water, color = NA) +
    geom_sf(data = wintri$worldmap, fill = "#b4b4b4", color = water) +
    geom_sf(data = wintri$link, color = "darkgreen") +
    geom_sf(data = wintri$places, color = "darkgreen", cex = 1.8) +
    geom_sf(data = wintri$places, color = "yellow", cex = 1) +
    coord_sf(datum = NULL) +
    theme_map() +
    theme(
        plot.title = element_text(hjust = 0.5)  # Center the title
    ) +
    ggtitle("Winkel tripel")
p3

# Mollweide
p4 <- p0 +
    coord_sf(crs = "+proj=moll") +
    ggtitle("Mollweide")
p4

# Bonne
p5 <- p0 +
    coord_sf(crs = st_crs("ESRI:54024")) +
    ggtitle("Bonne")
p5

# Orthographic
grid <- st_make_grid(bg, cellsize = c(0.5, 1), what = "polygons")

p6 <- ggplot() + 
  geom_sf(data = grid, fill = water, color = water) +
  geom_sf(data = worldmap, fill = "#b4b4b4", color = water) +
  geom_sf(data = link, color = "darkgreen") +
  geom_sf(data = places, color = "darkgreen", cex = 1.8) +
  geom_sf(data = places, color = "yellow", cex = 1) +
  theme_map() +
  theme(
    plot.title = element_text(hjust = 0.5),  # Center the title
  ) +
  coord_sf(crs = "+proj=ortho") +
  ggtitle("Orthographic")
p6

# All plots together
all_in_one <- (p1 | p2) /
              (p3 | p4) /
              (p5 | p6)

all_in_one

# Save the map as PNG file
img_name <- "19-projections.png"
ggsave(img_name, plot = all_in_one, width = 8, height = 12, dpi = 300,
  bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 22

# Load packages
library(dplyr)
library(sf)
library(terra)
library(rnaturalearth)
library(ggplot2)
library(tidyterra)
library(cowplot)
library(RColorBrewer)
library(magick)

# Download South American map (without Falkland Islands and French Guyana)
sam <- ne_countries(continent = "South America", scale = "small", type = "map_units",
    returnclass = "sf") %>%
  filter(!adm0_iso %in% c("B12", "GUF"))

# Download capital cities
sam_caps <- ne_download(type = "populated_places", category = "cultural",
  returnclass = "sf") %>%
  filter(ADM0NAME %in% sam$admin &
      FEATURECLA == "Admin-0 capital")

# Create a raster as template
empty_rast <- rast(ext(sam), resolution = 0.1)
crs(empty_rast) <- st_crs(sam)$wkt

# Rasterized countries
sam_r <- rasterize(vect(sam), empty_rast, field = "adm0_iso")

# Empty raster to collect results
dist_r <- sam_r
values(dist_r) <- NA

for (cid in unique(sam$adm0_iso)) {
  message("Processing: ", cid)
  # subset country polygon and its capital
  r <- mask(sam_r, sam_r == cid, inverse = TRUE, maskvalues = TRUE)
  p <- sam_caps %>%
    filter(ADM0_A3 == cid)
  # distance raster inside country
  d <- distance(r, vect(p))/1000
  d <- mask(d, r)
  # write it into the output raster
  dist_r <- cover(dist_r, d)
}

# Reclassify by a distance step
steps <- 250 # in km
maxd <- as.integer(global(dist_r, "max", na.rm = TRUE)["adm0_iso", ]/steps) + 1
breaks <- c(0:maxd)*steps
breaks <- cbind(head(breaks, -1), tail(breaks, -1), seq_along(breaks[-1]))

dist_class <- classify(
  dist_r,
  rcl = breaks,
  include.lowest = TRUE,
  right = FALSE
)

# set palette
the_colors <- colorRampPalette(brewer.pal(9, "YlOrRd"))(nrow(breaks))

# Plotting map
m <- ggplot() +
  geom_spatraster(data = as.factor(dist_class)) +
  scale_fill_manual(
    values = the_colors,
    na.value = NA,
  ) +
  geom_sf_text(data = sam_caps, aes(label = NAME_ES), angle = 45) +
  theme_map() +
  theme(legend.position = "none")
m

# Save the map as PNG file
img_name <- "22-naturalearth.png"
ggsave(img_name, plot = m, width = 8, height = 10, dpi = 300,
  bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 23

# Load packages
library(osmdata)
library(sf)
library(lwgeom)
library(dplyr)
library(ggplot2)
library(patchwork)
library(magick)

# Do query to osm
bbox <- st_bbox(c(xmin = 7.012, ymin = 50.761, xmax = 7.13, ymax = 50.861), crs = 4326)

roads <- opq(bbox = bbox) %>%
  add_osm_feature(key = "highway") %>%
  osmdata_sf()

roads <- roads$osm_lines %>%
  select(osm_id, name, highway, lanes, surface)

# Get bounding box and apply coordinates transformation
bbox <- st_as_sfc(bbox) %>%
  st_as_sf() %>%
  st_transform(crs = "EPSG:25832")

roads <- roads %>%
  st_transform(crs = st_crs(bbox)) %>%
  st_crop(y = bbox)

# Plot 1: raw data
m1 <- ggplot() +
  geom_sf(data = bbox, fill = "grey", color = NA) +
  geom_sf(data = roads, col = "black") +
  labs(title = "Road network") +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
  )
m1

# Split the background and calculate surface
blocks <- st_split(bbox, roads) %>%
  st_collection_extract(type = "POLYGON")
blocks$size <- as.numeric(st_area(blocks))/1000000 # Size in square kilometers

# Plot 2: histogram
m2 <- ggplot(data = blocks, aes(x = size)) +
  geom_histogram(bins = 200) +
  xlim(0, 0.05) +
  labs(title = "Size values distribution") +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
  )
m2

# Setting manual cut levels
blocks <- blocks %>%
  mutate(size_class = as.factor(cut(size,
      breaks = c(0, 0.01, 0.03, 0.25, max(blocks$size)),
      labels = FALSE)))

# Produce the map
m3 <- ggplot() +
  geom_sf(data = blocks, aes(fill = size_class), color = NA) +
  scale_fill_manual(
    values = c(
      "4" = "#44bd4aff",
      "3" = "#c9ee8cff",
      "2" = "#e68e0bff",
      "1" = "#c53408ff"
    ),
    name = "Size\nCategories"  # legend title
  ) +
  theme_minimal() +
  theme(
    axis.title = element_blank()
  ) +
  labs(title = "Classified blocks") +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
  )
m3

# Add cities to the map

# Some labels for References
cities <- rbind(
  Niederkassel = c(50.815267, 7.037569),
  Ranzel = c(50.832243, 7.030632),
  Rheidt = c(50.789460, 7.049923),
  Mondorf = c(50.779176, 7.066840),
  Bergheim = c(50.776256, 7.095724),
  Eschmar = c(50.792068, 7.109497),
  Kriegsdorf = c(50.808077, 7.102283),
  Uckendorf = c(50.820749, 7.064948),
  Libur = c(50.838601, 7.070700),
  Hersel = c(50.769734, 7.045189),
  Uedorf = c(50.783207, 7.032862),
  Widdig = c(50.794806, 7.024796),
  Urfeld = c(50.804954, 7.023736)
)
cities <- data.frame(name = rownames(cities), lon = cities[,2], lat = cities[,1])

# Cites to space
cities <- cities %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  st_transform(crs = st_crs(blocks))

# Produce map
m4 <- ggplot() +
  geom_sf(data = blocks, aes(fill = as.factor(size_class)), color = NA) +
  scale_fill_manual(
    values = c(
      "4" = "#44bd4aff",
      "3" = "#c9ee8cff",
      "2" = "#e68e0bff",
      "1" = "#c53408ff"
    ),
    name = "Size\nCategories"  # legend title
  ) +
  geom_sf_label(data = cities, aes(label = name), col = "black", alpha = 0.8) +
  theme_minimal() +
  labs(title = "Map with cities") +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    axis.title = element_blank()
  )
m4

# All plots together
all_in_one <- (m1 | m2) /
              (m3 | m4)

# Write the image
img_name <- "23-process.png"
ggsave(img_name, plot = all_in_one, height = 20, width = 16, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 24

# Load packages
library(dplyr)
library(stringr)
library(sf)
library(rnaturalearth)
library(ggplot2)
library(ggrepel)
library(magick)

# Download map of Chile
chile <- ne_countries(country = "Chile", scale = "medium", returnclass = "sf") |> 
  #st_make_valid() |> 
  st_crop(y = st_bbox(c(xmin = -77, ymin = -57, xmax = -65, ymax = -17)))

# Places listed manually
places <- c(
  #"Ahinco",
  "Antuco;-37.178;-71.620;Antu ko, agua de sol.",
  #"Arauco",
  "Batuco;-33.2238;-70.8062;Vatruko, agua de totora.",
  "Butaco;-37.74846;-72.64778;Fut<c3><a1> ko, estero grande.",
  "Calbuco;-41.7616;-73.1277;Kallf<c3><bc> ko, agua azul.",
  #"Catapilco",
  "Chacabuco;-33.02799;-70.69054;Estero del chacay.",
  "Cherquenco;-38.514;-72.026;Agua del cherc<c3><a1>n.",
  "Chanco;-35.682;-72.523;Cha<c3><b1> ko, aguas de cha<c3><b1>.",
  "Choshuenco;-39.8238;-72.0834;Aguas amarillas.",
  "Coihueco;-36.642;-71.840;Koiwe ko, Agua de coihue.",
  "Coinco;-34.207;-70.941;Coy<c3><ba>n ko, agua del arenal.",
  #"Colaco",
  "Collico;-39.7911;-73.2033;Aguas rojas.",
  "Coltauco;-34.007;-71.125;Koltraw ko, agua de renacuajos.",
  "Culenc<c3><b3>;-36.7086;-72.7377;Agua del cul<c3><a9>n.",
  #"Cumeico;;;Kumei ko, aguas buenas.",
  "Cunco;-38.9039;-72.0374;Agua clara.",
  "Curic<c3><b3>;-34.9580;-71.2395;Aguas negras.",
  #"Curi<c3><b1>anco",
  "Gualleco;-35.147;-71.955;Walleko, agua de hualles.",
  #"Huasco",
  "Lanco;-39.44746;-72.77636;Aguas tranquilas.",
  "Llico;-34.76566;-72.07456;Lli ko, desaguadero.",
  "Lumaco;-38.16415;-72.90630;Luma ko, agua de luma.",
  #"Malleco;;;Agua gredosa.",
  "Malloco;-33.60418;-70.85152;Mallo ko, agua de greda dorada.",
  "Melipeuco;-38.780;-71.675;Encuentro de cuatro aguas.",
  "Nantoco;-27.5241;-70.2679;Agua de pozo.",
  "Paillaco;-40.0576;-72.8709;Pailla ko, aguas tranquilas.",
  "Peleco;-36.96265;-72.67439;Peleko, agua con barro.",
  #"Pemuco",
  "Penco;-36.73434;-72.99368;Pen ko, divisando el agua.",
  "Perquenco;-38.376;-72.375;Perken ko, agua sanadora.",
  "Quidico;-38.24640;-73.48566;Agua solitaria.",
  "Quilaco;-37.68165;-72.00542;Quila ko, tres r<c3><ad>os.",
  "Ranco;-40.32063;-72.48148;Agua con oleaje.",
  "Renaico;-37.6645;-72.5832;Agua de cueva.",
  "Temuco;-38.7380;-72.5949;Temu ko, Agua de temu."
)

# Transform to data frame
places <- as.data.frame(do.call(rbind, str_split(places, pattern = ";")))
names(places) <- c("name", "latitude", "longitude", "description")

# Transform to sf
places <- places |>
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

# Add coordinates to the data frame for ggrepel
places <- places |> 
  mutate(X = st_coordinates(places)[ , 1], Y = st_coordinates(places)[ , 2])

# Do the plot
m <- ggplot() +
  geom_sf(data = chile, fill = "gray", color = NA) +
  geom_sf(data = places) +
  geom_text_repel(data = places, aes(x = X, y = Y, label = name), size = 4,
      max.overlaps = 100) +
  ylim(-42, -27) +
  theme(
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90)
  )
m

# Save the map as PNG file
img_name <- "24-names.png"
ggsave(img_name, plot = m, width = 8, height = 10, dpi = 300,
  bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 27

# Load packages
library(sf)
library(dplyr)
library(rmapshaper)
library(rnaturalearth)
library(units)
library(ggplot2)
library(RColorBrewer)
library(magick)

# Download South American map
sam <- ne_countries(continent = "South America", scale = "large",
    returnclass = "sf") |> 
  st_crop(y = c(xmin = -83.6, ymin = -57, xmax = -32.8, ymax = 13.25)) |>
  st_cast("POLYGON") |> 
  st_make_valid() |> 
  mutate(pid = row_number())

# Select only polygons with neighbours
sam_n <- sam |> 
  mutate(neighbours = lengths(st_touches(sam))) |> 
  filter(neighbours != 0)

# Get pairs
pairs <- st_touches(sam_n)

pairs_l <- list()
for(i in seq_along(pairs))
  pairs_l[[i]] <- cbind(i, pairs[[i]])
pairs_l <- do.call(rbind, pairs_l)
for (i in seq_along(pairs_l[ , 1]))
  pairs_l[i, ] <- pairs_l[i, order(pairs_l[i, ])]
pairs_l <- unique(as.data.frame(pairs_l))

# convert to pid
for (i in seq_along(pairs_l[1, ]))
  pairs_l[ , i] <- sam_n$pid[match(pairs_l[ , i], 1:nrow(sam_n))]

# Extract Borders Pairwise
sam_borders <- list()

for(i in seq_along(pairs_l[ , 1]))
  sam_borders[[i]] <- sam_n  |>
    filter(sam_n$pid %in% unlist(pairs_l[i, ])) |> 
    ms_simplify(keep = .2, keep_shapes = TRUE) |> 
    ms_innerlines() |> 
    as_tibble() |> 
    st_as_sf() |> 
    mutate(p1 = pairs_l[i, 1], p2 = pairs_l[i, 2])

sam_borders <- do.call(rbind, sam_borders) |> 
  mutate(lid = row_number()) |> 
  mutate(
    ctr1 = sam_n$adm0_a3[match(p1, sam_n$pid)],
    ctr2 = sam_n$adm0_a3[match(p2, sam_n$pid)]) |> 
  mutate(group = paste(ctr1, ctr2, sep = ".")) |> 
  group_by(group) |> 
  summarise(geometry = st_cast(st_combine(geometry), "MULTILINESTRING"))

# Add length and new identifier
sam_borders <- sam_borders |> 
  mutate(length = as.numeric(st_length(sam_borders))/1000, lid = row_number()) |> 
  filter(!grepl("SPI", group))

# Overview of borders
ggplot(data = sam_borders) +
  geom_sf(data = sam) +
  geom_sf(data = sam_borders, color = "red") +
  facet_wrap(~lid)

# Map with coloured lines according to their extension

continent <- "#770000ff"

m <- ggplot(data = sam_borders, aes(color = length)) +
  geom_sf(data = sam, fill = continent, color = continent) +
  geom_sf(data = sam_borders, size = 0.8) +
  scale_color_gradientn(colors = rev(brewer.pal(9, "YlOrRd")), name = "Length\n(km)") +
  theme_minimal()
m

# Save the map as PNG file
img_name <- "27-boundaries.png"
ggsave(img_name, plot = m, height = 9, width = 6, dpi = 300, bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)

Day 30

# Load packages
library(dplyr)
library(stringr)
library(sf)
library(elevatr)
library(terra)
library(tidyterra)
library(rnaturalearth)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(magick)

# Places listed manually
places <- c(
  # CHILE --------------------------------------------------------------------------------
  #"Ahinco",
  "Antuco;-37.178;-71.620;Antu ko, agua de sol.",
  #"Arauco",
  "Batuco;-33.2238;-70.8062;Vatruko, agua de totora.",
  "Batuco;-31.9347;-70.6119;",
  "Butaco;-37.74846;-72.64778;Fut<c3><a1> ko, estero grande.",
  "Calbuco;-41.7616;-73.1277;Kallf<c3><bc> ko, agua azul.",
  "Calbuco;-40.412769;-72.721596;",
  #"Catapilco",
  "Chacabuco;-33.02799;-70.69054;Chacay bucu, Estero del chacay.",
  "Cherquenco;-38.514;-72.026;Agua del cherc<c3><a1>n.",
  "Chanco;-35.682;-72.523;Cha<c3><b1> ko, aguas de cha<c3><b1>.",
  "Chanco;-42.771;-73.718;",
  "Chanco;-40.30964;-72.19957;",
  "Chanco;-38.275;-72.946;",
  "Chanco;-39.95785;-72.74837;",
  "Chanco;-40.42323;-73.32567;",
  "Choshuenco;-39.8238;-72.0834;Aguas amarillas.",
  "Coihueco;-36.642;-71.840;Koiwe ko, Agua de coihue.",
  "Coihueco;-38.07828;-73.30915;",
  "Coihueco;-38.42747;-71.89427;",
  "Coihueco;-40.9654;-73.2733;",
  "Coihueco;-35.78819;-72.43655;",
  "Coihueco;-38.50277;-72.57671;",
  "Coinco;-34.207;-70.941;Coy<c3><ba>n ko, agua del arenal.",
  #"Colaco",
  "Collico;-39.7911;-73.2033;Aguas rojas.",
  "Collico;-37.97377;-73.41940;",
  "Collico;-39.34333;-72.21987;",
  "Collico;-38.48104;-71.85265;",
  "Collico;-39.89258;-72.79232;",
  "Collico;-40.43421;-72.75288;",
  "Collico;-38.08120;-72.41110;",
  "Collico;-39.14893;-73.13247;",
  "Collico;-39.41862;-72.69306;",
  "Coltauco;-34.007;-71.125;Koltraw ko, agua de renacuajos.",
  "Culenc<c3><b3>;-34.79118;-71.69583;Agua del cul<c3><a9>n.",
  "Culenco;-36.72606;-72.73790;",
  "Culenco;-37.40616;-72.79562;",
  "Culenco;-36.12973;-72.57693;",
  "Culenco;-37.00999;-72.18151;",
  "Culenco;-36.30579;-72.28064;",
  #"Cumeico;;;Kumei ko, aguas buenas.",
  "Cunco;-38.9039;-72.0374;Agua clara.",
  "Cunco;-41.80216;-73.37734;",
  "Curic<c3><b3>;-34.9580;-71.2395;Aguas negras.",
  #"Curi<c3><b1>anco",
  "Gualleco;-35.147;-71.955;Walleko, agua de hualles.",
  #"Huasco",
  "Lanco;-39.44746;-72.77636;Aguas tranquilas.",
  "Llico;-34.76566;-72.07456;Lli ko, desaguadero.",
  "Llico;-37.19875;-73.56531;",
  "Lumaco;-38.16415;-72.90630;Luma ko, agua de luma.",
  "Lumaco;-39.714796;-72.438102;",
  "Lumaco;-39.06325;-72.70522;",
  "Lumaco;-39.98046;-73.28662;",
  "Lumaco;-40.01605;-72.79103;",
  #"Malleco;;;Agua gredosa.",
  "Malloco;-33.60418;-70.85152;Mallo ko, agua de greda dorada.",
  "Melipeuco;-38.780;-71.675;Encuentro de cuatro aguas.",
  #"Nantoco;-27.5241;-70.2679;Agua de pozo.",
  "Paillaco;-40.0576;-72.8709;Pailla ko, aguas tranquilas.",
  "Paillaco;-40.01605;-72.79103;",
  "Peleco;-36.96265;-72.67439;Peleko, agua con barro.",
  "Peleco;-37.87008;-73.37811;",
  #"Pemuco",
  "Penco;-36.73434;-72.99368;Pen ko, divisando el agua.",
  "Perquenco;-38.376;-72.375;Perken ko, agua sanadora.",
  "Quidico;-38.24640;-73.48566;Agua solitaria.",
  "Quidico;-37.35436;-73.60080;",
  "Quilaco;-37.68165;-72.00542;Quila ko, tres r<c3><ad>os.",
  "Quilaco;-37.35436;-73.60080;Agua de quila.",
  "Quilaco;-36.70627;-72.70984;",
  "Quilaco;-38.4124;-72.6272;",
  "Quilaco;-38.05237;-73.04318;",
  "Quilaco;-38.13481;-72.40022;",
  "Ranco;-40.32063;-72.48148;Agua con oleaje.",
  "Renaico;-37.6645;-72.5832;Agua de cueva.",
  "Temuco;-38.7380;-72.5949;Temu ko, Agua de temu.",
  # ARGENTINA ----------------------------------------------------------------------------
  "Angaco;-31.1893;-68.1139;Corrientes.",
  "Atreuc<c3><b3>;-37.1275;-63.8144;Aguas malas (salobres).",
  "Atreuco;-39.73010;-71.07828;",
  "Chapalc<c3><b3>;-38.64456;-63.09620;Agua barrosa.",
  "Claromec<c3><b3>;-38.85134;-60.06775;Tres arroyos.",
  "Comic<c3><b3>;-41.071336;-67.476976;Aguas escondidas.",
  "Cuchillo C<c3><b3>;-38.333333;-64.616667;Agua de chuchillos.",
  "Huinca Renanc<c3><b3>;-34.84107;-64.37490;Pozo de agua del cristiano.",
  "Huinganco;-37.15797;-70.62117;Agua de huing<c3><a1>n.",
  "Leubuc<c3><b3>;-36.87205;-63.03994;Lewf<c3><bc> ko, aguada del r<c3><ad>o.",
  #"Mach<c3><b3>nico;-40.332581;-71.417677;Lugar de p<c3><a1>ncora.",
  "Naic<c3><b3>;-36.92368;-64.40477;Manantial que baja.",
  "<c3><91>orquinc<c3><b3>;-41.8491;-70.9020;Agua de <c3><b1>orquin.",
  "Pehuen-C<c3><b3>;-39.002;-61.545;Agua de pehu<c3><a9>n.",
  "Realic<c3><b3>;-35.03752;-64.24633;Agua en forma de plato."
)

# Transform to data frame
places <- as.data.frame(do.call(rbind, str_split(places, pattern = ";")))
names(places) <- c("name", "latitude", "longitude", "description")

# Transform to sf
places <- places |>
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

# Add coordinates to the data frame for ggrepel
places <- places |> 
  mutate(X = st_coordinates(places)[ , 1], Y = st_coordinates(places)[ , 2])


# Define bounding box
bbox <- c(xmin = -77, xmax = -58, ymin = -43.5, ymax = -30.5)

# Load DEM
dem <- get_elev_raster(locations = st_sf(st_as_sfc(st_bbox(bbox, crs = 4326))),
    z = 7, clip = "locations")

# Load continent
sam <- ne_countries(continent = "South America", scale = "large",
    returnclass = "sf") |> 
  #filter(name %in% c("Chile", "Argentina")) |> 
  st_make_valid() |> 
  group_by(continent) |> 
  summarise(geometry = st_union(geometry))

# Mask dem with polygon
dem_masked <- mask(dem, sam)

# Calculate Hillshades
slope <- terrain(rast(dem_masked), "slope", unit = "radians")
aspect <- terrain(rast(dem_masked), "aspect", unit = "radians")
hill <- shade(slope, aspect, 30, 270)

# normalize names
#names(hill) <- "shades"

# Hillshading, but we need a palette
pal_greys <- hcl.colors(1000, "Grays")

# Main Map
m <- ggplot() +
  geom_sf(data = st_sf(st_as_sfc(st_bbox(bbox, crs = 4326))), fill = "#364891ff",
    color = NA) +
  geom_spatraster(data = hill, show.legend = FALSE) +
  scale_fill_gradientn(colors = pal_greys, na.value = NA) +
  geom_spatraster(data = rast(dem_masked), alpha = 0.7, show.legend = FALSE) +
  # scale_fill_gradientn(colors = grad_hypso, na.value = NA) +
  geom_sf(data = places, color = "white") +
  geom_text_repel(data = places, aes(x = X, y = Y, label = name),
    color = "white", size = 4, max.overlaps = 100) +
  xlim(bbox["xmin"], bbox["xmax"]) +
  ylim(bbox["ymin"], bbox["ymax"]) +
  theme(
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90)
  ) +
  theme_map()
m

# Minimap
bbox2 <- c(xmin = -82, xmax = -32, ymin = -60, ymax = -10)

mm <- ggplot() +
  # geom_sf(data = st_sf(st_as_sfc(st_bbox(bbox2, crs = 4326))), fill = "#364891ff",
  #   color = NA) +
  geom_sf(data = sam, fill = "gray", color = NA) +
  geom_sf(data = st_sf(st_as_sfc(st_bbox(bbox, crs = 4326))), fill = NA, color = "yellow",
      size = 1) +
  xlim(bbox2["xmin"], bbox2["xmax"]) +
  ylim(bbox2["ymin"], bbox2["ymax"]) +
  theme_map() +
  theme(
    panel.background = element_rect(fill = "#364891ff", color = NA)
  )
mm

# Final image
map <- ggdraw(m)  +
          draw_plot(mm, x = 0.72, y = 0.64, width = 0.2, height = 0.2)
map

# Save the map as PNG file
img_name <- "30-makeover.png"
ggsave(img_name, plot = map, width = 10, height = 10, dpi = 300,
  bg = "transparent")

# Trim transparent edges
img <- image_trim(image = image_read(img_name))
image_write(image = img, path = img_name)