Postal codes and tessellations | R-bloggers

Getting the data

First, let’s get the municipal boundaries and postal points. We use the 2024
release of both datasets.

library(giscoR)
library(sf)
library(dplyr)

municipalities <- gisco_get_lau(
  year = 2024, country = "SI", epsg = 3035
) |>
  select(GISCO_ID, LAU_NAME) |>
  st_make_valid()

raw_points <- gisco_get_postalcodes(
  year = 2024, country = "SI", epsg = 3035
)

We work in ETRS89 / LAEA Europe (EPSG:3035), a projected coordinate system
with coordinates in meters. This makes the 200 m and 50 m margins used later
meaningful. Voronoi uses Euclidean distances in this plane, rather than road
distances. The projection does not preserve every distance exactly either.

Before we start drawing, what have we actually downloaded?

postal_points <- raw_points |>
  mutate(POSTCODE = trimws(as.character(POSTCODE))) |>
  filter(!is.na(POSTCODE), POSTCODE != "") |>
  select(POSTCODE)

code_counts <- postal_points |>
  st_drop_geometry() |>
  count(POSTCODE, name = "points_per_code")

knitr::kable(data.frame(
  Indicator = c(
    "Municipalities", "Downloaded points", "Points with a postal code",
    "Distinct postal codes", "Codes represented by one point"
  ),
  Value = c(
    nrow(municipalities), nrow(raw_points), nrow(postal_points),
    nrow(code_counts), sum(code_counts$points_per_code == 1)
  )
))

These counts give us a first look at the file, but they are not a quality check.
They do not establish whether every code is present or whether the coordinates
are correct. With one point per code, we will be drawing areas around
representatives rather than estimating them from individual addresses.

Two points at exactly the same coordinates cannot generate separate Voronoi
cells. We sort by postal code and keep the first point at each location, as in
the original workflow. This is a reproducible tie-breaking rule, not evidence
that the retained code is correct. The counts below show how many records we
remove and how many locations have conflicting codes.

postal_points <- postal_points |> arrange(POSTCODE)

conflicts <- postal_points |>
  group_by(geometry) |>
  summarise(codes = n_distinct(POSTCODE), .groups = "drop") |>
  filter(codes > 1)

duplicates <- duplicated(postal_points$geometry)

knitr::kable(data.frame(
  Indicator = c("Duplicate locations removed", "Locations with conflicting codes"),
  Value = c(sum(duplicates), nrow(conflicts))
))

postal_points <- postal_points[!duplicates, ]

Starting with Ljubljana

Let’s pick Ljubljana and keep the points inside it. We use st_within() so that
a point exactly on a shared boundary is not assigned to both neighboring
municipalities.

city <- municipalities |> filter(LAU_NAME == "Ljubljana")

city_points <- st_filter(postal_points, city, .predicate = st_within)

base_map <- ggplot() +
  geom_sf(data = city, fill = "#e8e4d8", colour = ink, linewidth = 0.4)

gisco_attr <- gisco_attributions()

base_map +
  geom_sf(data = city_points, colour = "#bd5137", size = 2) +
  labs(
    title = "Ljubljana: the starting point",
    subtitle = "A municipal boundary and points with postal codes",
    caption = gisco_attr
  )

Building the cells

Each cell contains locations at least as close to its generating point as to
any other point. Cell boundaries occur where distances to two points are equal.
For this first demonstration, we keep all interior points, without the margins used in
the national workflow.

cells <- city_points |>
  st_union() |>
  st_voronoi(envelope = st_as_sfc(st_bbox(city))) |>
  st_collection_extract("POLYGON") |>
  st_sf(geometry = _)

base_map +
  geom_sf(data = cells, fill = NA, colour = "#bd5137", linewidth = 0.5) +
  geom_sf(data = city_points, colour = ink, size = 2) +
  labs(
    title = "An area of proximity for every point",
    subtitle = "Cells still extend beyond the municipal boundary",
    caption = "Own elaboration using data from Eurostat / GISCO, 2024"
  )

st_union() gathers the points so we build one diagram. One detail to keep in
mind: envelope supplies a working extent, but does not clip the cells to
the municipality. We need an intersection for that.

Clipping and recovering postal codes

It would be convenient if cells came out in the same order as the input points,
but we should not rely on that. Instead, we recover each cell’s code with a
spatial relationship before clipping. Each cell must contain exactly one
generating point.

point_matches <- st_intersects(cells, city_points)

cells$POSTCODE <- city_points$POSTCODE[unlist(point_matches)]

clipped_geometry <- st_intersection(st_geometry(cells), st_geometry(city))

city_cells <- st_sf(
  POSTCODE = cells$POSTCODE[attr(clipped_geometry, "idx")[, 1]],
  geometry = clipped_geometry
)

