---
title: "Creating Neighbours using sf objects"
author: "Roger Bivand"
output:
  html_document:
    toc: true
    toc_float:
      collapsed: false
      smooth_scroll: false
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Creating Neighbours using sf objects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Creating Neighbours using sf objects

## Introduction

This vignette tracks the legacy nb vignette, which was based on part of the first (2008) edition of ASDAR. It adds hints to the code in the nb vignette to use the sf vector representation instead of the sp vector representation to create neighbour objects. . 

## Summary

This is a summary of the results below:

- In general, if you need to reproduce results from using `"Spatial"` objects in **spdep**, coerce sf objects to sp objects before constructing neighbour objects (particularly if polygon centroids are used for point representation).

- However, for new work, you should use `"sf"` objects read in using **sf**.

- From **spdep** 1.1-7, a number of steps have been taken to choose more efficient approaches especially for larger data sets, using functions in **sf** and the approximate nearest neighbour (ANN) implementation in **dbscan** rather than **RANN**.

## Data set

We'll use the whole NY 8 county set of boundaries, as they challenge the implementations more than just the Syracuse subset. The description of the input geometries from ADSAR is: New York leukemia: used and documented extensively in Waller and Gotway (2004) and with data formerly made available in Chap. 9 of `http://web1.sph.emory.edu/users/lwaller/ch9index.htm`; the data import process is described in the help file of NY_data in spdep; geometries downloaded from the CIESIN server at [https://sedac.ciesin.columbia.edu/data/set/acrp-boundary-1992/data-download](https://sedac.ciesin.columbia.edu/data/set/acrp-boundary-1992/data-download), file /pub/census/usa/tiger/ny/bna_st/t8_36.zip, and extensively edited; a zip archive NY_data.zip of shapefiles and a GAL format neighbours list is on the book website. Further, the zipfile is now at: [a new location requiring login]
(https://sedac.ciesin.columbia.edu/ftpsite/pub/census/usa/tiger/ny/bna_st/t8_36.zip). The object listw_NY is directly imported from nyadjwts.dbf on the Waller & Gotway (2004) chapter 9 website.

The version of the New York 8 counties geometries used in ASDAR and included as a shapefile in spdep was converted from the original BNA file using an external utility program to convert to MapInfo format and converted on from there using GDAL 1.4.1 (the OGR BNA driver was not then available; it entered OGR at 1.5.0, release at the end of 2007), and contains invalid geometries. What was found in mid-2007 was that included villages were in/excluded by in-out umbilical cords to the boundary of the enclosing tract, when the underlying BNA file was first converted to MapInfo (holes could not exist then).

Here we will use a GPKG file created as follows (rgdal could also be used with the same output; GDAL here is built with GEOS, so the BNA vector driver will use geometry tests: The BNA driver supports reading of polygons with holes or lakes. It determines what is a hole or a lake only from geometrical analysis (inclusion, non-intersection tests) and ignores completely the notion of polygon winding (whether the polygon edges are described clockwise or counter-clockwise). GDAL must be built with GEOS enabled to make geometry test work.):

```{r, eval=FALSE}
library(sf)
sf_bna <- st_read("t8_36.bna", stringsAsFactors=FALSE)
table(st_is_valid(sf_bna))
sf_bna$AREAKEY <- gsub("\\.", "", sf_bna$Primary.ID)
data(NY_data, package="spData")
key <- as.character(nydata$AREAKEY)
sf_bna1 <- sf_bna[match(key, sf_bna$AREAKEY), c("AREAKEY")]
sf_bna2 <- merge(sf_bna1, nydata, by="AREAKEY")
sf_bna2_utm18 <- st_transform(sf_bna2, "+proj=utm +zone=18 +datum=NAD27")
st_write(sf_bna2_utm18, "NY8_bna_utm18.gpkg")
```


## nb and listw objects (copied from the nb_igraph vignette)

Since the **spdep** package was created, *spatial weights* objects have been constructed as lists with three components and a few attributes, in old-style class `listw` objects. The first component of a `listw` object is an `nb` object, a list of `n` integer vectors, with at least a character vector `region.id` attribute with `n` unique values (like the `row.names` of a `data.frame` object); `n` is the number of spatial entities. Component `i` of this list contains the integer identifiers of the neighbours of `i` as a sorted vector with no duplication and values in `1:n`; if `i` has no neighbours, the component is a vector of length `1` with value `0L`. The `nb` object may contain an attribute indicating whether it is symmetric or not, that is whether `i` is a neighbour of `j` implies that `j` is a neighbour of `i`. Some neighbour definitions are symmetric by construction, such as contiguities or distance thresholds, others are asymmetric, such as `k`-nearest neighbours. The `nb` object redundantly stores both `i`-`j` and `j`-`i` links.

The second component of a `listw` object is a list of `n` numeric vectors, each of the same length as the corresponding non-zero vectors in the `nb`object. These give the values of the spatial weights for each `i`-`j` neighbour pair. It is often the case that while the neighbours are symmetric by construction, the weights are not, as for example when weights are *row-standardised* by dividing each row of input weights by the count of neighbours or cardinality of the neighbour set of `i`. In the `nb2listw`function, it is also possible to pass through general weights, such as inverse distances, shares of boundary lengths and so on.

The third component of a `listw` object records the `style` of the weights as a character code, with `"B"` for binary weights taking values zero or one (only one is recorded), `"W"` for row-standardised weights, and so on. In order to subset `listw` objects, knowledge of the `style` may be necessary.

## Comparison of sp and sf approaches

First some housekeeping and setup to permit this vignette to be built when packages are missing or out-of-date:

```{r, echo=FALSE}
rv <- R.Version()
dothis <- FALSE
if (rv$major > "3" || (rv$major == "3" && !(rv$minor >= "3.0"))) dothis=TRUE
```

```{r, echo=dothis, eval=dothis}
if (!suppressPackageStartupMessages(require(sf, quietly=TRUE))) {
  message("install the sf package")
  dothis <- FALSE
}
if (dothis) sf_extSoftVersion()
```

Let us read the GPKG file with valid geometries in to 'sf' and 'sp' objects:

```{r, echo=dothis, eval=dothis}
NY8_sf <- st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData"), quiet=TRUE)
table(st_is_valid(NY8_sf))
```

### Contiguity neighbours for polygon support

Here we first generate a queen contiguity nb object using the legacy spdep approach. This first either uses a pre-computed list of vectors of probable neighbours or finds intersecting bounding boxes internally. Then the points on the boundaries of each set of polygons making up an observation are checked for a distance less than snap to any of the points of the set of polygons making up an observation included in the set of candidate neighbours. Because contiguity is symmetric, only i to j contiguities are tested. A queen contiguity is found as soon as one point matches, a rook contiguity as soon as two points match:

```{r, echo=dothis, eval=dothis}
suppressPackageStartupMessages(library(spdep))
reps <- 10
eps <- sqrt(.Machine$double.eps)
system.time(for(i in 1:reps) NY8_sf_1_nb <- poly2nb(NY8_sf, queen=TRUE, snap=eps))/reps
```
Using STR trees to check the intersection of envelopes (bounding boxes) is much faster than the internal method in poly2nb for large n (the default threshold is `small_n=500`). From **spdep** 1.1-7, use is made of GEOS through **sf** to find candidate neighbours. Because contiguity is symmetric by definition, `foundInBox=` only requires intersections for higher indices, leading to a slight overhead to remove duplicates, as `st_intersects()` reports both `i j` ans `j i` relationships.

```{r, echo=dothis, eval=dothis}
system.time(for(i in 1:reps) NY8_sf_1_nb_STR <- poly2nb(NY8_sf, queen=TRUE, snap=eps, small_n=250))/reps
```
In this case, n is small, so the extra work in GEOS has a higher overhead than the simple, legacy approach.

```{r, echo=dothis, eval=dothis}
NY8_sf_1_nb
```

```{r, echo=dothis, eval=dothis}
all.equal(NY8_sf_1_nb, NY8_sf_1_nb_STR, check.attributes=FALSE)
```

spdep::poly2nb uses two heuristics, first to find candidate neighbours from intersecting bounding boxes, and second to use the symmetry of the relationship to halve the number of remaining tests. This means that performance is linear in n, but with overhead for identifying candidates, and back-filling symmetric neighbours. Further, spdep::poly2nb stops searching for queen contiguity as soon as the first neighbour point is found within snap distance (if not identical, which is tested first); second neighbour point for rook contiguities. spdep::poly2nb was heavily optimised when written, as processor speed was a major constraint at that time. 

The addition of STR tree queries to identify candidates permits the construction of contiguous neighbour objects for quite large objects, for example the ZCTA shapefile with 33144 features. sf::st_read imports the data in less than a second and the polygon geometries are valid. Finding the candidate neighbours with `st_intersects( st_as_sfc( lapply( pl, function(x) {st_as_sfc(st_bbox(x))[[1]]} ) ) )` where `pl` is an `"sfc"` object takes about 5.6 s, `lapply` to remove duplicate `i j / j i` relationships, and a total run time of about 17 s. Setting the `small_n=` threshold above 33144 gives the legacy approach, about 110 s. The legacy `snap=` argument works in the same way with **sf** or **sp** objects.


### Contiguity neighbours from invalid polygons

Next, we explore a further possible source of differences in neighbour object reproduction, using the original version of the tract boundaries used in ASDAR, but with some invalid geometries as mentioned earlier:

```{r, echo=dothis, eval=dothis}
NY8_sf_old <- st_read(system.file("shapes/NY8_utm18.shp", package="spData"), quiet=TRUE)
table(st_is_valid(NY8_sf_old))
```

We can see that there are a number of differences between the neighbour sets derived from the fully valid geometries and the older partly invalid set:

```{r, echo=dothis, eval=dothis}
try(NY8_sf_old_1_nb <- poly2nb(NY8_sf_old), silent = TRUE)
all.equal(NY8_sf_old_1_nb, NY8_sf_1_nb, check.attributes=FALSE)
```

spdep::poly2nb does not object to using invalid geometries, as it only uses the boundary points defining the polygons (as do STR tree construction and query functions, because valid bounding boxes can be constructed from invalid polygons).

Using the standard "trick", we can buffer by 0 to try to make the geometries valid:

```{r, echo=dothis, eval=dothis}
NY8_sf_old_buf <- st_buffer(NY8_sf_old, dist=0)
table(st_is_valid(NY8_sf_old_buf))
```
Hoverver, in doing this, we change the geometries, so the new sets of neighbours still differ from those made with the valid geometries in the same ways as before imposing validity:

```{r, echo=dothis, eval=dothis}
try(NY8_sf_old_1_nb_buf <- poly2nb(NY8_sf_old_buf), silent = TRUE)
all.equal(NY8_sf_old_1_nb_buf, NY8_sf_1_nb, check.attributes=FALSE)
```

The neighbour sets are the same for the old boundaries with or without imposing validity:

```{r, echo=dothis, eval=dothis}
all.equal(NY8_sf_old_1_nb_buf, NY8_sf_old_1_nb, check.attributes=FALSE)
```

## Point-based neighbours

### Finding points for polygon objects

`knearneigh()` and `dnearneigh()` require point support, so decisions must be taken about how to place the point in the areal object. We can use `st_centroid()` to get the centroids, using the `of_largest_polygon=TRUE` argument to make sure that the centroid is that of the largest polygon id the observation is made up of more than one external ring:

```{r, echo=dothis, eval=dothis}
NY8_ct_sf <- st_centroid(st_geometry(NY8_sf), of_largest_polygon=TRUE)
```

or `st_point_on_surface()` which guarantees that the point will fall on the surface of a member polygon:

```{r, echo=dothis, eval=dothis}
NY8_pos_sf <- st_point_on_surface(st_geometry(NY8_sf))
```

or indeed taking the centre of the largest inscribed circle (the function returns a radius line segment, so we choose the central point, not the point on the circle):

```{r, echo=dothis, eval=dothis}
if (unname(sf_extSoftVersion()["GEOS"] >= "3.9.0")) 
    NY8_cic_sf <- st_cast(st_inscribed_circle(st_geometry(NY8_sf), nQuadSegs=0), "POINT")[(1:(2*nrow(NY8_sf)) %% 2) != 0]
```

We need to check whether coordinates are planar or not:

```{r, echo=dothis, eval=dothis}
st_is_longlat(NY8_ct_sf)
```

### Graph-based neighbours


From this, we can check the graph-based neighbours (planar coordinates only):

```{r, echo=dothis, eval=dothis}
suppressPackageStartupMessages(require(deldir))
NY84_nb <- tri2nb(NY8_ct_sf)
if (require(dbscan, quietly=TRUE)) {
  NY85_nb <- graph2nb(soi.graph(NY84_nb, NY8_ct_sf))
} else NY85_nb <- NULL
NY86_nb <- graph2nb(gabrielneigh(NY8_ct_sf))
NY87_nb <- graph2nb(relativeneigh(NY8_ct_sf))
```

### K-nearest neighbours

K-nearest neighbours use the coordinate matrices, and can handle Great Circle distances, but this is not demonstrated here, as the data set used is planar, in which case `dbscan::kNN()` in 2D or 3D building a kd-tree is used:

```{r, echo=dothis, eval=dothis}
system.time(for (i in 1:reps) NY88_nb_sf <- knn2nb(knearneigh(NY8_ct_sf, k=1)))/reps
```
Legacy code may be used omitting the kd-tree:

```{r, echo=dothis, eval=dothis}
system.time(for (i in 1:reps) NY89_nb_sf <- knn2nb(knearneigh(NY8_ct_sf, k=1, use_kd_tree=FALSE)))/reps
```

### Distance neighbours

Distance neighbours need a threshold - `nbdists` shows the maximum distance to first nearest neighbour:

```{r, echo=dothis, eval=dothis}
dsts <- unlist(nbdists(NY88_nb_sf, NY8_ct_sf))
summary(dsts)
max_1nn <- max(dsts)
```

`dnearneigh` can also handle Great Circle distances, but this is not demonstrated here, as the data set used is planar:

```{r, echo=dothis, eval=dothis}
system.time(for (i in 1:reps) NY810_nb <- dnearneigh(NY8_ct_sf, d1=0, d2=0.75*max_1nn))/reps
```
By default, the function uses `dbscan::frNN()` to build a kd-tree in 2D or 3D which is then used to find distance neighbours. For small n, the argument `use_kd_tree=FALSE` may speed up computation a little by reverting to legacy code not building a kd-tree first, but in general the differences are so small that the user will not notice:

```{r, echo=dothis, eval=dothis}
system.time(for (i in 1:reps) NY811_nb <- dnearneigh(NY8_ct_sf, d1=0, d2=0.75*max_1nn, use_kd_tree=FALSE))/reps
```
