KDE for Population Modelling

Notes on using Kernel Density Estimation for modelling population density

Packages

First, I am loading all packages necessary for the following analysis.

library(raster)# deal with raster
library(sf) # spatial class
library(mapview) # interactive map viewing
library(tmap) # cartography
library(tmaptools) # little helpers
library(dplyr) # manipulate
library(esri2sf) # get data from the server
library(spatstat) # spatial statistics
library(sp) 
library(rgdal)
library(maptools)
library(spex)

Data

Data Zones Centroids

url_dz =  "http://sedsh127.sedsh.gov.uk/arcgis/rest/services/ScotGov/StatisticalUnits/MapServer/4"

dz_cntr = esri2sf(url_dz, outFields=c("TotPop2011")) %>% st_transform(27700)

## [1] "Feature Layer"
## [1] "esriGeometryPoint"

plot(dz_cntr)

NHS Health Boards bounadaries - above method does not work for polygons in this situation

hb_pol = st_read("data/SG_NHS_HealthBoards_2019.shp")

## Reading layer `SG_NHS_HealthBoards_2019' from data source `/Users/michal/Documents/kdepop/data/SG_NHS_HealthBoards_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 7564.996 ymin: 530356.4 xmax: 468812.4 ymax: 1218625
## projected CRS:  OSGB 1936 / British National Grid

plot(hb_pol)

Area of Interest

Here I will select NHS Lothian Health Board and clip the Data Zones centroids to its boundaries

Filter

hb_lothian = hb_pol %>% 
  filter(HBName == "Lothian") %>% 
  select(-Shape_Leng, -Shape_Area)

plot(hb_lothian)

Subset

dz_lothian = dz_cntr[hb_lothian,]

plot(dz_lothian)

## Binning

GRID

rec_grid = st_make_grid(dz_lothian, cellsize = 500, square = TRUE)

rec_grid_sf = st_sf(rec_grid)

rec_grid_sf <- rec_grid_sf %>% mutate(id = row_number())

rec_grid_join = st_join(rec_grid_sf, dz_lothian)

rec_grid_sum = rec_grid_join %>%
  group_by(id) %>%
  summarize(
    dz_count = n(),
    pop_sum  = sum(TotPop2011)
  ) %>% 
  filter(!is.na(pop_sum))

## `summarise()` ungrouping output (override with `.groups` argument)

Map

tmap_mode(mode = c("view"))

## tmap mode set to interactive viewing

tm_shape(rec_grid_sum) +
  tm_fill(col = "pop_sum",
          style = "jenks",
          convert2density = TRUE) + 
  tm_borders(alpha=.8, col = "white")

ESRI - Point Density (SPatial Analyst) - focal operation with movable window e.g. 1 point / (9 cells * 100m2 area per cell) = 0.0011 points per square meter (within a 3x3 search window)

https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/point-density.htm

Kernel Density Estimation

Definition

Kernel Density Estimation (KDE) - method of analyzing the first-order of intensity of point pattern

• the kernel
• the bandwidth
• the edge effect

https://mathisonian.github.io/kde/

Library(spatstat)

Convert sf class to ppp class - DZ Centroids

dz_lothian.sp <- as_Spatial(dz_lothian)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in CRS definition