city_areas <- city_cells |>
  group_by(POSTCODE) |>
  summarise(.groups = "drop") |>
  st_make_valid()

city_colours <- setNames(
  hcl.colors(nrow(city_areas), palette = "Dark 3"),
  sort(city_areas$POSTCODE)
)
city_labels <- st_sf(
  POSTCODE = city_areas$POSTCODE,
  geometry = st_point_on_surface(st_geometry(city_areas))
)

ggplot(city_areas) +
  geom_sf(aes(fill = POSTCODE), colour = background, linewidth = 0.6) +
  geom_sf(
    data = city_points, shape = 21, fill = ink,
    colour = background, size = 2.4
  ) +
  geom_sf_text(
    data = city_labels,
    aes(label = POSTCODE), colour = ink, fontface = "bold", size = 3.5
  ) +
  scale_fill_manual(values = city_colours, guide = "none") +
  coord_sf(datum = NA) +
  labs(
    title = "From points to estimated areas",
    subtitle = "Cells clipped and merged by postal code",
    caption = "Teaching example · Not official postal boundaries"
  )

We intersect the geometry columns, then use the intersection’s idx matrix to
recover each cell’s postal code. This keeps the attribute assignment explicit.
Label positions are also calculated from geometries alone.

Grouping by postal code with summarise() merges the geometries. With many
addresses sharing a code, their cells would dissolve into one geometry, possibly
with several separate parts.

From one municipality to the whole country

Ljubljana is done. Now let’s repeat the operation for each component of each
municipality. Municipalities can have separate pieces: normalizing to
MULTIPOLYGON and then casting to POLYGON lets us handle them independently,
without automatically extending a code to a piece with no points.

The municipality ID and name apply to every component, so we mark those
attributes as spatially constant before splitting the geometries.

components <- municipalities |>
  st_set_agr("constant") |>
  st_cast("MULTIPOLYGON") |>
  st_cast("POLYGON") |>
  mutate(component_id = row_number())

Why use an inward margin?

In the Spanish project, we first select addresses within an inward margin of
200 m. If none are found, we try 50 m. This reduces the influence of
points close to a boundary, where small mismatches between layers may affect
their assignment.

There is no universal rule behind these numbers. Goerlich describes a 20 m
margin. The 200 m and 50 m choices come from my implementation. We keep them
here to explain that workflow, although sparse postal points make it easier to
exclude useful information. Comparing other margins would be a useful next step.

select_points <- function(component, postal_points) 
  for (margin in c(200, 50)) 
    interior <- st_buffer(st_geometry(component), -margin)

    if (all(st_is_empty(interior))) next

    selected_points <- st_filter(postal_points, interior,
      .predicate = st_within
    )

    if (nrow(selected_points) > 0) 
      return(selected_points)
    
  

  postal_points[0, ]

The margin only selects points: the resulting areas cover the original
component, not the smaller buffered geometry. If both attempts fail, we leave
the component unassigned. We do not borrow the nearest point from outside the
municipality.

Three cases, one function

There are three cases to handle. With no points, we keep the component and
assign NA. If all selected points share a code, we assign it to the whole
component. We only need Voronoi when several codes are present: no need to split
an area that already has a single assignment.

build_areas <- function(component, postal_points) >
    st_collection_extract("POLYGON") 

result <- lapply(seq_len(nrow(components)), function(i) 
  build_areas(components[i, ], postal_points)
) |>
  bind_rows()

Handling unassigned areas

Can we fill any of the gaps? A component without points may belong to a
municipality whose other components have information. We recover its code only
if the municipality has one distinct recovered code. This is an imputation
from the available data, not proof that the municipality really has just one
code.

unique_codes <- result |>
  st_drop_geometry() |>
  filter(!is.na(POSTCODE)) |>
  distinct(GISCO_ID, POSTCODE) |>
  add_count(GISCO_ID) |>
  filter(n == 1)

unassigned_before <- sum(is.na(result$POSTCODE))

index <- match(result$GISCO_ID, unique_codes$GISCO_ID)

recoverable <- is.na(result$POSTCODE) & !is.na(index)

result$POSTCODE[recoverable] <- unique_codes$POSTCODE[index[recoverable]]

Before drawing the final map, let’s tidy up the geometries and keep only
surfaces. Repairs can produce collections containing residual lines or points,
which do not belong in an area layer. Dissolving by postal code also removes
municipal boundaries between areas sharing that code.

