Estimation of Population from Pycnophylatic Interpolation | R Documentation |
Given a SpatialGridDataFrame
of population estimates and a set of polygons,
compute a population estimate based on Tobler's pycnophylactic interpolation
algorithm for each zone. The result is a vector.
estimate.pycno(sgdf, spdf)
sgdf |
A |
spdf |
A |
Takes the estimate of population density for each pixel, checks which polygon each pixel is in, and aggregates them. Accuracy depends on the scale of pixels in the initial interpolation.
A vector in which each each pixel set at the estimated population aggregation to each zone in spdf
.
Pycnophylatic interpolation has the property that the sum of the estimated values associated with all of the pixels in any polygon equals the supplied population for that polygon. A further property is that all pixel values are greater than or equal to zero. The method is generally used to obtain pixel-based population estimates when total populations for a set of irregular polygons (eg. counties) are known.
Chris Brunsdon
Tobler, W.R. (1979) Smooth Pycnophylactic Interpolation for Geographical Regions. Journal of the American Statistical Association, v74(367) pp. 519-530.
pycno
library(sp)
# Read in data for North Carolina as a SpatialPolygonsDataFrame
#nc.sids <- readShapeSpatial(system.file("shapes/sids.shp", package="maptools")[1],
# IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
nc.sids <- as(sf::st_read(system.file("shape/nc.shp", package="sf")), "Spatial")
row.names(nc.sids) <- as.character(nc.sids$FIPSNO)
# Compute the pycnophylactic surface for 1974 births as a SpatialGridDataFrame
# Note probably shouldn't really base grid cells on Lat/Long coordinates
# This example just serves to illustrate the use of the functions
births74 <- pycno(nc.sids,nc.sids$BIR74,0.05,converge=1)
# Create a new 'blocky' set of zones
#blocks <- gUnionCascaded(nc.sids,1*(coordinates(nc.sids)[,2] > 36) +
# 2*(coordinates(nc.sids)[,1] > -80))
crds <- sf::st_coordinates(sf::st_centroid(sf::st_geometry(sf::st_as_sf(nc.sids)),
of_largest_polygon = TRUE))
block_ID <- 1*(crds[,2] > 36) + 2*(crds[,1] > -80)
temp <- sf::st_as_sf(nc.sids)
temp$block_ID <- block_ID
blocks <- as(aggregate(temp, by=list(temp$block_ID), head, n=1), "Spatial")
# Plot the blocky zones
plot(blocks)
# Aggregate data to them
estimates <- estimate.pycno(births74,blocks)
# Write the estimates on to the map
text(coordinates(blocks),as.character(estimates))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.