# Estimating densities from a point pattern

statistics
point pattern
Author

Thierry Onkelinx

Published

June 28, 2017

In this example we focus on a set of 10450 coordinates in a small area. The goal is to estimate the local density of points, expressed as the number of point per unit area. The raw coordinates are given in WGS84 (EPSG:4326), which is a geodetic coordinate system. That is not suited for calculating distances, so we need to re-project the points into a local projected coordinate system. In this case we use Lambert72 (EPSG:3170). Next we calculate the density. To visualise the density, we have to transform the results back in to WGS84.

The data used in this example is real data by centred to a different location for privacy reasons. The dataset is available on GitHub.

First we must read the data into R. Plotting the raw data helps to check errors in the data (Figure 1). The easiest way to plot it, is to convert the dataframe into a simple features object.

``````library(leaflet)
library(sf)
crs_wgs84 <- st_crs(4326)
read.delim("data/points.txt", sep = " ") |>
st_as_sf(coords = c("lon", "lat"), crs = crs_wgs84) -> points
leaflet(points) |>

The original coordinates are in the WGS84 projection. We need to project them onto a plane by using an appropriate local coordinate system, in this case Lambert72 (EPSG code 31370).

``````crs_lambert <- st_crs(31370)
st_transform(points, crs = crs_lambert) |>
st_coordinates() |>
as.data.frame() -> points_lambert``````

Once we have the points into a projected coordinate system, we can calculate the densities. We start by defining a grid. `cellsize` is the dimension of the square grid cell in the units of the projected coordinate system. Meters in case of Lambert72. The boundaries of the grid are defined using `pretty()`, which turns a vector of numbers into a “pretty” vector with rounded numbers. The combination of the boundaries and the cell size determine the number of grid cells `n` in each dimension. `diff()` calculates the difference between to adjacent numbers of a vector. The density is calculated with `MASS::kde2d()` based on the vectors with the longitude and latitude, the number of grid cells in each dimension and the boundaries of the grid. This returns the grid as a list with elements `x` (a vector of longitude coordinates of the centroids), `y` (a vector of latitude coordinates of the centroids) and `z` (a matrix with densities). The values in `z` are densities for the ‘average’ point per unit area. When we multiply the value `z` with the area of the grid cell and sum all of them we get 1. So if we multiple `z` with the number of points we get the density of the points per unit area.

We use `dplyr::mutate()` to convert it into a `data.frame`. The last two steps convert the centroids into a set of coordinates for square polygons.

``````library(MASS)
library(dplyr)
xlim <- range(pretty(points_lambert\$X)) + c(-100, 100)
ylim <- range(pretty(points_lambert\$Y)) + c(-100, 100)
n <- c(diff(xlim), diff(ylim)) / params\$cellsize + 1
dens <- kde2d(
x = points_lambert\$X, y = points_lambert\$Y, n = n, lims = c(xlim, ylim)
)
dx <- diff(dens\$x[1:2])
dy <- diff(dens\$y[1:2])
sum(dens\$z * dx * dy)``````
`` 1``
``````expand.grid(lon = dens\$x, lat = dens\$y) |>
mutate(
density = as.vector(dens\$z) * length(points_lambert),
id = seq_along(.data\$density)
) -> dens``````

In order to visualise the result, we have to do a few step. First convert the dataframe into a simple object with points. Then convert the points into a grid of polygons. Next we attach the estimated density as counts per grid cell. Finally we re-project the coordinates back to WGS84. Then we can display the grid with a web based background image (Figure 2).

``````dens |>
st_as_sf(coords = c("lon", "lat"), crs = crs_lambert) |>
st_make_grid(
cellsize = c(dx, dy),
offset = c(min(dens\$lon) - dx / 2, min(dens\$lat) - dy / 2)
) |>
st_as_sf() |>
mutate(
density = dens\$density * params\$cellsize ^ 2
) |>
st_transform(crs = crs_wgs84) -> dens_wgs84``````

`leaflet` requires a predefined function with a colour pallet. We use `leaflet::colorNumeric()` to get a continuous pallet. Setting `stroke = FALSE` removes the borders of the polygon. `fillOpacity` sets the transparency of the polygons.

``````pal <- colorNumeric(
palette = rev(rainbow(100, start = 0, end = .7)),
domain = c(0, dens_wgs84\$density)
)
map <- leaflet(dens_wgs84) |>
addPolygons(color = ~pal(density), stroke = FALSE, fillOpacity = 0.5) |>
addLegend(pal = pal, values = ~density)
map``````

## Session info

These R packages were used to create this post.

