Merge pull request #164 from r-spatial/poly2nb_snap

Poly2nb snap

NAMESPACE

@@ -84,7 +84,8 @@ export(autocov_dist)
 
 export(set.VerboseOption, get.VerboseOption, set.ZeroPolicyOption,
         get.ZeroPolicyOption, get.SubgraphOption, set.SubgraphOption,
-        get.SubgraphCeiling, set.SubgraphCeiling)
+        get.SubgraphCeiling, set.SubgraphCeiling, get.NoNeighbourOption,
+        set.NoNeighbourOption)
 export(set.mcOption, get.mcOption, set.coresOption, get.coresOption,
         set.ClusterOption, get.ClusterOption)
 

NEWS.md

@@ -1,5 +1,9 @@
 # Version 1.3-6 (development)
 
+* #162 add option for no-neighbour checking for `poly2nb` - default report whether no-neighbour observations are present
+
+* #162 change the default `snap=` argument to `poly2nb` to 10mm 
+
 * Condition on forthcoming `tmap` 4
 
 * #160 handle `n.comp.nb` delay in `print.nb` and elsewhere when the total number of neighbours is large

R/AAA.R

@@ -15,6 +15,7 @@ assign("listw_is_CsparseMatrix", FALSE, envir = .spdepOptions)
 assign("cluster", NULL, envir = .spdepOptions)
 assign("report_nb_subgraphs", TRUE, envir = .spdepOptions)
 assign("nb_subgraphs_N+E", 100000L, envir = .spdepOptions)
