GIS, a transformational approach

This quick post rounds out the GIS transformations from area data to other spatial types. It turns out—to my relief, if not my surprise—that a package is available that implements the suggested transformation from areas to field data in the table on page 26 of Geographic Information Analysis, making this one of the simpler posts in the series so far.

Code

library(sf)
library(terra)
library(predicts)
library(dplyr)
library(stringr)
library(patchwork)
library(ggplot2)
library(ggspatial)

theme <- theme_void() +
  theme(legend.position = "inside",
        legend.position.inside = c(0.5, 0.075),
        legend.direction = "horizontal",
        legend.title.position = "top",
        legend.key.width = unit(0.8, "cm"),
        legend.key.height = unit(0.5, "cm"),
        legend.text = element_text(angle = 45, hjust = 1),
        legend.background = element_rect(fill = "#ffffff80", colour = NA), 
        panel.background = element_rect(fill = NA),
        panel.ontop = TRUE)

From areas…

As always, we need some data. It’s surprisingly tricky to find an area of New Zealand where the population density is high(ish) and somewhat variable over an extensive area such that the result of density smoothing are likely to be interesting. Best I can do is some old census data by 100 or so Census Area Units for central Auckland back in 2006.

Code

context <- st_read("ak-context.gpkg")
polygons <- st_read("ak-2006-ethnicity-and-tb.gpkg") |> 
  mutate(
    area_km2 = (st_area(geom) |> units::drop_units()) / 1e6,
    pop_density_km2 = population / area_km2)

Here’s what those data look like.

Code

g1 <- ggplot() +
  geom_sf(data = context, colour = NA) + 
  geom_sf(data = polygons, aes(fill = pop_density_km2), colour = NA) +
  scale_fill_distiller("Population density sq.km",
                       palette = "Purples", direction = 1) +
  coord_sf(expand = FALSE) +
  annotation_scale() +
  theme
g1

Figure 1: The SA2 population data, by count, and by density

There are a lot of people (many of them students) in apartments living at somewhat urban densities in the central city area, but the rest of this fairly extensive area (about 155 sq.km) is fairly suburban. Perhaps recent developments in Auckland are set to change that, although New Zealand has resisted densification for so long I wouldn’t bet the house on it.

… To fields

Anyway, we are here to create a pycnophylactic-ally smoothed surface from this. Pycnophylactic smoothing seems to me like a word that Waldo Tobler made up for a bet. In the original paper¹ he suggests that the term means ‘mass preserving’, although the prefix ‘pycno-’ apparently refers to thickness or density in the Greek from whence it came. Either way it’s easier to think of it as meaning ‘volume-preserving’, as the volume under the surface produced by this transformation is equal to that of prisms extruded to a height given by the areal density of the countable phenomenon of interest. This is a desirable property for any smoothing process to have as it means that the sum across all cells in the resulting raster surface will match the total across all regions in the source data.

Anyway, the mechanics of the process consist of repeated local averaging and rescaling to ensure that the quantity attributed to each region after each round of averaging matches the starting quantity. Repeated application of averaging and rescaling sees the population of denser regions migrate towards those regions’ centres and of less dense regions towards more densely populated neighbouring regions.

A function to perform this iterative process is provided by the predicts package. The pycnophy function seems like a bit of an afterthought in the package, which is mostly concerned with species distribution models. Running the function is straightforward. Convergence can take a while if you go for a very small cell size. In this case 100m seems about right given the extent of the data.

Code

pop_raster <- polygons |>
  as("SpatVector") |>
  rast(res = 100)

pycno <- pop_raster |>
  pycnophy(as(polygons, "SpatVector"), "population") |>
  terra::as.data.frame(xy = TRUE) |>
  rename(population_est = layer)

That’s it. Here’s a map of the result. Showing the boundaries of the regions makes clear how volume has migrated towards the centre of dense areas and away from the centre of low density areas.

Code

g2 <- ggplot() +
  geom_sf(data = context, colour = NA) +
  geom_raster(data = pycno, aes(x = x, y = y, fill = population_est)) +
  scale_fill_distiller("Smoothed population per ha.",
                       palette = "Purples", direction = 1) +
  geom_sf(data = polygons, fill = NA, lwd = 0.1) +
  coord_sf(expand = FALSE) +
  theme
g2

Figure 2: Pycnophylactic smoothed map of population

And we can confirm that the smoothing has preserved volume by summing the population from the source regions, and from the smoothed surface.

Code

polygons$population |> sum()

[1] 409941

Code

pycno$population_est |> sum()

[1] 409941

Hurrah! They are exactly equal.

There’s not a lot more to be said about this, other than perhaps to acknowledge yet another contribution by Waldo Tobler. I can’t help adding, in passing that in the opening paragraph of the paper he states that, “[A] common fact of geography is that places influence each other”,² a rather less snappy take on his eponymous ‘first law’.

Footnotes

Tobler WR. 1979. Smooth pycnophylactic interpolation for geographical regions. Journal of the American Statistical Association *74**(367) 519–530.↩︎
ibid. page 519↩︎