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) |>
addTiles() |>
addCircleMarkers(radius = 1)Estimating densities from a point pattern
statistics
point pattern
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.
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_lambertOnce 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] 1expand.grid(lon = dens$x, lat = dens$y) |>
mutate(
density = as.vector(dens$z) * length(points_lambert),
id = seq_along(.data$density)
) -> densIn 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_wgs84leaflet 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) |>
addTiles() |>
addPolygons(color = ~pal(density), stroke = FALSE, fillOpacity = 0.5) |>
addLegend(pal = pal, values = ~density)
mapSession 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 [1] CRAN (R 4.3.0)
callr 3.7.3 2022-11-02 [1] CRAN (R 4.3.0)
class 7.3-22 2023-05-03 [4] CRAN (R 4.3.1)
classInt 0.4-9 2023-02-28 [1] CRAN (R 4.3.1)
cli 3.6.1 2023-03-23 [1] CRAN (R 4.3.0)
codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.0)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0)
crosstalk 1.2.0 2021-11-04 [1] CRAN (R 4.3.0)
DBI 1.1.3 2022-06-18 [1] CRAN (R 4.3.0)
digest 0.6.32 2023-06-26 [1] CRAN (R 4.3.1)
dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.3.0)
e1071 1.7-13 2023-02-01 [1] CRAN (R 4.3.1)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.0)
evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0)
fansi 1.0.4 2023-01-22 [1] CRAN (R 4.3.0)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0)
fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.3.0)
htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.3.0)
jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.1)
KernSmooth 2.23-21 2023-05-03 [1] CRAN (R 4.3.0)
knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0)
lattice 0.21-8 2023-04-05 [4] CRAN (R 4.3.0)
leafem 0.2.0 2022-04-16 [1] CRAN (R 4.3.1)
leaflet * 2.1.2 2023-03-10 [1] CRAN (R 4.3.0)
lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.3.0)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
mapview 2.11.0 2022-04-16 [1] CRAN (R 4.3.1)
MASS * 7.3-60 2023-05-04 [4] CRAN (R 4.3.1)
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0)
pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
processx 3.8.2 2023-06-30 [1] CRAN (R 4.3.1)
proxy 0.4-27 2022-06-09 [1] CRAN (R 4.3.1)
ps 1.7.5 2023-04-18 [1] CRAN (R 4.3.0)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
raster 3.6-23 2023-07-04 [1] CRAN (R 4.3.1)
Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.3.0)
rlang 1.1.1 2023-04-28 [1] CRAN (R 4.3.0)
rmarkdown 2.23 2023-07-01 [1] CRAN (R 4.3.1)
rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.3.0)
satellite 1.0.4 2021-10-12 [1] CRAN (R 4.3.1)
scales 1.2.1 2022-08-20 [1] CRAN (R 4.3.0)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
sf * 1.0-13 2023-05-24 [1] CRAN (R 4.3.0)
sp 2.0-0 2023-06-22 [1] CRAN (R 4.3.1)
terra 1.7-39 2023-06-23 [3] CRAN (R 4.3.1)
tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0)
units 0.8-2 2023-04-27 [1] CRAN (R 4.3.0)
utf8 1.2.3 2023-01-31 [1] CRAN (R 4.3.0)
vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.3.0)
webshot 0.5.5 2023-06-26 [1] CRAN (R 4.3.1)
xfun 0.39 2023-04-20 [1] CRAN (R 4.3.0)
yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0)
[1] /home/thierry/R/x86_64-pc-linux-gnu-library/4.3
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
──────────────────────────────────────────────────────────────────────────────