Mapping building use with a hexagonal grid

I needed a compact way to show the composition of building uses across Spain without pixel‑level clutter. We will aggregate 100 m building‑use rasters to a 20 km hexagonal grid and visualize the mix of agricultural, industrial, and commercial uses with overlapping proportional symbols blended by multiplication. The Canary Islands are shown as an inset.

Tip Hexagons reduce directional bias compared to squares and improve the perception of regional patterns when summarising rasters.

Packages

pkgs <- c("terra","sf","tidyverse","fs","mapSpain","ggforce",
          "patchwork","scales","ggblend","rmapshaper","classInt")
for(p in pkgs) if(!require(p, character.only = TRUE)) install.packages(p)

library(terra)
library(sf)
library(tidyverse)
library(fs)
library(mapSpain)
library(ggforce)
library(patchwork)
library(scales)
library(ggblend)
library(rmapshaper)
library(classInt)

Data preparation

Study area

We obtain autonomous community boundaries, project to EPSG:25830, and exclude the Canary Islands for the main map.

prov <- esp_get_ccaa(moveCAN = FALSE) |>
  st_transform(25830) |>
  filter(ccaa.shortname.ga != "Canarias")

Input rasters and reclassification

We read the building use rasters from 2020, treat zeros as missing, and aggregate layers to build the comercial category (commerce + services + offices). The data come from our 2023 paper; see 10.5194/essd-2023-53 for details. Download data here, and for Canary Island here.

# import raster layers  with regexp if all years are downloaded
landuse <- dir_ls("spain", regexp = "2020.tif$") |> rast()

landuse[landuse == 0] <- NA

# subset and aggregate uses
comercial <- sum(subset(landuse, c(2, 4, 5)), na.rm = TRUE)

# joint raster layers new aggregate use with agriculture and industrial
landuse <- c(subset(landuse, c(1, 2)), comercial)
names(landuse) <- c("agriculture", "industrial", "comercial")

landuse 

class       : SpatRaster 
size        : 9998, 11927, 3  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : -40877.9, 1151822, 3880869, 4880669  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
sources     : spat_397864b84d08_14712_zxzxx9jyyx6UnnI.tif  (2 layers) 
              memory  
names       : agriculture, industrial, comercial 
min values  :           1,          1,         1 
max values  :          58,         44,        64 

Hexagonal grid

We create a 20 km hex grid that covers the study area and add an id. To build the spatial units used for aggregation, we rely on st_make_grid(), a very flexible function from the sf package that allows us to generate regular grids directly from any input geometry. By default, the function creates square cells, but by setting square = FALSE we instead request hexagonal polygons, which are often preferred in spatial aggregation because they reduce directional bias and provide a more uniform neighborhood structure. The argument cellsize = units::set_units(20, km) specifies the diameter (edge‑to‑edge distance) of each hexagon, ensuring all cells have comparable area across the entire extent. We then assign each hexagon a unique identifier, which is later used to join aggregated raster values back to the grid.

grid <- st_make_grid(prov, cellsize = units::set_units(20, km), 
                     square = FALSE) |>
  st_as_sf() |>
  mutate(poly_id = row_number())

Extract raster values by hexagon

We extract raster values for each hexagon and compute category totals.

extr_pixels <- terra::extract(landuse, vect(grid))

abudf <- extr_pixels |>
  pivot_longer(2:4, names_to = "use", values_to = "nbu") |>
  group_by(ID, use) |>
  summarise(nbu = sum(nbu, na.rm = TRUE), .groups = "drop")

Compute within‑hexagon proportions

In this step we rejoin the extracted raster data to the grid. To make the map speak about composition rather than sheer magnitude, we convert absolute counts to relative shares.

grid <- left_join(grid, abudf, by = join_by(poly_id == ID)) |>
           filter(nbu != 0) |>
            group_by(poly_id) |>
             mutate(nbu_rel = nbu * 100 / sum(nbu, na.rm = TRUE)) |>
              ungroup()

We define discrete classes for proportional symbol sizes.

cl <- c(0, 15, 30, 60, 80, 100)

Administrative context

We get internal borders from the grid by assigning each hexagon to the community with the largest overlap.

# regional boundaries
ccaa <- esp_get_ccaa() |>
  st_transform(25830)

# intersection between grid and boundaries
lim <- st_intersection(grid, ccaa)
inter <- lim |> mutate(area = as.numeric(st_area(lim)))

# Which hexagon corresponds to which region?
asig <- inter |>
  arrange(poly_id, desc(area)) |>
  group_by(poly_id) |>
  slice_head(n = 1) |>
  ungroup() |>
  select(poly_id, ine.ccaa.name)

grid <- left_join(grid, st_drop_geometry(asig), by = "poly_id")

# dissolve geometry by region
grid_lim <- group_by(grid, ine.ccaa.name) |>
                 summarise()

plot(grid_lim)

Background hexes and locator

We create background hexes to indicate the full coverage and a small segment to locate the Canary Islands.

