This vignette tracks the legacy nb vignette, which was based on partof the first (2008) edition of ASDAR. It adds hints to the code in thenb vignette to use the sf vector representation instead of the sp vectorrepresentation to create neighbour objects. .
This is a summary of the results below:
In general, if you need to reproduce results from using"Spatial" objects inspdep, coerce sfobjects to sp objects before constructing neighbour objects(particularly if polygon centroids are used for pointrepresentation).
However, for new work, you should use"sf" objectsread in usingsf.
Fromspdep 1.1-7, a number of steps have beentaken to choose more efficient approaches especially for larger datasets, using functions insf and the approximate nearestneighbour (ANN) implementation indbscan rather thanRANN.
We’ll use the whole NY 8 county set of boundaries, as they challengethe implementations more than just the Syracuse subset. The descriptionof the input geometries from ADSAR is: New York leukemia: used anddocumented extensively in Waller and Gotway (2004) and with dataformerly made available in Chap. 9 ofhttp://web1.sph.emory.edu/users/lwaller/ch9index.htm; thedata import process is described in the help file of NY_data inspData; geometries downloaded from the CIESIN serverformerly atsedac.ciesin.columbia.edu/data/set/acrp-boundary-1992/data-download, nowpossiblyearthdata.nasa.gov/data/catalog/sedac-ciesin-sedac-acrp-1992-bf-1.00,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 listis on the book website. Further, the zipfile is now at a new locationrequiring login. The object listw_NY is directly imported fromnyadjwts.dbf on the Waller & Gotway (2004) chapter 9 website.
The version of the New York 8 counties geometries used in ASDAR andincluded as a shapefile in spdep was converted from the original BNAfile using an external utility program to convert to MapInfo format andconverted on from there using GDAL 1.4.1 (the OGR BNA driver was notthen 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 thatincluded villages were in/excluded by in-out umbilical cords to theboundary of the enclosing tract, when the underlying BNA file was firstconverted to MapInfo (holes could not exist then).
Here we will use a GPKG file created as follows (rgdal could also beused with the same output; GDAL prior to GDAL 3.3 built with GEOS, sothe BNA vector driver will use geometry tests: The BNA driver supportsreading of polygons with holes or lakes. It determines what is a hole ora lake only from geometrical analysis (inclusion, non-intersectiontests) and ignores completely the notion of polygon winding (whether thepolygon edges are described clockwise or counter-clockwise). GDAL mustbe built with GEOS enabled to make geometry test work.):
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")Since thespdep package was created,spatialweights objects have been constructed as lists with threecomponents and a few attributes, in old-style classlistwobjects. The first component of alistw object is annb object, a list ofn integer vectors, withat least a character vectorregion.id attribute withn unique values (like therow.names of adata.frame object);n is the number of spatialentities. Componenti of this list contains the integeridentifiers of the neighbours ofi as a sorted vector withno duplication and values in1:n; ifi has noneighbours, the component is a vector of length1 withvalue0L. Thenb object may contain anattribute indicating whether it is symmetric or not, that is whetheri is a neighbour ofj implies thatj is a neighbour ofi. Some neighbourdefinitions are symmetric by construction, such as contiguities ordistance thresholds, others are asymmetric, such ask-nearest neighbours. Thenb objectredundantly stores bothi-j andj-i links.
The second component of alistw object is a list ofn numeric vectors, each of the same length as thecorresponding non-zero vectors in thenbobject. These givethe values of the spatial weights for eachi-jneighbour pair. It is often the case that while the neighbours aresymmetric by construction, the weights are not, as for example whenweights arerow-standardised by dividing each row of inputweights by the count of neighbours or cardinality of the neighbour setofi. In thenb2listwfunction, it is alsopossible to pass through general weights, such as inverse distances,shares of boundary lengths and so on.
The third component of alistw object records thestyle of the weights as a character code, with"B" for binary weights taking values zero or one (only oneis recorded),"W" for row-standardised weights, and so on.In order to subsetlistw objects, knowledge of thestyle may be necessary.
First some housekeeping and setup to permit this vignette to be builtwhen packages are missing or out-of-date:
if (!suppressPackageStartupMessages(require(sf, quietly=TRUE))) { message("install the sf package") dothis <- FALSE}if (dothis) sf_extSoftVersion()## GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ ## "3.14.0" "3.11.3" "9.6.2" "true" "true" "9.6.2"Let us read the GPKG file with valid geometries in to ‘sf’ and ‘sp’objects:
NY8_sf <- st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData"), quiet=TRUE)table(st_is_valid(NY8_sf))## ## TRUE ## 281Here we first generate a queen contiguity nb object using the legacyspdep approach. This first either uses a pre-computed list of vectors ofprobable neighbours or finds intersecting bounding boxes internally.Then the points on the boundaries of each set of polygons making up anobservation are checked for a distance less than snap to any of thepoints of the set of polygons making up an observation included in theset of candidate neighbours. Because contiguity is symmetric, only i toj contiguities are tested. A queen contiguity is found as soon as onepoint matches, a rook contiguity as soon as two points match:
suppressPackageStartupMessages(library(spdep))reps <- 10eps <- sqrt(.Machine$double.eps)system.time(for(i in 1:reps) NY8_sf_1_nb <- poly2nb(NY8_sf, queen=TRUE, snap=eps))/reps## user system elapsed ## 0.1067 0.0026 0.1097Using spatial indices to check intersection of polygons is muchfaster than the legacy method in poly2nb. Fromspdep1.1-7, use is made of GEOS throughsf to find candidateneighbours whenfoundInBox=NULL, the default value. Becausecontiguity is symmetric by definition,foundInBox= onlyrequires intersections for higher indices, leading to a slight overheadto remove duplicates, asst_intersects() reports bothi j ansj i relationships. Asst_intersects() does not report whether neighbours arequeen orrook, a further step is needed todistinguish the two cases.
NY8_sf_1_nb## Neighbour list object:## Number of regions: 281 ## Number of nonzero links: 1632 ## Percentage nonzero weights: 2.066843 ## Average number of links: 5.807829spdep::poly2nb uses two heuristics, first to find candidateneighbours from intersecting polygons (st_intersects()),and second to use the symmetry of the relationship to halve the numberof remaining tests. This means that performance is linear in n, but withoverhead for identifying candidates, and back-filling symmetricneighbours. Further,spdep::poly2nb() stops searching forqueen contiguity as soon as the first neighbour point is found withinsnap distance (if not identical, which is tested first); a secondneighbour point indicates rook contiguities. For details of alternativesfor spherical geometries, see section @ref(spher-poly2nb) below.
Next, we explore a further possible source of differences inneighbour object reproduction, using the original version of the tractboundaries used in ASDAR, but with some invalid geometries as mentionedearlier (NY8_utm18.gpkg was created from the original ESRIShapefile used in ASDAR):
if (packageVersion("spData") >= "2.3.2") { NY8_sf_old <- sf::st_read(system.file("shapes/NY8_utm18.gpkg", package="spData"))} else { NY8_sf_old <- sf::st_read(system.file("shapes/NY8_utm18.shp", package="spData"))}## Reading layer `NY8_utm18' from data source ## `/home/rsb/lib/r_libs/spData/shapes/NY8_utm18.gpkg' using driver `GPKG'## Simple feature collection with 281 features and 17 fields## Geometry type: POLYGON## Dimension: XY## Bounding box: xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545## Projected CRS: WGS 84 / UTM zone 18Ntable(st_is_valid(NY8_sf_old))## ## FALSE TRUE ## 5 276We can see that there are a number of differences between theneighbour sets derived from the fully valid geometries and the olderpartly invalid set:
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)## [1] "Component 57: Numeric: lengths (4, 5) differ" ## [2] "Component 58: Numeric: lengths (5, 6) differ" ## [3] "Component 66: Numeric: lengths (7, 11) differ"## [4] "Component 73: Numeric: lengths (4, 5) differ" ## [5] "Component 260: Numeric: lengths (8, 9) differ"Usingst_make_valid() to make the geometries valid:
NY8_sf_old_val <- st_make_valid(NY8_sf_old, dist=0)table(st_is_valid(NY8_sf_old_val))## ## TRUE ## 281we also see that the geometry type of the geometry columnchanges:
class(st_geometry(NY8_sf_old))## [1] "sfc_POLYGON" "sfc"class(st_geometry(NY8_sf_old_val))## [1] "sfc_GEOMETRY" "sfc"and checking the"sfg" objects, two now have objects ofdifferent topological dimensions.
table(sapply(st_geometry(NY8_sf_old_val), function(x) class(x)[[2]]))## ## MULTIPOLYGON POLYGON ## 3 278This can be remedied usingst_collection_extract() toget the polygon objects:
NY8_sf_old_val <- st_collection_extract(NY8_sf_old_val, "POLYGON")table(sapply(st_geometry(NY8_sf_old_val), function(x) class(x)[[2]]))## ## MULTIPOLYGON ## 281However, in making the geometries valid, we change the geometries, sothe new sets of neighbours still differ from those made with the validgeometries in the same ways as before imposing validity:
try(NY8_sf_old_1_nb_val <- poly2nb(NY8_sf_old_val), silent = TRUE)all.equal(NY8_sf_old_1_nb_val, NY8_sf_1_nb, check.attributes=FALSE)## [1] "Component 57: Numeric: lengths (4, 5) differ" ## [2] "Component 58: Numeric: lengths (5, 6) differ" ## [3] "Component 66: Numeric: lengths (7, 11) differ"## [4] "Component 73: Numeric: lengths (4, 5) differ" ## [5] "Component 260: Numeric: lengths (8, 9) differ"The neighbour sets are the same for the old boundaries with orwithout imposing validity:
all.equal(NY8_sf_old_1_nb_val, NY8_sf_old_1_nb, check.attributes=FALSE)## [1] TRUEknearneigh() anddnearneigh() require pointsupport, so decisions must be taken about how to place the point in theareal object. We can usest_centroid() to get thecentroids, using theof_largest_polygon=TRUE argument tomake sure that the centroid is that of the largest polygon id theobservation is made up of more than one external ring:
NY8_ct_sf <- st_centroid(st_geometry(NY8_sf), of_largest_polygon=TRUE)orst_point_on_surface() which guarantees that the pointwill fall on the surface of a member polygon:
NY8_pos_sf <- st_point_on_surface(st_geometry(NY8_sf))or indeed taking the centre of the largest inscribed circle (thefunction returns a radius line segment, so we choose the central point,not the point on the circle):
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:
st_is_longlat(NY8_ct_sf)## [1] FALSEFrom this, we can check the graph-based neighbours (planarcoordinates only):
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## Warning in graph2nb(soi.graph(NY84_nb, NY8_ct_sf)): neighbour object has 2 sub-graphsNY86_nb <- graph2nb(gabrielneigh(NY8_ct_sf))NY87_nb <- graph2nb(relativeneigh(NY8_ct_sf))K-nearest neighbours use the coordinate matrices, and can handleGreat Circle distances, but this is not demonstrated here, as the dataset used is planar, in which casedbscan::kNN() in 2D or 3Dbuilding a kd-tree is used:
system.time(for (i in 1:reps) suppressWarnings(NY88_nb_sf <- knn2nb(knearneigh(NY8_ct_sf, k=1))))/reps## user system elapsed ## 0.0255 0.0009 0.0265Legacy code may be used omitting the kd-tree:
system.time(for (i in 1:reps) suppressWarnings(NY89_nb_sf <- knn2nb(knearneigh(NY8_ct_sf, k=1, use_kd_tree=FALSE))))/reps## user system elapsed ## 0.0265 0.0005 0.0271Distance neighbours need a threshold -nbdists shows themaximum distance to first nearest neighbour:
dsts <- unlist(nbdists(NY88_nb_sf, NY8_ct_sf))summary(dsts)## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 82.85 912.85 1801.11 3441.04 4461.26 17033.11max_1nn <- max(dsts)dnearneigh can also handle Great Circle distances, butthis is not demonstrated here, as the data set used is planar:
system.time(for (i in 1:reps) suppressWarnings(NY810_nb <- dnearneigh(NY8_ct_sf, d1=0, d2=0.75*max_1nn)))/reps## user system elapsed ## 0.0630 0.0016 0.0649By default, the function usesdbscan::frNN() to build akd-tree in 2D or 3D which is then used to find distance neighbours. Forsmall n, the argumentuse_kd_tree=FALSE may speed upcomputation a little by reverting to legacy code not building a kd-treefirst, but in general the differences are so small that the user willnot notice:
system.time(for (i in 1:reps) suppressWarnings(NY811_nb <- dnearneigh(NY8_ct_sf, d1=0, d2=0.75*max_1nn, use_kd_tree=FALSE)))/reps## user system elapsed ## 0.0490 0.0013 0.0504Spherical point-based neighbours may be found using Great Circledistances. These have been used for many years, but the switch ofsf 1.0-0 to uses2 by default hasopened up new opportunities where spatial indexing on the sphere mayhelp.
pts_ll <- st_transform(NY8_ct_sf, "OGC:CRS84")st_is_longlat(pts_ll)## [1] TRUEIf the input geometries are in geographical coordinates, andsf_use_s2() isTRUE,knearneigh()will use spatially indexed points ands2::s2_closest_edges() (seehttps://github.com/r-spatial/s2/issues/125#issuecomment-860107442)
(old_use_s2 <- sf_use_s2())## [1] TRUEand performs well with also with larger data sets:
sf_use_s2(TRUE)system.time(for (i in 1:reps) pts_ll1_nb <- knn2nb(knearneigh(pts_ll, k=6)))/reps## user system elapsed ## 0.0389 0.0000 0.0392For this smaller data set, the legacy approach without spatialindexing is adequate, but slows down as the number of observationsincreases:
sf_use_s2(FALSE)## Spherical geometry (s2) switched offsystem.time(for (i in 1:reps) pts_ll2_nb <- knn2nb(knearneigh(pts_ll, k=6)))/reps## user system elapsed ## 0.0276 0.0000 0.0277The WGS84 ellipsoid Great Circle distances differ a very little fromthes2 spherical distances, yielding output that herediverges for two tract centroids:
all.equal(pts_ll1_nb, pts_ll2_nb, check.attributes=FALSE)## [1] "Component 52: Mean relative difference: 1.466667" ## [2] "Component 124: Mean relative difference: 0.0251046"pts_ll1_nb[[52]]## [1] 15 38 48 49 50 53pts_ll2_nb[[52]]## [1] 37 38 48 49 50 53pts_ll1_nb[[124]]## [1] 117 122 123 125 133 134pts_ll2_nb[[124]]## [1] 116 117 123 125 133 134sf_use_s2(old_use_s2)## Spherical geometry (s2) switched onDistance neighbours are more problematic. Whilenbdists() works well withs2 sphericalcoordinates, none of the tried adaptations fordnearneigh()work adequately yet. An argumentuse_s2= is set to TRUE ifs2 > 1.0-7, usings2::s2_closest_edges() or the legacy brute-force approach,then only calculating distances fromi tojand copying those toj toi through symmetry.The distance metric is alway"km".
max_1nn_ll <- max(unlist(nbdists(knn2nb(knearneigh(pts_ll, k=1)), pts_ll)))## Warning in knn2nb(knearneigh(pts_ll, k = 1)): neighbour object has 62 sub-graphsargs(dnearneigh)## function (x, d1, d2, row.names = NULL, longlat = NULL, bounds = c("GE", ## "LE"), use_kd_tree = TRUE, symtest = FALSE, use_s2 = packageVersion("s2") > ## "1.0.7", k = 200, dwithin = TRUE) ## NULLIf we permits2 methods to run, without otherarguments set, ands2 > 1.0-7,s2::s2_dwithin_matrix() is run:
if (packageVersion("s2") > "1.0.7") { system.time(for (i in 1:(reps/5)) suppressWarnings(pts_ll3_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn_ll)))/(reps/5)}## user system elapsed ## 0.0715 0.0000 0.0715Alternatively, spherical distances can be used withdwithin=FALSE ands2::s2_closest_edges();although running in similar time,s2::s2_closest_edges()depends on the additionalk= argument, which, if mis-set,may miss valid neighbours:
system.time(for (i in 1:(reps/5)) suppressWarnings(pts_ll5_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn_ll, dwithin=FALSE)))/(reps/5)## user system elapsed ## 0.063 0.000 0.063if (packageVersion("s2") > "1.0.7") all.equal(pts_ll3_nb, pts_ll5_nb, check.attributes=FALSE)## [1] TRUEUsings2::s2_closest_edges() respectsd1 > 0 without requiring a second pass in R, so isfaster thans2::s2_dwithin_matrix():
if (packageVersion("s2") > "1.0.7") { system.time(for (i in 1:(reps/5)) suppressWarnings(pts_ll3a_nb <- dnearneigh(pts_ll, d1=5, d2=0.75*max_1nn_ll, dwithin=FALSE)))/(reps/5)}## user system elapsed ## 0.0530 0.0000 0.0525Usings2::s2_dwithin_matrix() requires a second pass,one for the lower bound, another for the upper bound, and a setdifference operation to find neighbours in the distance band:
if (packageVersion("s2") > "1.0.7") { system.time(for (i in 1:(reps/5)) suppressWarnings(pts_ll5a_nb <- dnearneigh(pts_ll, d1=5, d2=0.75*max_1nn_ll)))/(reps/5)}## user system elapsed ## 0.0980 0.0000 0.0985if (packageVersion("s2") > "1.0.7") all.equal(pts_ll3a_nb, pts_ll5a_nb, check.attributes=FALSE)## [1] TRUESettinguse_s2=FALSE falls back to the legacy version,which uses symmetry to reduce time:
system.time(for (i in 1:reps) suppressWarnings(pts_ll6_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn_ll, use_s2=FALSE)))/reps## user system elapsed ## 0.0374 0.0001 0.0376Minor differences may occur between the legacy ellipsoid ands2 spherical approaches:
all.equal(pts_ll5_nb, pts_ll6_nb, check.attributes=FALSE)## [1] "Component 20: Numeric: lengths (6, 5) differ" ## [2] "Component 28: Numeric: lengths (7, 6) differ" ## [3] "Component 112: Numeric: lengths (109, 108) differ"## [4] "Component 116: Numeric: lengths (109, 108) differ"## [5] "Component 122: Numeric: lengths (105, 106) differ"## [6] "Component 123: Numeric: lengths (108, 107) differ"## [7] "Component 130: Numeric: lengths (108, 109) differ"## [8] "Component 134: Numeric: lengths (106, 105) differ"## [9] "Component 158: Numeric: lengths (101, 102) differ"## [10] "Component 165: Numeric: lengths (101, 102) differ"## [11] "Component 168: Numeric: lengths (101, 102) differ"## [12] "Component 179: Numeric: lengths (89, 90) differ" ## [13] "Component 180: Numeric: lengths (96, 97) differ" ## [14] "Component 188: Numeric: lengths (46, 47) differ" ## [15] "Component 189: Numeric: lengths (55, 56) differ" ## [16] "Component 196: Numeric: lengths (47, 46) differ" ## [17] "Component 210: Numeric: lengths (106, 104) differ"## [18] "Component 226: Numeric: lengths (88, 87) differ" ## [19] "Component 229: Numeric: lengths (55, 53) differ" ## [20] "Component 235: Numeric: lengths (40, 39) differ" ## [21] "Component 237: Numeric: lengths (14, 15) differ" ## [22] "Component 245: Numeric: lengths (16, 15) differ"system.time(for (i in 1:reps) suppressWarnings(pts_ll6a_nb <- dnearneigh(pts_ll, d1=5, d2=0.75*max_1nn_ll, use_s2=FALSE)))/reps## user system elapsed ## 0.0302 0.0000 0.0303if (packageVersion("s2") > "1.0.7") all.equal(pts_ll5a_nb, pts_ll6a_nb, check.attributes=FALSE)## [1] "Component 20: Numeric: lengths (6, 5) differ" ## [2] "Component 28: Numeric: lengths (7, 6) differ" ## [3] "Component 112: Numeric: lengths (62, 61) differ" ## [4] "Component 113: Numeric: lengths (62, 63) differ" ## [5] "Component 116: Numeric: lengths (56, 55) differ" ## [6] "Component 119: Numeric: lengths (68, 69) differ" ## [7] "Component 122: Numeric: lengths (43, 44) differ" ## [8] "Component 123: Numeric: lengths (50, 49) differ" ## [9] "Component 128: Numeric: lengths (65, 64) differ" ## [10] "Component 130: Numeric: lengths (61, 63) differ" ## [11] "Component 132: Numeric: lengths (45, 46) differ" ## [12] "Component 134: Numeric: lengths (46, 45) differ" ## [13] "Component 136: Numeric: lengths (61, 62) differ" ## [14] "Component 147: Numeric: lengths (50, 51) differ" ## [15] "Component 154: Numeric: lengths (56, 57) differ" ## [16] "Component 158: Numeric: lengths (49, 50) differ" ## [17] "Component 165: Mean relative difference: 0.02823018"## [18] "Component 168: Numeric: lengths (54, 56) differ" ## [19] "Component 179: Numeric: lengths (77, 78) differ" ## [20] "Component 180: Numeric: lengths (85, 86) differ" ## [21] "Component 188: Numeric: lengths (39, 40) differ" ## [22] "Component 189: Numeric: lengths (48, 49) differ" ## [23] "Component 196: Numeric: lengths (45, 44) differ" ## [24] "Component 210: Numeric: lengths (68, 66) differ" ## [25] "Component 226: Numeric: lengths (82, 81) differ" ## [26] "Component 229: Numeric: lengths (48, 46) differ" ## [27] "Component 235: Numeric: lengths (38, 37) differ" ## [28] "Component 237: Numeric: lengths (14, 15) differ" ## [29] "Component 245: Numeric: lengths (15, 14) differ"It also turns out that whensf functions are used tofind contiguity neighbours,s2 spatial indexingfunctionality is accessed in finding candidate neighbours inintersecting geometries.
NY8_sf_ll <- st_transform(NY8_sf, "OGC:CRS84")st_is_longlat(NY8_sf_ll)## [1] TRUEThe timings are a little slower whenst_intersects()hands off geometry predicates tos2_intersects_matrix(),but the results are the same, and because spatial indexing is used, thisscales well for larger data sets:
sf_use_s2(TRUE)system.time(for (i in 1:reps) NY8_sf_1_nb_ll <- poly2nb(NY8_sf_ll, queen=TRUE, snap=eps))/reps## user system elapsed ## 0.1651 0.0018 0.1674all.equal(NY8_sf_1_nb, NY8_sf_1_nb_ll, check.attributes=FALSE)## [1] TRUE