keep_surfaces <- function(x) 
  geometry <- st_make_valid(st_geometry(x))

  geometry <- lapply(geometry, function(part) 
    if (!inherits(part, c("POLYGON", "MULTIPOLYGON", "GEOMETRYCOLLECTION"))) 
      return(st_multipolygon())
    

    surfaces <- suppressWarnings(st_collection_extract(st_sfc(part), "POLYGON"))

    if (length(surfaces) == 0L )

  st_geometry(x) <- st_sfc(geometry, crs = st_crs(x))

  x[!st_is_empty(x), ]


unassigned <- result |>
  filter(is.na(POSTCODE)) |>
  keep_surfaces()

postal_areas <- result |>
  filter(!is.na(POSTCODE)) |>
  group_by(POSTCODE) |>
  summarise(.groups = "drop") |>
  keep_surfaces() |>
  st_cast("MULTIPOLYGON")

The complete map

For the national overview, let’s color the areas by the first digit of their
code. This keeps the map readable, without claiming that we have validated a
classification of postal regions. Gray is reserved for unassigned components,
with their own legend entry.

postal_colours <- c(
  "1" = "#007F7F", "2" = "#E87500", "3" = "#3B6FB6", "4" = "#A64E9A",
  "5" = "#649A32", "6" = "#D9AF00", "8" = "#B94350", "9" = "#6255A4",
  "Unassigned" = "#727A80"
)

postal_groups <- sort(unique(substr(postal_areas$POSTCODE, 1, 1)))

ggplot() +
  geom_sf(data = municipalities, fill = "#d5d2ca", colour = NA) +
  geom_sf(
    data = postal_areas,
    aes(fill = substr(POSTCODE, 1, 1)), colour = background, linewidth = 0.15
  ) +
  geom_sf(
    data = unassigned, aes(fill = "Unassigned"),
    colour = background, linewidth = 0.15
  ) +
  scale_fill_manual(
    values = postal_colours, name = "Postal code group",
    breaks = c(postal_groups, "Unassigned"),
    labels = c(paste0(postal_groups, "xxx"), "Unassigned")
  ) +
  coord_sf(datum = NA) +
  labs(
    title = "Slovenia, point by point",
    subtitle = "Proximity areas clipped to municipal boundaries",
    caption = paste("Eurostat / GISCO, 2024 · Own elaboration",
      "Dark grey: unassigned areas · Estimated boundaries, not official",
      sep = "\n"
    )
  )

Checking the result

The map is ready, but let’s check the geometry before going further. Are the
polygons valid? Do assigned and unassigned areas cover the territory? Are there
any appreciable overlaps? We allow one square meter for numerical differences in
the geometry operations.

territory <- st_union(st_geometry(municipalities))

coverage <- st_union(c(st_geometry(postal_areas), st_geometry(unassigned)))

difference_m2 <- as.numeric(sum(st_area(st_sym_difference(territory, coverage))))

overlap_m2 <- as.numeric(
  sum(st_area(postal_areas)) + sum(st_area(unassigned)) - sum(st_area(coverage))
)

knitr::kable(data.frame(
  Indicator = c(
    "Codes with an estimated area", "Components assigned by municipal recovery",
    "Unassigned records", "Unassigned area (km²)",
    "Coverage difference (m²)", "Overlap (m²)"
  ),
  Value = round(c(
    nrow(postal_areas), sum(recoverable), nrow(unassigned),
    as.numeric(sum(st_area(unassigned))) / 1e6, difference_m2,
    overlap_m2
  ), 4)
))

These checks verify the geometry, not postal accuracy. We could cover an
entire country with valid polygons and still have a poor approximation of its
delivery areas.

Back to Spain

We now have the full sequence: points, proximity areas, clipping and merging by
postal code. The same steps work with much more detailed addresses, such as the
CartoCiudad data in my Spanish postal codes
project.

The input data matters as much as the algorithm. With one representative point
per code, Voronoi mainly maps distances between representatives. With address
points, it can better approximate where each code’s addresses are concentrated.
It still knows nothing about delivery routes, physical barriers or official
boundaries.

We have not compared these results with an independent postal source for
Slovenia. Before using them for analysis, we would need to review coverage,
coordinates, the relationship with municipalities and sensitivity to the
margins. For now, we have a working example that shows how the polygons are
built, one step at a time.

References

• Pokojski, W. & Pokojska, P. (2018). Voronoi diagrams – inventor, method,
  applications. Polish Cartographical Review, 50(3), 141–150.
• Goerlich, F. J. (2022). Elaboración de un mapa de códigos postales de España
  con recursos libres. Cómo evitar pagar por disponer de información de
  referencia. Working Papers Ivie, no. 2022-3. València: Ivie.

  Updated online version:

• Sánchez, J. I. (2020, November 13). “Así hice el mapa de los códigos
  postales de España con SIG y datos abiertos.” NoSoloSIG.
• Hernangómez, D. (2026). Códigos postales de España [Dataset].

Related