Postal codes and tessellations

Wait 5 sec.

[This article was first published on One World | R packages, maps and open-source projects - R Project, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.From postal points to areas using Voronoi polygons 13 min.I recently published my project Postal codes ofSpain, with estimated postalareas, an interactive map and GeoPackage downloads. But there is something wehave not looked at yet: how do those address points become polygons? The answeris… Voronoi polygons.What is a Voronoi polygon?Imagine a few shops on a flat map. Which one is closest to your house? Ask thatquestion for every location and give each shop its own color. The result is aVoronoi diagram, also known as a Voronoi tessellation, Voronoidecomposition, Voronoi partition or Dirichlet tessellation. Thistechnique divides a plane into regions based on the distance to a given set ofpoints.In our example, each shop is a site, or generating point. Its Voronoicell, also called a Thiessen polygon, contains the locations at least asclose to that site as to any other. Along a shared edge, two sites areequally close. Where several cells meet, there may be a tie between more sites.With ordinary straight-line distance, the edges follow perpendicular bisectorsbetween sites. Some cells extend indefinitely, so maps usually show them clippedto a study area.The name honors Georgy Voronoi, who studied these partitions in the earlytwentieth century. The idea has earlier roots, including Dirichlet’snineteenth-century work (hence the name Dirichlet tessellation). Pokojski andPokojska (2018) review the history andgeographical applications.These partitions are useful for exploring proximity to shops or public services,assigning areas to weather stations and studying spatial patterns. For example,Thiessen polygons can supply area weights when estimating rainfall from raingauges. The underlying question is always about the chosen distance: the nearesthospital in a straight line may not be the quickest one to reach.A tiny example in RLet’s start with five points in an imaginary square, including one at itscenter, using just sf and ggplot2. The coordinates have arbitrary unitsand no geographical coordinate reference system.library(sf)library(ggplot2)# Define a common theme for the maps.background select(POSTCODE)code_counts 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) ))) Indicator Value Municipalities 212 Downloaded points 557 Points with a postal code 557 Distinct postal codes 557 Codes represented by one point 557 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 coordinatesare correct. With one point per code, we will be drawing areas aroundrepresentatives rather than estimating them from individual addresses.Two points at exactly the same coordinates cannot generate separate Voronoicells. We sort by postal code and keep the first point at each location, as inthe original workflow. This is a reproducible tie-breaking rule, not evidencethat the retained code is correct. The counts below show how many records weremove and how many locations have conflicting codes.postal_points arrange(POSTCODE)conflicts group_by(geometry) |> summarise(codes = n_distinct(POSTCODE), .groups = "drop") |> filter(codes > 1)duplicates