In JeremyGelb/spNetwork: Spatial Analysis on Network

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
backup_option <- options()
base_wd <- getwd()

Introduction

The goal of this vignette is to present the differences between the three Network Kernel Density Estimation methods implemented in the package spNetwork. They are presented with 3D rendering to get a better understanding of their differences.

For a more detailed description, the reader should read the provided references.

Let us start with the definition of a simple situation with one event and a small network. The spatial coordinates are given in metres.

options(rgl.useNULL=TRUE)
library(spNetwork)
library(sf)
library(rgl)
library(tmap)
library(dbscan)

# we define a set of lines
wkt_lines <- c(
    "LINESTRING (0 0, 1 0)", "LINESTRING (1 0, 2 0)", "LINESTRING (2 0, 3 0)",
    "LINESTRING (0 1, 1 1)", "LINESTRING (1 1, 2 1)", "LINESTRING (2 1, 3 1)",
    "LINESTRING (0 2, 1 2)", "LINESTRING (1 2, 2 2)", "LINESTRING (2 2, 3 2)",
    "LINESTRING (0 3, 1 3)", "LINESTRING (1 3, 2 3)", "LINESTRING (2 3, 3 3)",
    "LINESTRING (0 0, 0 1)", "LINESTRING (0 1, 0 2)", "LINESTRING (0 2, 0 3)",
    "LINESTRING (1 0, 1 1)", "LINESTRING (1 1, 1 2)", "LINESTRING (1 2, 1 3)",
    "LINESTRING (2 0, 2 1)", "LINESTRING (2 1, 2 2)", "LINESTRING (2 2, 2 3)",
    "LINESTRING (3 0, 3 1)", "LINESTRING (3 1, 3 2)", "LINESTRING (3 2, 3 3)"
    )

# and create a spatial lines dataframe
linesdf <- data.frame(wkt = wkt_lines,
                      id = paste("l",1:length(wkt_lines),sep=""))

all_lines <- st_as_sf(linesdf, wkt = "wkt")

# we define now one event
event <- data.frame(x=c(1),
                    y=c(1.5))
event <- st_as_sf(event, coords = c("x","y"))

# and map the situation
tm_shape(all_lines) + 
  tm_lines("black") + 
  tm_shape(event) + 
  tm_dots("red", size = 0.2)

d3_plot_situation <- function(lines, events, pt_samples, densities, scales){

  open3d(scale=scales)

  densities <- round(densities,7)

  ## finding for each event its closest samples
  XY1 <- st_coordinates(pt_samples)
  XY2 <- st_coordinates(events)
  idx <- dbscan::kNN(x = XY1, query = XY2, k = 1)$id
  #idx <- knnx.index(XY1, query = XY2, k = 1)

  events$dens <- densities[idx]
  events$x <- XY2[,1]
  events$y <- XY2[,2]
  eidx <- do.call(c,lapply(1:nrow(events), function(i){c(i,i)}))
  vert_lines <- st_drop_geometry(events[eidx,])
  vert_lines$dens <- ifelse(1:nrow(vert_lines)%%2 == 0, vert_lines$dens,0)


  ## plotting the situation
  #line_coords <- do.call(rbind,unlist(coordinates(lines), recursive = FALSE))
  line_coords <- st_coordinates(lines)
  sample_coords <- st_coordinates(pt_samples)

  segments3d(x=line_coords[,1],
             y=line_coords[,2],
             z=rep(0, nrow(line_coords)))

  segments3d(x=vert_lines$x,
             y=vert_lines$y,
             z=vert_lines$dens,
             col = "red")

  coords_events <- st_coordinates(events)

  points3d(
    x = coords_events[,1],
    y = coords_events[,2],
    z = rep(0,nrow(event)),
    col = "red",
    size = 5
  )

  points3d(
    x = sample_coords[,1],
    y = sample_coords[,2],
    z = densities,
    col = "blue"
  )

  axes3d()
  title3d(xlab="X",ylab="Y",zlab="Z")
}

knitr::knit_hooks$set(webgl = hook_webgl)

The simple method

The first method was proposed by @xie2008kernel. The differences with a classical KDE are:

• the events are snapped on a network
• the distances between sampling points and events are network distances (i.e. not Euclidean distance)

To visualize the result, we calculate here the density values and plot them in 3D. We sample the kernel values at each one centimetres on the lines.

samples_pts <- lines_points_along(all_lines,0.01)

simple_kernel <- nkde(all_lines, event, w = 1,
                  samples = samples_pts, kernel_name = "quartic", 
                  bw = 2, method = "simple", div = "bw", 
                  digits = 3, tol = 0.001, grid_shape = c(1,1),
                  check = FALSE,
                  verbose = FALSE)

