Opening a spatial subset with {sf}

Day 3 and 4 of 30DayMapChallenge
R
30DayMapChallenge
datavisualization
spatial
Author

Michaël

Published

2022-11-04

Modified

2024-04-12

A photo of A view just outside the village of Arsos, Limassol

Wine villages – CC BY-NC by Lefteris Katsouromallis

Day 3 and 4 of 30DayMapChallenge: « polygons » and « green » (previously).

Intersecting an area of interest with a layer at opening time

The CORINE Landcover dataset is distributed as a geopackage weighting more than 8 Go. To limit the memory used when we only work on a subset, we can clip it at opening time. Here we will map the Cyprus Island:

library(dplyr)
library(ggplot2)
library(stringr)
library(sf)
library(rnaturalearth)
library(glue)

# Using the contour of Cyprus (enlarged) from naturalearth to clip
bb <- ne_countries(scale = "medium", country = "cyprus", returnclass = "sf") |> 
  st_transform("EPSG:3035") |> 
  st_buffer(90000) |> 
  pull(geometry) |> 
  st_as_text()

# Corine Land Cover 
# download from https://land.copernicus.eu/pan-european/corine-land-cover/clc2018
# (registration required)
# passing the bounding area
cyprus_clc <- read_sf("data/U2018_CLC2018_V2020_20u1.gpkg", query = glue("
  SELECT * 
  FROM U2018_CLC2018_V2020_20u1
  WHERE st_intersects(Shape, st_polygonfromtext('{bb}'))"))

legend_colors <- list("TRUE"  = "mediumaquamarine",
                      "FALSE" = "grey90")

legend_labs <- list("TRUE"  = "forest",
                    "FALSE" = "other")

# exclude sea (code 523)
# and classify on forest (codes 3xx)
cyprus_clc |> 
  filter(Code_18 != "523") |> 
  ggplot() +
  geom_sf(aes(fill = str_detect(Code_18, "^3"),
              color = str_detect(Code_18, "^3"))) +
  scale_fill_manual(name= "type",
                    values = legend_colors,
                    labels = legend_labs) +
  scale_color_manual(name = "type",
                     values = legend_colors,
                     labels = legend_labs) +
  labs(title = "Cyprus",
       subtitle = "Landcover",
       caption = glue("data: Copernicus CLC 2018
                      projection LAEA
                      r.iresmi.net {Sys.Date()}")) +
  theme_minimal() +
  theme(legend.position = "bottom",
               plot.caption = element_text(size = 7))

A map of Cyprus forests

Cyprus landcover