+assign("report_nb_noneighs", TRUE, envir = .spdepOptions)
 setOldClass(c("listw"))
 
 .onLoad <- function(lib, pkg) {

R/poly2nb.R

@@ -1,12 +1,12 @@
-# Copyright 2001-2021 by Roger Bivand 
+# Copyright 2001-2024 by Roger Bivand 
 #
 #
 # Modified by Micah Altman 2010
 
 
 
-poly2nb <- function(pl, row.names=NULL, snap=sqrt(.Machine$double.eps),
-	queen=TRUE, useC=TRUE, foundInBox=NULL) {
+poly2nb <- function(pl, row.names=NULL, snap=NULL, queen=TRUE, useC=TRUE,
+        foundInBox=NULL) {
         verbose <- get("verbose", envir = .spdepOptions)
         .ptime_start <- proc.time()
         sf <- NULL
@@ -48,12 +48,38 @@ poly2nb <- function(pl, row.names=NULL, snap=sqrt(.Machine$double.eps),
 			else regid <- row.names
 		}
 	}
-        if (snap < 0) snap <- abs(snap)
-#        if (snap < .Machine$double.eps) {
-#            bsnap <- .Machine$double.eps
-#        } else { 
-#            bsnap <- snap
-#        }
+        if (!is.null(snap)) {
+            stopifnot(is.numeric(snap))
+            stopifnot(is.finite(snap))
+            stopifnot(length(snap) == 1L)
+            if (snap < 0) snap <- abs(snap)
+        } else {
+            if (sf) {
+                paras <- sf::st_crs(pl, parameters=TRUE)
+                if (length(paras) == 0L) {
+                    snap <- sqrt(.Machine$double.eps)
+                } else {
+                    if (paras$IsGeographic) {
+                        snap <- 9e-8
+                    } else {
+                        tenmm <- units::set_units(0.01, "metre")
+                        if (grepl("metre", paras$units_gdal)) {
+                            snap <- as.numeric(tenmm)
+                        } else {
+                            snap0 <- try(units::set_units(tenmm, paras$ud_unit,
+                                mode="standard"), silent=TRUE)
+                            if (inherits(snap0, "try-error")) {
+                                snap <- sqrt(.Machine$double.eps)
+                            } else {
+                                snap <- as.numeric(snap0)
+                            }
+                        }
+                    }
+                }
+            } else {
+                snap <- sqrt(.Machine$double.eps)
+            }
+        }
         vbsnap <- c(-snap, snap)
         if (verbose) cat("handle IDs:", (proc.time() - .ptime_start)[3], "\n")
         .ptime_start <- proc.time()
@@ -199,12 +225,17 @@ poly2nb <- function(pl, row.names=NULL, snap=sqrt(.Machine$double.eps),
 	attr(ans, "call") <- match.call()
 	if (queen) attr(ans, "type") <- "queen"
 	else attr(ans, "type") <- "rook"
+	attr(ans, "snap") <- snap
 	ans <- sym.attr.nb(ans)
-        NE <- n + sum(card(ans))
+        cans <- card(ans)
+        if (get.NoNeighbourOption()) {
+            if (any(cans == 0L)) warning("some observations have no neighbours;\nif this seems unexpected, try increasing the snap argument.")
+        }
+        NE <- n + sum(cans)
         if (get.SubgraphOption() && get.SubgraphCeiling() > NE) {
           ncomp <- n.comp.nb(ans)
           attr(ans, "ncomp") <- ncomp
-          if (ncomp$nc > 1) warning("neighbour object has ", ncomp$nc, " sub-graphs")
+          if (ncomp$nc > 1) warning("neighbour object has ", ncomp$nc, " sub-graphs;\nif this sub-graph count seems unexpected, try increasing the snap argument.")
         }
         if (verbose) cat("done:", (proc.time() - .ptime_start)[3], "\n")
         .ptime_start <- proc.time()

R/spChkOption.R

@@ -29,6 +29,17 @@ set.VerboseOption <- function(check) {
 	invisible(res)
 }
 
+get.NoNeighbourOption <- function() {
+	get("report_nb_noneighs", envir = .spdepOptions)
+}
+
+set.NoNeighbourOption <- function(check) {
+	if (!is.logical(check)) stop ("logical argument required")
+	res <- get("report_nb_noneighs", envir = .spdepOptions)
+	assign("report_nb_noneighs", check, envir = .spdepOptions)
+	invisible(res)
+}
+
 get.SubgraphOption <- function() {
 	get("report_nb_subgraphs", envir = .spdepOptions)
 }

inst/tinytest/test_poly2nb_snap.R

@@ -6,3 +6,9 @@ expect_true(isTRUE(all.equal(sum(card(poly2nb(p_sp))), 0L)))
 expect_true(isTRUE(all.equal(sum(card(poly2nb(p))), 0L)))
 expect_true(isTRUE(all.equal(sum(card(poly2nb(p_sp, snap=5))), 2L)))
 expect_true(isTRUE(all.equal(sum(card(poly2nb(p, snap=5))), 2L)))
+nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet=TRUE)
+nc_lcc_ft <- st_transform(nc, "EPSG:32019")
+nc_lcc_m <- st_transform(nc, "EPSG:3348")
+expect_equal(attr(poly2nb(nc), "snap"), 9e-8)
+expect_equal(attr(poly2nb(nc_lcc_m), "snap"), 0.01)
+expect_equal(attr(poly2nb(nc_lcc_ft), "snap"), 0.0328083333333333)

inst/tinytest/test_subgraph_warning.R

@@ -6,6 +6,7 @@ col_geoms[21] <- st_buffer(col_geoms[21], dist=-0.05)
 st_geometry(columbus) <- col_geoms
 set.SubgraphOption(FALSE)
 expect_false(get.SubgraphOption())
+set.NoNeighbourOption(FALSE)
 expect_silent(nb <- poly2nb(columbus))
 set.SubgraphOption(TRUE)
 expect_true(get.SubgraphOption())

man/poly2nb.Rd

@@ -5,14 +5,13 @@
 \description{
 The function builds a neighbours list based on regions with contiguous boundaries, that is sharing one or more boundary point. The current function is in part interpreted and may run slowly for many regions or detailed boundaries, but from 0.2-16 should not fail because of lack of memory when single polygons are built of very many border coordinates.}
 \usage{
-poly2nb(pl, row.names = NULL, snap=sqrt(.Machine$double.eps),
- queen=TRUE, useC=TRUE, foundInBox=NULL)
+poly2nb(pl, row.names = NULL, snap=NULL, queen=TRUE, useC=TRUE, foundInBox=NULL)
 }
 
 \arguments{
   \item{pl}{list of polygons of class extending \code{SpatialPolygons}, or an \code{sf} or \code{sfc} object containing non-empty (multi-)polygon objects}
   \item{row.names}{character vector of region ids to be added to the neighbours list as attribute \code{region.id}, default \code{seq(1, nrow(x))}; if \code{pl} has \code{row.names}, they are used instead of the default sequence.}
-  \item{snap}{boundary points less than \code{snap} distance apart are considered to indicate contiguity; used both to find candidate and actual neighbours for planar geometries, but only actual neighbours for spherical geometries, as spherical spatial indexing itself injects some fuzzyness.}
+  \item{snap}{boundary points less than \code{snap} distance apart are considered to indicate contiguity; used both to find candidate and actual neighbours for planar geometries, but only actual neighbours for spherical geometries, as spherical spatial indexing itself injects some fuzzyness. If not set, for all \code{SpatialPolygons} objects, the default is as before \code{sqrt(.Machine$double.eps)}, with this value also used for \code{sf} objects with no coordinate reference system. For \code{sf} objects with a defined coordinate reference system, the default value is \code{1e-7} for geographical coordinates (approximately 10mm), is 10mm where projected coordinates are in metre units, and is converted from 10mm to the same distance in the units of the coordinates. Should the conversion fail, \code{snap} reverts to \code{sqrt(.Machine$double.eps)}.}
   \item{queen}{if TRUE, a single shared boundary point meets the contiguity condition, if FALSE, more than one shared point is required; note that more than one shared boundary point does not necessarily mean a shared boundary line}
   \item{useC}{default TRUE, doing the work loop in C, may be set to false to revert to R code calling two C functions in an \code{n*k} work loop, where \code{k} is the average number of candidate neighbours}
   \item{foundInBox}{default NULL using R code or \code{st_intersects()} to generate candidate neighbours (using \code{snap=} if the geometries are not spherical); if not NULL (for legacy purposes) a list of length \code{(n-1)} with integer vectors of candidate neighbours \code{(j > i)} (as created by the \code{poly_findInBoxGEOS} function in \pkg{rgeos} for clean polygons)}

man/set.spChkOption.Rd

@@ -9,6 +9,8 @@
 \alias{get.SubgraphOption}
 \alias{set.SubgraphCeiling}
 \alias{get.SubgraphCeiling}
+\alias{set.NoNeighbourOption}
+\alias{get.NoNeighbourOption}
 \alias{set.ZeroPolicyOption}
 \alias{get.ZeroPolicyOption}
 \alias{set.listw_is_CsparseMatrix_Option}
@@ -31,6 +33,8 @@ set.SubgraphOption(check)
 get.SubgraphOption()
 set.SubgraphCeiling(value)
 get.SubgraphCeiling()
+set.NoNeighbourOption(check)
+get.NoNeighbourOption()
 set.listw_is_CsparseMatrix_Option(check)
 get.listw_is_CsparseMatrix_Option()
 }