Let us see this in 3D !

d3_plot_situation(all_lines, event, samples_pts, simple_kernel, scales = c(1,1,3))

(The function used to display the 3d plot is available in the source code of the vignette)

As one can see, NKDE is only determined by the distance between the sampling points and the event. This is a problem, because at intersections, the event's mass is multiplied by the number of edges at that intersection. In more mathematical words, this is not a true kernel because it does not integrate to 1 on its domain.

This method remains useful for two reasons:

• With big datasets, it might be useful to use this simple method to do a primary investigation (quick calculation).
• In a purely geographical view, this method is intuitive. In the case of crime analysis for example, one could argue that the strength of the event should not be affected by intersections on the network. In that case, the kernel function can be seen as a distance decaying function.

What to retain about the simple method

| Pros | Cons | |:-------------------:|:--------------------------:| |quick to calculate | biased (not a true kernel) | |intuitive | overestimate the densities | |continuous | |

The discontinuous method

To overcome the limits of the previous method, @okabe2012spatial have proposed the discontinuous NKDE, which has been extended by @sugihara2010simple for cases where the network has cycles shorter than the bandwidth.

discontinuous_kernel <- nkde(all_lines, event, w = 1,
                  samples = samples_pts, kernel_name = "quartic",
                  bw = 2, method = "discontinuous", div = "bw",
                  digits = 3, tol = 0.001, grid_shape = c(1,1),
                  check = FALSE,
                  verbose = FALSE)

d3_plot_situation(all_lines, event, samples_pts, discontinuous_kernel, scales = c(1,1,3))

As one can see, the values of the NKDE are split at intersections to avoid the multiplication of the mass observed in the simple version. However, this NKDE is discontinuous, which is counter-intuitive, and leads to sharp differences between density values in the network, and could be problematic in networks with numerous and close intersections.

What to retain about the discontinuous method

| Pros | Cons | |:-------------------:|:--------------------------:| |quick to calculate | counter-intuitive | |unbiased | discontinuous | | | contrasted results |

The continuous method

Finally, the continuous method takes the best of the two worlds: it adjusts the values of the NKDE at intersections to ensure that it integrates to one on its domain, and applies a backward correction to force the density values to be continuous. This process is accomplished by a recursive function, described in the book Spatial Analysis Along Networks [@okabe2012spatial]. This function is very time consuming, so it might be necessary to stop it when the recursion is too deep. Considering that the kernel density is divided at each intersection, stopping the function at deep level 15 produces results that are almost equal to the true values.

continuous_kernel <- nkde(all_lines, event, w = 1,
                  samples = samples_pts, kernel_name = "quartic",
                  bw = 2, method = "continuous", div = "bw",
                  digits = 3, tol = 0.001, grid_shape = c(1,1),
                  check = FALSE,
                  verbose = FALSE)

d3_plot_situation(all_lines, event, samples_pts, continuous_kernel, scales = c(1,1,3))

As one can see, the values of the NKDE are continuous, and the density values close to the events have been adjusted. This method produces smoother results than the discontinuous method.

what to retain about the discontinuous method

| Pros | Cons | |:-------------------:|:--------------------------:| | unbiased | long calculation time | | continuous | smoother values |

If we add a second event, their mass will add up when they share some range.

# we define now two events
events <- data.frame(x=c(1,2),
                    y=c(1.5,1.5))
events <- st_as_sf(events, coords = c("x","y"))

continuous_kernel <- nkde(all_lines, events, w = 1,
                  samples = samples_pts, kernel_name = "quartic",
                  bw = 1.3, method = "continuous", div = "bw",
                  digits = 3, tol = 0.001, grid_shape = c(1,1),
                  check = FALSE,
                  verbose = FALSE)

d3_plot_situation(all_lines, events, samples_pts, continuous_kernel, scales = c(1,1,3))

The larger the bandwidth, the more common links the events share in their range.

continuous_kernel <- nkde(all_lines, events, w = 1,
                  samples = samples_pts, kernel_name = "quartic",
                  bw = 1.6, method = "continuous", div = "bw",
                  digits = 3, tol = 0.001, grid_shape = c(1,1),
                  check = FALSE,
                  verbose = FALSE)

d3_plot_situation(all_lines, events, samples_pts, continuous_kernel, scales = c(1,1,3))

# reset all the user parameters
options(backup_option)
setwd(base_wd)
oldpar <- par(mfrow = c(1,2))
par(oldpar)

References

JeremyGelb/spNetwork documentation built on May 24, 2024, 7:23 p.m.