grid_aux <- st_make_grid(prov, cellsize = units::set_units(20, km), square = FALSE)
plot(grid_aux)

missing_hex <- st_intersects(grid_aux, filter(prov, nuts2.name != "Canarias"))
miss <- grid_aux[lengths(missing_hex) != 0] |>
  st_as_sf() |>
  dplyr::mutate(id = dplyr::row_number()) |>
  dplyr::filter(!id %in% c(1465, 899, 1012, 743)) # exclude hexagons of super small islands 

plot(miss)

# canary separation line for the inset map
can_box <- st_linestring(matrix(c(-8.238606, 36.756331,
                                  -6.2,      35.532012),
                                nrow = 2, byrow = TRUE)) |>
  st_sfc(crs = 4326) |>
  st_transform(25830)

Constructing the map

Additive legend

We draw three overlapping circles as a reminder of the color mapping for each use. To build the additive legend, we use geom_circle() from ggforce, which lets us draw perfect circles by specifying a center (x0, y0) and a radius (r). This makes it easy to place three uniformly sized circles representing the three land‑use categories. Since the circles intentionally overlap, we apply blend("multiply") from the ´ggblend´ package, a blending mode that darkens the intersecting areas. This visually mimics how the mixed-use symbols behave on the main map: the more overlap, the stronger the combined color—reinforcing the idea of composition rather than isolated categories.

Since blend() operates on the previous graphical layer, and here we need to blend multiple layers together, all the geom_circle() calls must be wrapped inside a list. Passing this list through the pipe to blend("multiply") instructs to treat the group of layers as a single blended unit, ensuring that the overlapping colours combine correctly in the legend—just as they do in the main map.

# color and labels
cols <- c(agricultural = "#f5d300", comercial = "#08f7fe", industrial = "#ff1493")

# artificial data.frame for circles 
legend_df <- data.frame(
  x = c(-0.4, 0.4, 0.0),
  y = c(0.0, 0.0, 0.35),
  r = c(0.55, 0.55, 0.55),
  key = c("agricultural", "comercial", "industrial"),
  stringsAsFactors = FALSE
)

# angle position for labels radians!!
angles <- c(agricultural = 5*pi/4, comercial = 2*pi, industrial = 3*pi/4)
label_vec <- c(agricultural = "Agricultural", comercial = "Commercial", industrial = "Industrial")

# calculate x and y for each label
pos_labels <- within(legend_df, {
  key <- as.character(key)
  theta <- angles[key]
  lab1_x <- x + r * cos(theta)
  lab1_y <- y + r * sin(theta)
  lab2_x <- x - r * cos(theta)
  lab2_y <- y - r * sin(theta)
  lab1 <- label_vec[key]
  lab2 <- label_vec[key]
})

# legend
legend_plot <- ggplot() +
  list(
    geom_circle(data = subset(legend_df, key == "agricultural"),
                         aes(x0 = x, y0 = y, r = r),
                         fill = cols["agricultural"], color = NA, alpha = 0.85),
    geom_circle(data = subset(legend_df, key == "comercial"),
                         aes(x0 = x, y0 = y, r = r),
                         fill = cols["comercial"], color = NA, alpha = 0.85),
    geom_circle(data = subset(legend_df, key == "industrial"),
                         aes(x0 = x, y0 = y, r = r),
                         fill = cols["industrial"], color = NA, alpha = 0.85)
  ) |>
  blend("multiply") +
  geom_label(data = pos_labels, aes(lab1_x, lab1_y, label = lab1),
             fill = alpha("white", .5), size = 5, linewidth = 0) +
  coord_equal(xlim = c(-1.3, 1.3), ylim = c(-0.6, 1.2), expand = FALSE, clip = "off") +
  theme_void(base_family = "sans")

legend_plot

Main map

We plot three proportional‑symbol layers (one per use) at the centroids of each hexagon and blend them using multiplication in the same way as for the legend. To ensure that the proportional symbols remain visually coherent within the map, the size of the circles must be constrained so that even the largest symbol does not exceed the width of a hexagon. Since the symbols are plotted at the centroid of each cell, oversized circles would spill over into neighbouring hexagons, creating visual clutter and making category comparisons harder to read.

