--- title: "`{pycnogrid}`: Flexible pycnophylactic interpolation to discrete global and local grid systems" author: - name: "Christopher D. Higgins" orcid: "0000-0002-3551-7750" email: "cd.higgins@utoronto.ca" affiliations: "University of Toronto" vignette: > %\VignetteIndexEntry{Getting started with `{pycnogrid}`} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} knitr: opts_chunk: collapse: true comment: '#>' bibliography: references.bib abstract: | Geographic data aggregated to administrative and statistical reporting units may not align with the geographic frame associated with particular research questions. Such units can also present the scale and zoning issues central to the modifiable areal unit problem (MAUP). In response, researchers and practitioners may turn to areal interpolation methods that can be used to transfer data from one geometry to another, with Tobler's pycnophylactic interpolation method distinguished by its emphasis on mass or volume preservation and spatial smoothing. Moreover, alternative hierarchical discrete global grid systems (DGGSs) that preserve area or shape or local grid systems present interesting opportunities for specifying alternative spatial support systems that have the potential to be more robust to the MAUP. Building on existing implementations of Tobler's pycnophylactic interpolation, the {pycnogrid} package offers a more flexible implementation of the algorithm that extends it beyond regular lattices to support a range of grid types, including the H3, A5, S2, and ISEA DGGSs as well as rasters and other local grids. Built on contemporary simple features and grounded in the foundations of spatial analysis, {pycnogrid} provides a tool for creating gridded, spatially smooth, and mass-preserving representations of aggregate extensive data for use in further analyses. --- ```{r} #| label: setup #| message: false #| include: false library(dplyr) library(sf) library(tmap) library(pycnogrid) ``` ## Introduction Spatial data capturing the socioeconomic, demographic, environmental, health, travel, and other characteristics of different places has become increasingly abundant. However, much of these data are released in a spatially aggregated form, represented through administrative and statistical reporting units such as census tracts, postal codes, or traffic analysis zones. While these spatial units facilitate data collection, privacy protection, and dissemination, they may not align with the geographic frame associated with particular research questions or analyses. This mismatch lies at the heart of the modifiable areal unit problem (MAUP) [@openshaw1979], which recognizes that observed spatial patterns and statistical relationships may vary according to the zoning system used to aggregate data. This includes the scale effect, related to the size of geographic units, and the zoning effect related to how the polygons are drawn. In response, researchers and practitioners may look to methods for transferring or interpolating information from source zones to more standardized target geographies while preserving known source zone totals. At the same time, the development of area or shape preserving hierarchical discrete global grid systems (DGGSs) and alternative local grid options presents an opportunity to extend traditional interpolation methods to accommodate different grid types. Among areal interpolation methods for extensive data, @tobler1979's pycnophylactic interpolation is distinguished by its emphasis on mass or volume preservation and spatial smoothing. Rather than assuming constant distributions of the source zone values across the target cells, pycnophylactic interpolation iteratively smooths the estimated target values while enforcing the preservation of the source zone totals in the system. With sufficiently fine target geometries, the result is a continuous spatial surface that preserves known aggregate quantities while reducing the abrupt discontinuities across zones introduced by administrative boundaries. In doing so, the method acknowledges that spatial aggregation obscures local variation and seeks to recover a *plausible* representation of the spatial phenomena at a finer alternative level of geography. The resulting surface embodies the intuition that nearby locations are expected to exhibit greater similarity than more distant locations [@tobler1970; @tobler2004]. Several implementations of Tobler's pycnophylactic interpolation exist for popular computing languages, including `{pycno}` [@pycno2025] for R and within the `{tobler}` package [@eliknaap2026] for Python, which is part of the larger Python Spatial Analysis Library (PySAL) [@sergiorey2026]. These existing implementations have largely followed Tobler's initial formulation of pycnophylactic interpolation from source zones to target raster cells. For example, the `pycno()` function from the `{pycno}` package converts source polygons to a regular grid and initializes a density surface by assigning a uniform density to all grid cells within each source zone. The density surface is then iteratively smoothed across neighbouring grid cells while correction steps ensure that the integrated target values within each source polygon remain equal to their original observed total. Built on the legacy `{sp}` framework, the package requires `SpatialPolygonsDataFrame` inputs and outputs a `SpatialGridDataFrame`. The PySAL `tobler` implementation follows a similar raster-based philosophy but performs smoothing through convolution-based operations on an array, with repeated correction steps used to preserve source-zone totals in the target cells. Building on these foundations, the primary contribution of `{pycnogrid}` is not a new pycnophylactic algorithm for users of the R computing language, but rather a more flexible implementation that extends pycnophylactic interpolation beyond regular raster lattices to a range of DGGS spatial grids. Capitalizing on the contemporary simple features `{sf}` [@pebesma2018] library and the neighbour relationships and spatial weights matrices that underpin much of modern spatial analysis and econometrics [@anselin1988] (through the `{spdep}` [@bivand2022] and `{sfdep}` [@sfdep2024] packages), `{pycnogrid}` implements pycnophylactic interpolation as a generalized neighbourhood-based smoothing problem. Rather than restricting smoothing to a regular raster lattice, values are smoothed using the neighbourhood structure defined by the underlying grid, allowing the method to be applied consistently across grid systems including H3, A5, S2, and various ISEA apertures, as well as traditional rasters and other local grids. This flexibility enables users to select a grid system whose geometric properties are better aligned with the requirements of subsequent analyses. ## Example Workflow The primary function in {pycnogrid} is `to_grid()`, which takes the following parameters: ```{{r}} #| label: to_grid #| echo: true pcynogrid::to_grid( source, value_col, id_col = NULL, grid_type = c("h3", "a5", "s2", "isea3h", "isea4h", "raster"), resolution, cell_inclusion = c("intersect", "centroid"), cell_allocation = c("area", "centroid"), nb_order = 1, max_iter = 500, tolerance = 1e-4, include_self = TRUE, missing_policy = c("abort", "warn", "ignore") ) ``` where - `source` is the source {sf} polygon layer containing the totals to be interpolated. Given the geometry calculations performed by the tool, only inputs with projected coordinate reference systems are accepted. - `value_col` is the column in `source` containing the count variable to be smoothed and preserved - `id_col` is an optional column uniquely identifying each source polygon, if omitted, an internal identifier is created - `grid_type` specifies the target grid system. Supported options are H3, A5, S2, ISEA grids with aperture-3 and 4, and raster-derived polygon cells - `resolution` controls the size of the target grid cells; its interpretation depends on the selected grid type - `cell_inclusion` defines how candidate grid cells are selected for interpolation. With `"intersect"`, cells are included if they intersect a source polygon. With `"centroid"`, cells are included only when their centroid falls inside a source polygon. - `cell_allocation` defines how source totals are allocated to grid cells. With `"area"`, values are allocated in proportion to the area of overlap between source polygons and grid cells. With `"centroid"`, each grid cell is assigned to the source polygon containing its centroid. - `nb_order` specifies the neighbourhood order used during smoothing. A value of 1 uses immediately adjacent cells, while larger values extend the smoothing neighbourhood outwards from a given cell. - `max_iter` sets the maximum number of smoothing iterations. If set to 0, the function returns the initial allocation without iterative smoothing. - `tolerance` defines the convergence threshold. Iteration stops when the relative change in estimated cell densities falls below this value. - `include_self` controls whether each cell includes its own current value when calculating the neighbourhood mean during smoothing. - `missing_policy` determines how the function handles source polygons that receive no target grid cells, which might arise due to a mismatch in source polygon sizes and target grid cell resolutions. "abort" stops with an error, "warn" returns a warning, and "ignore" proceeds silently. Together, these parameters define three main stages of the workflow: construction of the target grid, allocation of source totals to grid cells, and iterative pycnophylactic smoothing. The grid parameters determine the spatial support of the output surface, the allocation parameters determine the initial mass-preserving estimate, and the smoothing parameters control how values are redistributed across neighbouring cells while preserving the original source-zone totals. Additional functions, such as `to_h3()` and `to_a5()`, are grid-specific wrappers around `to_grid()`. The properties of the different grid systems are discussed further below. ### Source Zone Data The first part of a pycnophylactic workflow is to obtain more aggregate source data to be smoothed to a finer-resolution target grid. Pycnophylactic interpolation is designed for non-negative extensive variables -- quantities such as counts that can meaningfully be divided among smaller areas and whose total should be preserved. Due to the area and intersection calculations required for pycnophylactic smoothing, the source zones must be provided to the tool in a projected coordinate reference system. For an example case, 2020 population data was obtained for census tracts in New York City using the {tidycensus} [@tidycensus2026] package. A sub-sample of this data covering an area of Lower Manhattan is included in the package as `nyc_ct_small` (@fig-sample_population_count). The census tracts in this area contain a total of `r nyc_ct_small |> st_drop_geometry() |> summarize(sum(populationE)) |> prettyNum(big.mark = ",")` people. ```{r} #| label: fig-sample_population_count #| fig-cap: "Sample population count data for Lower Manhattan" #| warning: false #| echo: false tm_shape(nyc_ct_small) + tm_fill( fill = "populationE", fill.scale = tm_scale_intervals( values = "viridis", style = "jenks", n = 6, label.format = tm_label_format(big.num.abbr = c("K" = 3)) ), fill_alpha = .75, fill.legend = tm_legend(title = "Population Count", frame = FALSE) ) + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout( frame = FALSE, bg.color = "grey90" ) ``` ### Target Grid Systems The choice of a target grid system is not merely a computational consideration as different grid systems prioritize different geometric and analytical properties, including equal-area representation, shape preservation, neighbourhood structure, hierarchical indexing, and scalable spatial computation (@tbl-grids). ```{r} #| label: tbl-grids #| tbl-cap: "Supported target grid systems and their principal geometric and hierarchical properties" #| echo: false tibble::tribble( ~grid_system, ~cell_geometry, ~area_property, ~extent_and_hierarchy, ~subdivision_or_aperture, ~neighbour_topology, ~analytical_implication, "H3", "Mostly hexagons, with twelve pentagons", "Approximately equal-area", "Global discrete grid with a nested hierarchy", "Aperture-7 hierarchy with alternating cell orientation across resolutions", "Usually six edge-neighbours; pentagons have five", "Strong indexing and neighbourhood structure, particularly useful for mobility, accessibility, and spatial aggregation", "ISEA3H, ISEA4H", # ISEA7H", "Mostly hexagons, with twelve pentagons on an icosahedral projection", "Equal-area", "Global discrete grid with a nested hierarchy", "Aperture-3 and Aperture-4 hierarchies, #, and Aperture-7 hierarchies", "Usually six edge-neighbours; pentagons have five", "Equal-area hexagonal option with relatively fine-grained hierarchical scaling between resolutions", #"ISEA4H", #"Mostly hexagons, with twelve pentagons on an icosahedral projection", #"Equal-area", #"Global discrete grid with a nested hierarchy", #"Aperture-4 hierarchy", #"Usually six edge-neighbours; pentagons have five", #"Equal-area hexagonal option with a different resolution progression and nesting structure from aperture-3 systems", #"ISEA7H", #"Mostly hexagons, with twelve pentagons on an icosahedral projection", #"Equal-area", #"Global discrete grid with a nested hierarchy", #"Aperture-7 hierarchy", #"Usually six edge-neighbours; pentagons have five", #"Equal-area hexagonal option with resolution scaling broadly comparable to H3 but without H3's specific indexing system", #"ISEA4/3H", #"Mostly hexagons, with twelve pentagons on an icosahedral projection", #"Equal-area", #"Global discrete grid with a nested hierarchy", #"Mixed aperture-4 and aperture-3 hierarchy", #"Usually six edge-neighbours; pentagons have five", #"Allows a mixed resolution progression where the scale change between successive levels follows the selected aperture sequence", "A5", "Equal-area pentagonal cells", "Equal-area", "Global discrete grid with a hierarchical structure", "Five-way initial refinement, followed by four-way refinement", "Predominantly five edge-neighbours", "Provides a global equal-area alternative with strong indexing, although its pentagonal geometry may produce more directional variation than hexagonal grids", "S2", "Quadrilateral cells projected from the faces of a cube", "Not equal-area", "Global discrete grid with a nested hierarchy", "Quadtree subdivision: each cell has four children", "Variable topology across cube-face boundaries", "Strong global indexing and web-mapping infrastructure, but cell areas vary substantially across locations", "Local Raster", "Rectangular cells in a projected CRS", "Equal-area only when constructed in an appropriate projected CRS", "Usually regional or local; hierarchy is not intrinsic", "Resolution specified directly in map units", "Four- or eight-neighbour structure, depending on rook or queen contiguity", "Simple and familiar benchmark, but geometric properties depend on the selected projection and resolution", "Local Hex", "Regular hexagonal cells in a projected CRS", "Equal-area only when constructed in an appropriate projected CRS", "Usually regional or local; hierarchy is not intrinsic", "Resolution specified directly in map units", "Usually six edge-neighbours in a complete tessellation", "Local hexagonal support with six symmetric first-order neighbours, but geometric properties depend on the selected projection and resolution" ) |> gt::gt() |> gt::cols_label( grid_system = "Grid System", cell_geometry = "Cell Geometry", area_property = "Area Property", extent_and_hierarchy = "Extent and Hierarchy", subdivision_or_aperture = "Subdivision / Aperture", neighbour_topology = "Neighbour Topology", analytical_implication = "Analytical Implication" ) |> gt::tab_options( table.width = gt::pct(100), table.font.size = 12 ) ``` Sample discrete global grids generated at different resolution levels to cover the census tracts for Lower Manhattan from @fig-sample_population_count are shown in @fig-sample_dggs_grids. Among the DGGSs, the H3 grid system is a global-scale spatial indexing system developed by Uber to support routing and mobility analytics. It is a hierarchical hexagonal (mostly -- twelve pentagonal cells are required to accommodate the topology of a spherical surface) DGGS built on an icosahedral projection of the Earth with 16 resolution levels. H3 offers scalable hierarchical spatial indexing as each cell at a given level of the hierarchy can be sub-divided into a set of 7 child sells, although the cell tesselations rotate slightly at each resolution level. Compared to raster cells, the hexagonal tesselation also offers consistent neighbourhood relationships and more isotropic traversal pathways between neighbouring cells with approximately equal distances to all first-order neighbours. Because H3 is built on a gnomonic projection, hexagons tend to preserve their shape in projected spatial analytical workflows. However, because of this projection, H3 cells distort globally and are not equal area. While they may be approximately equal at the urban scale, cells compared at the global scale can differ in area quite significantly. Support for H3 in `{pycnogrid}` is offered through the `{h3o}` [@h3o2025] package. The Icosahedral Snyder Equal Area (ISEA) grid system partitions the globe into equal area hexagons with 31 supported resolution levels. Similar to H3, twelve pentagonal cells are required. In contrast to H3, where cells are approximately equal area, ISEA grids prioritize identical cell areas with more aligned nesting. This comes at the cost of shape preservation, with ISEA hexagons appearing more elongated than H3 when represented in conventional geographic or local projected coordinate systems. The ISEA grid supports several different apertures, including Aperture-3, 4, 7, and 4/3 mixed. At Aperture-3, parent cells divide into 3 child cells, while at Aperture-7 the ISEA grid behaves similar to H3 with each cell splitting into 7 across resolution levels. Support for ISEA grids is currently limited to Aperture-3 and 4 and offered through the `{hexify}` [@hexify2026] package. The A5 grid system is a recent DGGS that is explicitly designed to be equal area across the globe. Based on pentagons and 31 different resolution levels, A5 cells at a given resolution sub-divide into 5 child cells. The primary motivation behind A5 is spatial analysis as the equal area properties of the grid system greatly simplify comparisons of densities, rates, and other aggregated quantities over space. However, the subdivision of parent into child cells is not regular with cell shapes becoming more elongated at higher resolutions, meaning the distance between cell centroids is not uniform. With these slightly elongated pentagonal shapes, A5 cells are more anisotropic than H3 cells, although the extent to which impacts traversal pathways across the equal area cell topology has not yet studied. Support for the A5 grid system is offered through the `{a5R}` [@a5R2026] package. S2 was developed by Google for global-scale spatial indexing and geometric computation. S2 begins by projecting the Earth onto the flat planes of a cube. Within a 31-level resolution hierarchy, each cell sub-divides into 4 child cells. While use of H3 and A5 might be motivated by neighbourhood regularity and equal-area representation, respectively, S2 emphasizes hierarchical spatial indexing and efficient spherical geometry operations. Consequently, the use of S2 within {pycnogrid} may be attractive when interpolated data are intended for integration with large-scale geospatial databases, distributed computing environments, or workflows that already rely on S2 indexing. Support for S2 in `{pycnogrid}` is offered through the `{s2}` [@s22025] package. ```{r} #| label: fig-sample_dggs_grids #| fig-cap: "Sample DGGS hierarchies at different resolutions" #| fig-height: 6 #| fig-width: 9 #| warning: false #| echo: false sample_dggs_grids <- bind_rows( h3_grids <- purrr::map_dfr( list(8, 9, 10), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "h3", resolution = .x, cell_inclusion = "intersect" ) |> mutate(grid_type = "h3", cell_resolution = .x, zorder = cell_resolution) ), isea3h_grids <- purrr::map_dfr( list(17, 18, 19), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "isea3h", resolution = .x, cell_inclusion = "intersect" ) |> mutate(grid_type = "isea3h", cell_resolution = .x, zorder = cell_resolution) ), isea4h_grids <- purrr::map_dfr( list(13, 14, 15), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "isea4h", resolution = .x, cell_inclusion = "intersect" ) |> mutate(grid_type = "isea4h", cell_resolution = .x, zorder = cell_resolution) ), a5_grids <- purrr::map_dfr( list(13, 14, 15), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "a5", resolution = .x, cell_inclusion = "intersect" ) |> mutate(grid_type = "a5", cell_resolution = .x, zorder = cell_resolution) ), s2_grids <- purrr::map_dfr( list(14, 15, 16), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "s2", resolution = .x, cell_inclusion = "intersect" ) |> mutate(grid_type = "s2", cell_resolution = .x, zorder = cell_resolution) ) ) |> mutate(grid_type = forcats::fct_relevel(grid_type, "h3", "isea3h", "isea4h", "a5", "s2")) |> arrange(desc(zorder)) tm_shape(sample_dggs_grids) + tm_lines( col = "cell_resolution", col.scale = tm_scale_categorical(), col.free = TRUE, col.legend = tm_legend(title = "resolution", frame = FALSE) ) + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_facets(by = "grid_type", ncol = 3, free.coords = FALSE) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout( panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90" ) ``` Of the supported local grid types in @fig-sample_local_grids, raster grids are the most familiar, representing geographic space as a regular lattice of equally sized cells. Their simplicity has made them the dominant representation for continuous spatial phenomena and the traditional support for pycnophylactic interpolation. However, raster grids are inherently planar, which requires projection choices, and neighbourhood relationships depend on user-specified contiguity definitions (e.g., rook or queen adjacency). In terms of indexing, raster cells are not uniquely identified by a global index and any hierarchical relationships among cells of different resolutions would have to be manually specified. The neighbourhood structure is also anisotropic, with greater distances required to traverse cells diagonally than horizontally or vertically. Support for rasters in `{pycnogrid}` is offered through the `{terra}` [@terra2026] package. Finally, users can also generate local regular hexagonal tessellations. Unlike the global DGGS options, these grids are defined within the coordinate reference system of the input data and have no intrinsic global index or hierarchical structure. Their area and shape properties therefore depend on selecting an appropriate projected CRS, while their alignment and resolution are specified directly in the map units of that CRS. Built on the {sf} [@pebesma2018] package, this option provides a familiar local hexagonal support for analyses focused on a particular region or study area. ```{r} #| label: fig-sample_local_grids #| fig-cap: "Sample local grid hierarchies at different resolutions" #| fig-height: 3 #| fig-width: 9 #| message: false #| echo: false sample_local_grids <- dplyr::bind_rows( raster_grids <- purrr::map_dfr( list(1000, 500, 250), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "raster", resolution = .x, cell_inclusion = "intersect" ) |> mutate( grid_type = "raster", cell_resolution = .x, zorder = 1 / cell_resolution ) ), hex_grids <- purrr::map_dfr( list(1000, 500, 250), ~ pycnogrid:::prepare_target_cells( source_polys = nyc_ct_small, grid_type = "hex", resolution = .x, cell_inclusion = "intersect" ) |> mutate( grid_type = "hex", cell_resolution = .x, zorder = 1 / cell_resolution ) ) ) |> mutate(grid_type = forcats::fct_relevel(grid_type, "raster", "hex")) |> arrange(desc(zorder)) tm_shape(sample_local_grids) + tm_lines( col = "cell_resolution", col.scale = tm_scale_categorical(), col.free = TRUE, col.legend = tm_legend(title = "resolution", frame = FALSE) ) + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_facets_wrap(by = "grid_type", #ncol = 2, free.coords = FALSE) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout( panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90" ) ``` ### Pycnophylactic Interpolation With source data collected, the pycnophylactic smoothing can now be run. In this example case, an H3 grid at a resolution of 10 is used: ```{r} #| label: run to_grid pycno_nyc_ct_small <- nyc_ct_small |> pycnogrid::to_grid( value_col = populationE, grid_type = "h3", resolution = 10 ) ``` The results of this interpolation are shown in @fig-pycno_nyc_ct_small below: ```{r} #| label: fig-pycno_nyc_ct_small #| fig-cap: "Census tract population counts interpolated to an H3 grid" #| echo: false tm_shape(pycno_nyc_ct_small) + tm_fill( fill = "pycno_populationE", fill.scale = tm_scale_intervals( values = "viridis", style = "jenks", n = 6, label.format = tm_label_format(big.num.abbr = c("K" = 3)) ), fill_alpha = .75, fill.legend = tm_legend(title = "Interpolated \nPopulation", frame = FALSE) ) + tm_shape(nyc_ct_small) + tm_lines(col = "grey45") + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout( frame = FALSE, bg.color = "grey90" ) ``` The results for other grid types can be seen in @fig-pycno_nyc_ct_small_resolutions: ```{r} #| label: fig-pycno_nyc_ct_small_resolutions #| fig-cap: "Census tract population counts interpolated to all four grid types" #| fig-height: 8 #| fig-width: 10 #| echo: false grid_list <- list("h3", "isea3h", "isea4h", "a5", "s2", "raster", "hex") resolution_list <- list(10, 20, 16, 16, 16, 100, 100) pycno_nyc_ct_small_resolutions <- purrr::map2_dfr( grid_list, resolution_list, ~ to_grid( source = nyc_ct_small, value_col = "populationE", grid_type = .x, resolution = .y ) |> mutate(grid_type = .x, resolution = .y, facet_label = paste0(grid_type, " at resolution ", resolution)) ) tm_shape(pycno_nyc_ct_small_resolutions |> filter(grid_type != "h3" & grid_type != "isea4h" & grid_type != "hex")) + tm_fill( fill = "pycno_populationE", fill.scale = tm_scale_intervals( values = "viridis", style = "jenks", n = 6, label.format = tm_label_format(big.num.abbr = c("K" = 3)) ), fill_alpha = .75, fill.legend = tm_legend(title = "Interpolated \nPopulation", frame = FALSE), fill.free = TRUE ) + tm_facets(by = "facet_label", ncol = 2, free.coords = TRUE) + tm_shape(nyc_ct_small) + tm_lines(col = "grey45") + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout(panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90") ``` The summary statistics of the population variable for the different grid specifications is shown in @tbl-grid_descriptives. While the different grid types and resolutions necessarily produce different mean population values across the target cells, the total population spread over the target cells is intact. ```{r} #| label: tbl-grid_descriptives #| tbl-cap: "Population descriptive statistics for different grid types" #| echo: false pycno_nyc_ct_small_resolutions |> st_drop_geometry() |> group_by(grid_type) |> summarize( resolution = min(resolution), cell_count = n(), mean_pop = mean(pycno_populationE), sum_pop = sum(pycno_populationE) ) |> ungroup() |> mutate(grid_type = forcats::fct_relevel(grid_type, "h3", "isea3h", "isea4h", "a5", "s2", "raster", "hex")) |> arrange(grid_type) |> gt::gt() |> gt::cols_label( grid_type = "grid type", resolution = "grid resolution", cell_count = "grid cell count", mean_pop = "mean cell population", sum_pop = "total population" ) ``` The main customizations for `to_grid()` involve defining additional options to guide the pycnophylactic interpolation. First, the cell inclusion and allocation criteria impact how target cells are generated and how source values are allocated to them. - with `cell_inclusion = "intersect"` and `cell_allocation = "area"`, all cells intersecting the source polygons are retained, and source totals are allocated according to the area of overlap between source polygons and target cells. This is the most geographically complete and resolution robust representation of the source layer as the data from all source-target intersections are retained. Partially covered edge cells receive only the share of the source value associated with their covered area, avoiding the extrapolation of counts beyond source boundaries. These are the default settings. - with `cell_inclusion = "intersect"` and `cell_allocation = "centroid"`, all intersecting grid cells remain in the target grid, but source values are assigned only to cells whose centroids fall within a source polygon. Target cells whose centroids fall outside the source polygons do not receive any initial allocation of the source values, although the smoothing process may allocate some values to these cells. Compared to area-based allocation, this simpler allocation method is more sensitive to the grid resolution and the placement of target cell centroids. Values associated with any source polygons that have no assigned target cells, such as small or narrow polygons, will be omitted from the system. Moreover, by representing a source value over the full extent of a target cell, this allocation method can extend or extrapolate source values beyond the original source polygon boundaries. - with `cell_inclusion = "centroid"` and `cell_allocation = "area"`, only target cells whose centroids fall within the source layer are retained, but the retained cells receive source values according to their actual areas of overlap. While this avoids retaining cells that merely touch a source boundary, the resulting grid may not fully cover the source polygons, leading to a truncation in the spatial support offered by the target cells. The area-based allocation nevertheless helps to preserve source values in the system based on the areal overlap of source polygons and target cells. - with `cell_inclusion = "centroid"` and `cell_allocation = "centroid"`, target cells are both selected and assigned according to centroid location. Each retained cell is associated with the source polygon containing its centroid, producing the simplest and most discrete source-to-grid assignment with the strongest concentrations of source values within the target cells. However, this approach is the most sensitive to grid resolution and alignment relative to the source polygons, risking truncated spatial support of the source layers within the target cell tessellation, the extrapolation of values beyond the source extent, and the potential loss of source values if no cell centroids fall within the original polygon boundaries. The smoothed results using the default and other cell inclusion and allocation settings are shown in @fig-pycno_nyc_ct_small_combinations below. ```{r} #| label: fig-pycno_nyc_ct_small_combinations #| fig-cap: "Interpolated population counts with varying inclusion and allocation parameters" #| fig-height: 8 #| fig-width: 8 #| echo: false combinations <- tidyr::expand_grid( inclusion_list = list("intersect", "centroid"), allocation_list = list("area", "centroid") ) pycno_nyc_ct_small_combinations <- purrr::map2_dfr( combinations |> purrr::pluck("inclusion_list"), combinations |> purrr::pluck("allocation_list"), ~ to_grid( source = nyc_ct_small, value_col = "populationE", id_col = id, grid_type = "h3", resolution = 10, cell_inclusion = .x, cell_allocation = .y, nb_order = 1 ) |> mutate( inclusion = paste0("inclusion: ", .x), allocation = paste0("allocation: ", .y) ) ) |> mutate(inclusion = forcats::fct_relevel(inclusion, "inclusion: intersect", "inclusion: centroid")) tm_shape(pycno_nyc_ct_small_combinations) + tm_fill( fill = "pycno_populationE", fill.scale = tm_scale_intervals( values = "viridis", style = "jenks", n = 6, label.format = tm_label_format(big.num.abbr = c("K" = 3)) ), fill_alpha = .75, fill.legend = tm_legend(title = "Interpolated \nPopulation", frame = FALSE) ) + tm_facets_grid(rows = "inclusion", columns = "allocation", free.coords = FALSE) + tm_shape(nyc_ct_small) + tm_lines(col = "grey45") + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout(panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90") ``` The second customization involves setting the number of neighbours through the `nb_order` parameter. Here the default is `1`, which results in the spatial weights matrix accounting for spatial relationships between a given target cell and its first-order neighbours with Queen contiguity. Increasing the order of the neighbourhood includes target cells from a larger neighbourhood which this has the effect of further smoothing out the spatial patterns of the interpolated values in the target cells. This can be seen in @fig-pycno_nyc_ct_small_nb_order. ```{r} #| label: fig-pycno_nyc_ct_small_nb_order #| fig-cap: "Interpolated population counts with increasing neighbourhood order" #| fig-height: 8 #| fig-width: 8 #| echo: false neighbours_list <- seq(1, 4, by = 1) pycno_nyc_ct_small_neighbourss <- purrr::map_dfr( neighbours_list, ~ to_grid( source = nyc_ct_small, value_col = "populationE", grid_type = "h3", resolution = 10, nb_order = .x ) |> mutate(nb_order = .x, facet_label = paste0("nb_order = ", nb_order)) ) tm_shape(pycno_nyc_ct_small_neighbourss) + tm_fill( fill = "pycno_populationE", fill.scale = tm_scale_intervals( values = "viridis", style = "jenks", n = 6, label.format = tm_label_format(big.num.abbr = c("K" = 3)) ), fill_alpha = .75, fill.legend = tm_legend(title = "Interpolated \nPopulation", frame = FALSE), fill.free = FALSE ) + tm_facets(by = "facet_label", ncol = 2, free.coords = TRUE) + tm_shape(nyc_ct_small) + tm_lines(col = "grey45") + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout(panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90") ``` ## Detailed Example To explore how the smoothing algorithm works, the code below selects a subset of the `nyc_ct_small` census tracts and generates an H3 grid at `resolution = 9`. To limit the number of cells generated for this example, the `cell_inclusion` criteria is set to `"centroid"` so the target grid consists only of cells whose centroids fall within the source polygon geometries (@fig-nyc_ct_subset): ```{r} #| label: fig-nyc_ct_subset #| fig-cap: "Source zone subset and target cells with centroid-based inclusion criteria" #| echo: false subset_ids <- c( "36061005800", "36061008700", "36061009100" ) nyc_ct_subset <- nyc_ct |> filter(id %in% subset_ids) |> mutate(source_id = paste0("S", row_number())) prepped_source <- pycnogrid:::prep_source(source = nyc_ct_subset, value_col = "populationE", id_col = "id") source_polys <- prepped_source |> purrr::pluck("source") |> mutate(source_label_id = paste0("S", .sid)) source_values <- prepped_source |> purrr::pluck("source_values") input_total_original <- prepped_source |> purrr::pluck("input_total_original") input_total_represented <- prepped_source |> purrr::pluck("input_total_represented") target_cells <- pycnogrid:::prepare_target_cells( source_polys, grid_type = "h3", resolution = 9, cell_inclusion = "centroid" #cell_inclusion = "intersect" ) |> mutate(target_label_id = paste0("T", .tid)) tm_shape(source_polys) + tm_fill(col = "grey35", col_alpha = .8, fill = "grey", fill_alpha = .1) + tm_text("source_label_id", ymod = -1, angle = 330) + tm_shape(target_cells) + tm_fill(col = "blue", fill = "blue", fill_alpha = .1) + tm_text("target_label_id", col = "blue") + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout(panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90") ``` Within this scenario, the three source zones are indexed as $i = 1, \ldots, n = 3$ and have the following observed populations: $$ \mathbf{y} = \begin{bmatrix} 6445 \\ 4412 \\ 5712 \end{bmatrix} $$ with a total population $\sum_i y_i = 16569$. ### Areal Interpolation These source population values will be interpolated to the five target H3 grid cells indexed as $j = 1, \ldots, m = 5$.. The first step is to determine the spatial relationships between the source zones and target cells. Using the default `cell_allocation` based on `"area"`, this relationship is represented by the calculation of an areal overlap matrix $\mathbf{A}$, where each element $a_{ij}$ is the area of intersection between a source zone $i$ and target cell $j$: $$ a_{ij} = \mathrm{Area}(S_i \cap T_j) $$ For the example source zones and target cells, these overlap areas (in $m^2$) are: ```{r} #| echo: false #| warning: false # manually: inter_raw <- sf::st_intersection(source_polys, target_cells) |> dplyr::mutate(.weight = as.numeric(sf::st_area(geometry))) |> sf::st_drop_geometry() |> dplyr::filter(.weight > 0) |> dplyr::select(source_id, .tid, .weight) represented_sources <- inter_raw |> dplyr::distinct(source_id) source_index <- source_values |> dplyr::mutate(.sid = dplyr::row_number()) |> dplyr::select(source_id, .sid) inter <- inter_raw |> # semi join to drop any target cells that don't have any sources dplyr::semi_join(source_values, by = "source_id") |> # attach the numeric source index .sid. dplyr::left_join(source_index, by = "source_id") |> # Sum weights in case the same source-cell pair appears more than once, # which can happen after geometric operations involving multipart polygons # or multiple intersection fragments. dplyr::group_by(.sid, .tid, source_id) |> dplyr::summarise(.weight = sum(.weight, na.rm = TRUE), .groups = "drop") n_target <- nrow(target_cells) n_source <- nrow(source_values) A_matrix <- Matrix::sparseMatrix( i = inter$.sid, j = inter$.tid, x = inter$.weight, dims = c(n_source, n_target) ) # row and column normalize the A_matrix overlap matrix # need some sparse matrix trickery here for row and column sums row_sums <- Matrix::rowSums(A_matrix) W_matrix <- Matrix::Diagonal(x = 1 / pmax(row_sums, .Machine$double.eps), n = n_source) %*% A_matrix col_sums <- as.numeric(Matrix::colSums(A_matrix)) V_matrix <- A_matrix %*% Matrix::Diagonal(x = 1 / pmax(col_sums, .Machine$double.eps), n = n_target) # cross-reference with package function: source_cell_weights <- pycnogrid:::prepare_source_cell_weights( source_polys, source_values, target_cells, grid_type = "h3", resolution = 9, cell_allocation = "area", missing_policy = "warn", input_total_original ) #A_matrix == source_cell_weights$A_matrix #W_matrix == source_cell_weights$W_matrix #V_matrix == source_cell_weights$V_matrix ``` $$ \mathbf{A}= \begin{bmatrix} 98741.439 & 5558.235 & 2253.077 & 1419.009 & 36725.10 \\ 1506.071 & 77478.171 & 83178.843 & 0 & 0 \\ 0 & 3919.766 & 0 & 95050.580 & 38363.57 \end{bmatrix} $$ The rows of $\mathbf{A}$ sum to the total area of each source zone covered by target cells within the source-target intersection, that is $c_i = \sum_j a_{ij}$. With some grid cells not perfectly overlapping the source zones due to the use of centroid-based cell inclusion criteria in this example, the total covered area $c_i$ is less than the total area of a zone as defined by its geometry. Similarly, the column sums of $\mathbf{A}$, calculated as $c_j = \sum_i a_{ij}$, reflect the covered area of each target cell, or the portion of each target cell lying within the source polygons. These column sums will inform the construction of a column-normalized matrix $\mathbf{V}$ below. Next, the overlap areas for each source zone are converted to areal interpolation weights by row-standardization: $$ w_{ij} = \frac{a_{ij}}{c_i} $$ The resulting source zone to target cell weight matrix is: $$ \mathbf{W}= \begin{bmatrix} 0.682 & 0.038 & 0.016 & 0.010 & 0.254 \\ 0.009 & 0.478 & 0.513 & 0 & 0 \\ 0 & 0.029 & 0 & 0.692 & 0.279 \end{bmatrix} $$ Because each row of $\mathbf{W}$ sums to one (e.g. $\sum_j w_{ij}=1$), these weights ensure that the entire total from each source zone is distributed across the target cells while preserving the original total of the source values in the system. The initial target cell estimates are obtained by applying these weights to the source values in $\mathbf{y}$: ```{r} #| echo: false y <- source_index |> dplyr::left_join(source_values, by = "source_id") |> dplyr::arrange(.sid) |> dplyr::pull(source_value) # Calculate initial areal allocation x <- as.numeric(Matrix::t(W_matrix) %*% y) ``` $$ \mathbf{x}^{(0)} = \mathbf{W}^{\top}\mathbf{y} $$ which yields the areally-interpolated value estimates for each of the five target cells: $$ \mathbf{x}^{(0)} = \begin{bmatrix} 4439.057 \\ 2518.565 \\ 2363.417 \\ 4016.554 \\ 3231.407 \end{bmatrix} $$ From this, the source zone value total allocated within the system remains unchanged: $$ \sum_i^3 y_i = 16569 = \sum_j^5 x_j^{(0)} $$ Because pycnophylactic interpolation works with densities, these allocated totals are converted to density estimates by dividing each target cell's allocated population by its covered area in the source zone-target cell intersection $c_j$ (which are the column sums of $A$): $$ d_j^{(0)} = \frac{x_j^{(0)}}{c_j} $$ Resulting in the vector of initial target cell densities: ```{r} #| echo: false density <- x / col_sums ``` $$ \mathbf d^{(0)} = \begin{bmatrix} 0.0443 \\ 0.0290 \\ 0.0277 \\ 0.0416 \\ 0.0430 \end{bmatrix} $$ These densities represent the estimated population density over the portion of each target cell covered by the source polygons. Using area-based cell allocation, density is calculated for cells that extend beyond the source-zone boundaries using only the covered area represented in the source-target intersection. By extension, when a source zone is not fully covered by target cells, as is the case when using centroid-based cell inclusion criteria in this analysis, the incomplete spatial support results in the entire source value total for a zone being distributed across the cells that intersect it. Consequently, the procedure preserves source totals within the represented spatial support, even when that support differs from the original source geometry. ### Spatial Smoothing ```{r} #| echo: false S_matrix <- pycnogrid:::build_smoothing_matrix( target_cells = target_cells, nb_order = 1, include_self = TRUE ) ``` Once the initial density estimates have been obtained for the target cells, pycnophylactic smoothing proceeds iteratively. The first step is to define the neighbourhood relationships among target cells. Using a typical first-order Queen contiguity for the five target cells in this example, the neighbourhood structure is represented by the row-standardized sparse spatial weights matrix: $$ \mathbf{S} = \begin{bmatrix} 0.20 & 0.20 & 0.20 & 0.20 & 0.20 \\ 0.25 & 0.25 & 0.25 & 0.25 & . \\ 0.33 & 0.33 & 0.33 & . & . \\ 0.25 & 0.25 & . & 0.25 & 0.25 \\ 0.33 & . & . & 0.33 & 0.33 \end{bmatrix} $$ Each row in $\mathbf{S}$ describes the cells contributing to the smoothed density estimate of a target cell. Because the matrix is row-standardized, each row sums to one and the smoothing operation computes a local weighted average of neighbouring densities. The first smoothing step ($t+1$) involves multiplying the initial density vector $\mathbf{d^{(t)}}$ by the spatial weights matrix $\mathbf{S}$: ```{r} #| echo: false d_smooth <- as.numeric(S_matrix %*% density) ``` $$ \tilde{\mathbf{d}}^{(t + 1)} = \mathbf{S} \mathbf{d}^{(0)} $$ Here, the $\tilde{\mathbf{d}}$ denotes these are the smoothed density estimates before correction. This operation reduces abrupt differences between neighbouring cells by replacing each cell’s density with the average density of its neighbourhood. For target zone 1: $$ \tilde{d}_{1}^{(t+1)} = 0.20(0.0443) + 0.20(0.0290) + 0.20(0.0277) + 0.20(0.0416) + 0.20(0.0430) $$ $$ \tilde{d}_{1}^{(t+1)} = 0.0371 $$ However, smoothing alone does not always preserve the original source-zone totals. After smoothing, the implied source-zone total estimates are calculated by multiplying the 3 by 5 overlap matrix (which contains information on how much of each target cell's area is within the system defined by source zone-target cell spatial relationship) and the 5 by 1 column vector of smoothed densities: ```{r} #| echo: false source_est <- as.numeric(A_matrix %*% density) ``` $$ \tilde{\mathbf{y}}^{(t + 1)} = \mathbf{A} \tilde{\mathbf{d}}^{(t + 1)} $$ The resulting source values contained within the 3 by 1 column vector $\tilde{\mathbf{y}}^{(t + 1)}$ generally differ from their observed source totals ($\mathbf{y}$). To restore the original source-zone populations, a correction factor is calculated for each source zone: ```{r} #| echo: false r_ratio <- y / source_est ``` $$ \mathbf{r} = \frac{\mathbf{y}}{\tilde{\mathbf{y}}^{(t + 1)}} $$ Resulting in this column vector of ratios: $$ \mathbf{r}^{(t + 1)} = \begin{bmatrix} 1.034 \\ 0.957 \\ 0.998 \end{bmatrix} $$ The correction factor compares the observed source total to the total implied by the smoothed density surface. A value of $r_{i} = 1$ indicates that smoothing preserved the source-zone total exactly, so no adjustment is required. A value of $r_i > 1$ indicates that the smoothed density surface underestimates the observed population in source zone i, requiring the densities within that source zone to be increased. Conversely, a value of $r_i < 1$ indicates that smoothing overestimates the source-zone population, requiring the densities to be decreased. For example, $r_1^{(t+1)} = 1.034$ indicates that the smoothed density surface underestimates the population of source zone 1 by approximately 3.4% when aggregated back to the source geometry. Likewise, $r_2^{(t+1)} = 0.957$ indicates that the smoothed surface overestimates the population of source zone 2 by approximately 4.3%. These source-zone correction factors are then redistributed back to target cells according to their overlap relationships using a column standardized weights matrix $\mathbf{V}$, with $v_{ij} = \frac{a_{ij}}{c_j}$: ```{r} #| echo: false q_correction <- as.numeric(Matrix::t(V_matrix) %*% r_ratio) ``` $$ \mathbf{q}^{(t + 1)} = \mathbf{V}^\top \mathbf{r}^{(t + 1)} $$ From this, each target cell receives a weighted average of the correction factors from the source zones that it intersects. These correction factors can now be applied to adjust the smoothed densities for each target cell for the first iteration: ```{r} #| echo: false density_t_plus_one <- density * q_correction ``` $$ \mathbf{d}^{(t+1)} = \tilde{\mathbf d}^{(t+1)} \odot \mathbf q^{(t+1)} $$ Or, using scalar notation: $$ d_{j}^{(t+1)} = \tilde{d}_{j}^{(t+1)} q_{j}^{(t+1)} $$ In this way, pycnophylactic smoothing uses the original source zone polygons to ensure the total mass is preserved on the smoothed surface. The corrected densities form the starting point for the next iteration $t+2$. Smoothing and correction continue until the relative change in target-cell densities between successive iterations falls below a specified tolerance. In this case, the convergence criterion is defined by calculating the mean relative change in target cell densities across the current and previous iterations: ```{r} #| echo: false abs_change <- abs(density_t_plus_one - density) denom <- pmax(abs(density), 1e-12) last_error <- mean(abs_change / denom, na.rm = TRUE) ``` $$ \epsilon^{(t+1)} = \frac{1}{m} \sum_{j=1}^{m} \frac{\left| d_j^{(t+1)} - d_j^{(t)} \right|} {\max\!\left(d_j^{(t)},\delta\right)} $$ with $\delta=10^{-12}$ used to prevent any division by zero in the denominator. Iteration stops when $\epsilon^{(t+1)}$ is less than the convergence tolerance, which is set at a default of $10^{-4}$. For the present case, the error from the first iteration is `r round(last_error, 4)`, suggesting further rounds of smoothing are required to reach convergence. The final results of the interpolation are shown in @fig-pycno_nyc_ct_subset: ```{r} #| label: fig-pycno_nyc_ct_subset #| fig-cap: "Source zone and interpolated population values" #| fig-width: 8 #| fig-height: 4 #| echo: false subset_results <- dplyr::bind_rows(nyc_ct_subset |> to_grid( value_col = "populationE", id_col = id, grid_type = "h3", resolution = 9, #grid_type = "raster", #resolution = 275, cell_inclusion = "centroid", cell_allocation = "area", nb_order = 1 ) |> mutate( label_id = paste0("T", .tid), population = pycno_populationE, type = "target cells" ), source_polys |> mutate( population = .source_value, label_id = paste0("S", .sid), type = "source zones") ) subset_overlay <- dplyr::bind_rows( source_polys |> mutate( type = "target cells"), target_cells |> mutate(type = "source zones") ) tm_shape(subset_results) + tm_fill( fill = "population", fill.scale = tm_scale_ordinal( values = "viridis", label.format = tm_label_format(big.num.abbr = c("K" = 3)) ), fill_alpha = .75, fill.legend = tm_legend(title = "Population", frame = FALSE, position = tm_pos_in()), fill.free = TRUE ) + #tm_text("label_id") + tm_facets(by = "type", free.coords = FALSE) + tm_shape(subset_overlay) + tm_lines(col = "grey35", col_alpha = .8) + tm_facets(by = "type", free.coords = FALSE) + #tm_basemap("CartoDB.PositronNoLabels", alpha = .5) + tm_shape(nyc_ct) + tm_fill(fill = "grey97", zindex = 3) + tm_layout(panel.label.frame = FALSE, panel.label.bg = FALSE, frame = FALSE, bg.color = "grey90") ``` ## Conclusion `{pycnogrid}` offers a flexible implementation of Tobler's pycnophylactic interpolation for redistributing polygon-based count data to discrete global and raster grid systems. Building on existing implementations of the concept, the smoothing algorithm that underpins the package expresses the interpolation through the lens of neighbourhoods and spatial weights, enabling the pycnophylactic approach to be applied to a range of different grid systems. At present, these grid systems include H3, A5, S2, and ISEA grids, as well as traditional raster arrays. While the preservation of mass or volume inherent in the source count data is central to the pycnophylactic smoothing approach, a core aspect of `{pycnogrid}` is the option to allocate and preserve mass while accounting for source zone and target cell areal intersections. The package also makes explicit several analytical choices that are often implicit in interpolation workflows. Users may select how target cells are included, how source values are allocated to cells, the neighbourhood order used for smoothing, and whether cells contribute to their own neighbourhood mean. These choices affect the implied spatial support and degree of smoothing, and should be selected with reference to the intended downstream application rather than treated as purely technical defaults. While pycnophylactic interpolation produces a smooth, mass-preserving representation that is consistent with the available source totals and the chosen target-grid structure, the outputs are best understood as analytically-useful plausible estimates rather than as observed values at a finer level of geography. Future work could extend `{pycnogrid}` through ancillary-data weighting [e.g. @tapp2010], alternative smoothing operators, and support for additional area or shape preserving DGGSs. Systematic evaluations of how grid geometry, resolution, hierarchy, and smoothing assumptions affect the robustness of downstream spatial analyses to the MAUP would also be valuable. To that end, by bringing pycnophylactic interpolation to contemporary vector-based spatial workflows and multiple grid systems, `{pycnogrid}` provides a practical foundation for creating comparable, mass-preserving gridded representations of aggregated spatial data. ## References