class(dz_lothian.sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

dz_lothian.ppp <- as(dz_lothian.sp, "ppp")

class(dz_lothian.ppp)

## [1] "ppp"

Convert sf class to ppp class - HB Boundaries

hb_lothian.sp <- as_Spatial(hb_lothian)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in CRS definition

class(hb_lothian.sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

hb_lothian.ppp <- as(hb_lothian.sp, "owin")

class(hb_lothian.ppp)

## [1] "owin"

Plot the Point Pattern for Analysis

# assign boundaries to dz ppp

Window(dz_lothian.ppp) <- hb_lothian.ppp

# plot
plot(dz_lothian.ppp, main=NULL, cols=rgb(0,0,0,.2), pch=20)

Create KDE

kde_200 = density(dz_lothian.ppp, 
        sigma=200, # choose bandwith / diameter of the Kernel in the units your map is in
        eps = 100, # pixel resolution
        weights= dz_lothian.ppp$marks, # population 
        edge=TRUE, 
        varcov=NULL,
        at="pixels",
        eaveoneout=TRUE,
        adjust=1, 
        diggle=FALSE, 
        se=FALSE,
        kernel="gaussian", # choose kernel
        scalekernel=is.character(kernel), 
        positive=FALSE,
        verbose=TRUE)

kde_500 = density(dz_lothian.ppp, 
        sigma=500, # choose bandwith / diameter of the Kernel in the units your map is in
        eps = 100, # pixel resolution
        weights= dz_lothian.ppp$marks, # population
        edge=TRUE, 
        varcov=NULL,
        at="pixels",
        eaveoneout=TRUE,
        adjust=1, 
        diggle=FALSE, 
        se=FALSE,
        kernel="gaussian", # choose kernel
        scalekernel=is.character(kernel), 
        positive=FALSE,
        verbose=TRUE)

Plot KDE

plot(kde_200, main = "Sigma 200m")

plot(kde_500, main = "Sigma 500m")

ESRI - Kernel Density (Spatial Analyst)

https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/kernel-density.htm

Extract Clusters

kde_raster = raster(kde_500)

kde_pol = spex::polygonize(kde_raster)

kde_pol2 = kde_pol %>% 
   filter(layer > 0.001) %>% # DEFINE THRESHOLD
   summarize() 

kde_pol2 = kde_pol2 %>% st_set_crs(27700)

plot(kde_pol2)

ESRI - Raster to Polygon (Conversion)

https://pro.arcgis.com/en/pro-app/tool-reference/conversion/raster-to-polygon.htm

Settlements

# settlements
sett_url =  "http://sedsh127.sedsh.gov.uk/arcgis/rest/services/NRS/NRS/MapServer/9"

sett_pol = esri2sf(sett_url, outFields=c("name")) %>% st_transform(27700)

## [1] "Feature Layer"
## [1] "esriGeometryPolygon"

sett_lothian = sett_pol[hb_lothian,]

Map

tmap_mode(mode = c("view"))

## tmap mode set to interactive viewing

tm_shape(kde_pol2) +
  tm_fill(col = "red", alpha = 0.5) +
tm_shape(sett_lothian) +
  tm_fill(col = "blue", alpha = 0.5)

## Warning: The shape sett_lothian is invalid. See sf::st_is_valid

Finally, session info

sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] spex_0.7.1          maptools_1.0-2      rgdal_1.5-18       
##  [4] spatstat_1.64-1     rpart_4.1-15        nlme_3.1-148       
##  [7] spatstat.data_1.5-2 esri2sf_0.1.1       dplyr_1.0.2        
## [10] tmaptools_3.0       tmap_3.0            mapview_2.7.8      
## [13] sf_0.9-6            raster_3.1-5        sp_1.4-4           
## 
## loaded via a namespace (and not attached):
##  [1] magic_1.5-9             httr_1.4.2              jsonlite_1.7.2         
##  [4] viridisLite_0.3.0       splines_4.0.2           stats4_4.0.2           
##  [7] yaml_2.2.1              pillar_1.4.7            lattice_0.20-41        
## [10] glue_1.4.2              digest_0.6.27           RColorBrewer_1.1-2     
## [13] polyclip_1.10-0         leaflet.providers_1.9.0 colorspace_1.4-1       
## [16] htmltools_0.5.0         Matrix_1.2-18           XML_3.99-0.3           
## [19] pkgconfig_2.0.3         stars_0.4-1             purrr_0.3.4            
## [22] scales_1.1.1            webshot_0.5.2           tensor_1.5             
## [25] satellite_1.0.2         spatstat.utils_1.17-0   tibble_3.0.4           
## [28] mgcv_1.8-31             generics_0.1.0          ellipsis_0.3.1         
## [31] reproj_0.4.2            leafsync_0.1.0          magrittr_2.0.1         
## [34] crayon_1.3.4            deldir_0.2-3            evaluate_0.14          
## [37] crsmeta_0.3.0           foreign_0.8-80          lwgeom_0.2-3           
## [40] class_7.3-17            tools_4.0.2             lifecycle_0.2.0        
## [43] stringr_1.4.0.9000      munsell_0.5.0           compiler_4.0.2         
## [46] e1071_1.7-4             quadmesh_0.4.5          rlang_0.4.9            
## [49] classInt_0.4-3          units_0.6-7             grid_4.0.2             
## [52] dichromat_2.0-0         htmlwidgets_1.5.1       goftest_1.2-2          
## [55] crosstalk_1.1.0.1       proj4_1.0-10            leafem_0.1.1           
## [58] base64enc_0.1-3         rmarkdown_2.3           geometry_0.4.5         
## [61] codetools_0.2-16        curl_4.3                abind_1.4-5            
## [64] DBI_1.1.0               R6_2.5.0                knitr_1.29             
## [67] KernSmooth_2.23-17      stringi_1.5.3           parallel_4.0.2         
## [70] Rcpp_1.0.5              vctrs_0.3.5             png_0.1-7              
## [73] leaflet_2.0.3           tidyselect_1.1.0        xfun_0.16

 

---

A work by Michal Michalski