g <- ggplot() +
  geom_sf(data = miss, fill = "grey90", color = "white") +
  geom_sf(data = can_box, colour = "grey85") +
  geom_sf(data = ms_innerlines(grid_lim), colour = "white", linewidth = 1.1) +
  list(
    geom_sf(data = st_centroid(grid) |>
              dplyr::filter(use == "agriculture"),
            aes(color = use, size = cut(nbu_rel, cl))),
    geom_sf(data = st_centroid(grid) |>
              dplyr::filter(use == "comercial"),
            aes(color = use, size = cut(nbu_rel, cl))),
    geom_sf(data = st_centroid(grid) |>
              dplyr::filter(use == "industrial"),
            aes(color = use, size = cut(nbu_rel, cl)))
  ) |>
  ggblend::blend("multiply") +
  scale_color_manual(values = c("agriculture" = "#f5d300",
                                "comercial"   = "#08f7fe",
                                "industrial"  = "#ff1493"),
                     breaks = c("agriculture","comercial","industrial"),
                     labels = c("Agricultural","Commercial","Industrial")) +
  scale_size_manual(values = c(1, 2, 4, 6, 8) * 0.62,
                    labels = c("< 15","15–30","30–60","60–80","80–100 %")) +
  labs(size = "Share of buildings",
       title = stringr::str_wrap("Geography of building use in Spain | 2020", 60),
       subtitle = stringr::str_wrap(
         "Residential use is excluded. 'Commercial' merges public services and office uses to simplify spatial interpretation and comparability.",
         90),
       caption = "HISDAC-ES. Uhl et al. 2023. 10.5194/essd-2023-53") +
  guides(color = "none") +
  theme_void(base_family = "sans") +
  theme(legend.position = c(.83, .90),
        legend.title.position = "top",
        legend.text.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 18, face = "bold"),
        plot.margin = margin(20, 20, 20, 300))

Canary Islands inset

We process the Canary Islands separately (in REGCAN95 projection) and keep sizes and colours consistent.

landuse_can <- dir_ls("canary", regexp = "2020.tif$") |> rast()

landuse_can[landuse_can == 0] <- NA
comercial_can <- sum(subset(landuse_can, c(2, 4, 5)), na.rm = TRUE)
landuse_can <- c(subset(landuse_can, c(1, 3)), comercial_can)
names(landuse_can) <- c("agriculture", "industrial", "comercial")

lim_can <- esp_get_ccaa(moveCAN = FALSE) |>
             st_transform(4083) |>
              filter(ccaa.shortname.es == "Canarias") 

grid_can <- st_make_grid(lim_can, cellsize = units::set_units(20, km), square = FALSE) |>
  st_as_sf() |>
  mutate(poly_id = row_number())

extr_pixels <- terra::extract(landuse_can, vect(grid_can))

abudf_can <- extr_pixels |>
  pivot_longer(2:4, names_to = "use", values_to = "abu") |>
  group_by(ID, use) |>
  summarise(abu = sum(abu, na.rm = TRUE), .groups = "drop")

grid_can <- dplyr::left_join(grid_can, abudf_can, by = join_by(poly_id == ID)) |>
             filter(abu != 0) |>
              group_by(poly_id) |>
               mutate(abu_rel = abu * 100 / sum(abu, na.rm = TRUE)) |>
                ungroup()

map_can <- ggplot() +
  geom_sf(data = grid_can, fill = "grey90", color = "white") +
  list(
    geom_sf(data = st_centroid(grid_can) |>
              dplyr::filter(use == "agriculture"),
            aes(color = use, size = cut(abu_rel, cl))),
    geom_sf(data = st_centroid(grid_can) |>
              dplyr::filter(use == "comercial"),
            aes(color = use, size = cut(abu_rel, cl))),
    geom_sf(data = st_centroid(grid_can) |>
              dplyr::filter(use == "industrial"),
            aes(color = use, size = cut(abu_rel, cl)))
  ) |>
  ggblend::blend("multiply") +
  scale_color_manual(values = c("agriculture" = "#f5d300",
                                "comercial"   = "#08f7fe",
                                "industrial"  = "#ff1493"),
                     breaks = c("agriculture","comercial","industrial"),
                     labels = c("Agricultural","Commercial","Industrial")) +
  scale_size_manual(values = c(1, 2, 4, 6, 8) * 0.60) +
  labs(size = NULL, title = "Canary Islands") +
  guides(color = "none", size = "none") +
  theme_void(base_family = "sans") +
  theme(legend.position = c(.85, .87),
        legend.title.position = "top",
        legend.text.position = "bottom",
        legend.direction = "horizontal",
        plot.margin = margin(10, 10, 10, 10),
        plot.title = element_text(vjust = 1, hjust = .05, size = 11))

Compose and export

We insert the additive legend and the Canary Islands inset and export the final figure.

gfinal <- g +
  inset_element(legend_plot, .70, .80, .95, 1.3) +
  inset_element(map_can, 0, -0.2, .30, .34, align_to = "full")

ggsave("esp_buiuse_2020.png", gfinal,
       height = 13.3, width = 16, units = "in", bg = "white", dpi = 300)

What this map communicates

• Dominance: areas where one use exceeds ~60% are visible by symbol size and dominant colour.
• Interpretation caution: since the map shows the relative number of buildings by use, the categories reflect building functions, which do not necessarily correspond to the underlying economic activities of each area.
• Mixing: blending reveals coexistence of uses.
• Regional patterns: subtle borders help orientation without clutter.

Notes and good practices

• Keep rasters and vectors in a projected CRS before terra::extract().
• Compute centroids on projected geometries.
• Fixed size breaks (c(0, 15, 30, 60, 80, 100)) match legend labels
• Use exact=TRUE in extract() only if you need area weighting.