``sessioninfo::session_info()``
``````─ Session info ───────────────────────────────────────────────────────────────
setting  value
version  R version 4.3.1 (2023-06-16)
os       Ubuntu 22.04.3 LTS
system   x86_64, linux-gnu
ui       X11
language nl_BE:nl
collate  nl_BE.UTF-8
ctype    nl_BE.UTF-8
tz       Europe/Brussels
date     2023-09-02
pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
package     * version date (UTC) lib source
base64enc     0.1-3   2015-07-28  CRAN (R 4.3.0)
callr         3.7.3   2022-11-02  CRAN (R 4.3.0)
class         7.3-22  2023-05-03  CRAN (R 4.3.1)
classInt      0.4-9   2023-02-28  CRAN (R 4.3.1)
cli           3.6.1   2023-03-23  CRAN (R 4.3.0)
codetools     0.2-19  2023-02-01  CRAN (R 4.3.0)
colorspace    2.1-0   2023-01-23  CRAN (R 4.3.0)
crosstalk     1.2.0   2021-11-04  CRAN (R 4.3.0)
DBI           1.1.3   2022-06-18  CRAN (R 4.3.0)
digest        0.6.32  2023-06-26  CRAN (R 4.3.1)
dplyr       * 1.1.2   2023-04-20  CRAN (R 4.3.0)
e1071         1.7-13  2023-02-01  CRAN (R 4.3.1)
ellipsis      0.3.2   2021-04-29  CRAN (R 4.3.0)
evaluate      0.21    2023-05-05  CRAN (R 4.3.0)
fansi         1.0.4   2023-01-22  CRAN (R 4.3.0)
farver        2.1.1   2022-07-06  CRAN (R 4.3.0)
fastmap       1.1.1   2023-02-24  CRAN (R 4.3.0)
generics      0.1.3   2022-07-05  CRAN (R 4.3.0)
glue          1.6.2   2022-02-24  CRAN (R 4.3.0)
htmltools     0.5.5   2023-03-23  CRAN (R 4.3.0)
htmlwidgets   1.6.2   2023-03-17  CRAN (R 4.3.0)
jsonlite      1.8.7   2023-06-29  CRAN (R 4.3.1)
KernSmooth    2.23-21 2023-05-03  CRAN (R 4.3.0)
knitr         1.43    2023-05-25  CRAN (R 4.3.0)
lattice       0.21-8  2023-04-05  CRAN (R 4.3.0)
leafem        0.2.0   2022-04-16  CRAN (R 4.3.1)
leaflet     * 2.1.2   2023-03-10  CRAN (R 4.3.0)
lifecycle     1.0.3   2022-10-07  CRAN (R 4.3.0)
magrittr      2.0.3   2022-03-30  CRAN (R 4.3.0)
mapview       2.11.0  2022-04-16  CRAN (R 4.3.1)
MASS        * 7.3-60  2023-05-04  CRAN (R 4.3.1)
munsell       0.5.0   2018-06-12  CRAN (R 4.3.0)
pillar        1.9.0   2023-03-22  CRAN (R 4.3.0)
pkgconfig     2.0.3   2019-09-22  CRAN (R 4.3.0)
png           0.1-8   2022-11-29  CRAN (R 4.3.0)
processx      3.8.2   2023-06-30  CRAN (R 4.3.1)
proxy         0.4-27  2022-06-09  CRAN (R 4.3.1)
ps            1.7.5   2023-04-18  CRAN (R 4.3.0)
R6            2.5.1   2021-08-19  CRAN (R 4.3.0)
raster        3.6-23  2023-07-04  CRAN (R 4.3.1)
Rcpp          1.0.10  2023-01-22  CRAN (R 4.3.0)
rlang         1.1.1   2023-04-28  CRAN (R 4.3.0)
rmarkdown     2.23    2023-07-01  CRAN (R 4.3.1)
rstudioapi    0.14    2022-08-22  CRAN (R 4.3.0)
satellite     1.0.4   2021-10-12  CRAN (R 4.3.1)
scales        1.2.1   2022-08-20  CRAN (R 4.3.0)
sessioninfo   1.2.2   2021-12-06  CRAN (R 4.3.0)
sf          * 1.0-13  2023-05-24  CRAN (R 4.3.0)
sp            2.0-0   2023-06-22  CRAN (R 4.3.1)
terra         1.7-39  2023-06-23  CRAN (R 4.3.1)
tibble        3.2.1   2023-03-20  CRAN (R 4.3.0)
tidyselect    1.2.0   2022-10-10  CRAN (R 4.3.0)
units         0.8-2   2023-04-27  CRAN (R 4.3.0)
utf8          1.2.3   2023-01-31  CRAN (R 4.3.0)
vctrs         0.6.3   2023-06-14  CRAN (R 4.3.0)
webshot       0.5.5   2023-06-26  CRAN (R 4.3.1)
xfun          0.39    2023-04-20  CRAN (R 4.3.0)
yaml          2.3.7   2023-01-23  CRAN (R 4.3.0)

 /home/thierry/R/x86_64-pc-linux-gnu-library/4.3
 /usr/local/lib/R/site-library
 /usr/lib/R/site-library
 /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────``````