Question: determination of neighbors when using snap

[JosiahParry]

spdep/man/poly2nb.Rd

Line 14 in 4cb8ba9

\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)}.}

I'm trying to understand a bit more clearly how neighbors are identified when using a snap tolerance. I understand that the spatial index bounding box search is increased by snap, however, when matching neighbors, I'm not too sure I understand how we determine if 2 neighbors are adjacent.

Is my reading of the code in

spdep/src/polypoly.c

Lines 182 to 206 in 4cb8ba9

	for (i=0; i<(nn-1); i++) {
	R_CheckUserInterrupt();
	li = Rf_length(VECTOR_ELT(i_findInBox, i));
	nrsi = NRS[i];
	for (j=0; j<li; j++) {
	jj = INTEGER_POINTER(VECTOR_ELT(i_findInBox, i))[j] - ROFFSET;
	jhit = spOverlapC(bb1[i], bb2[i], bb3[i], bb4[i], bb1[jj],
	bb2[jj], bb3[jj], bb4[jj]);
	if (jhit > 0) {
	khit = 0;
	nrsj = NRS[jj];
	if (nrsi > 0 && nrsj > 0){
	khit = polypolyC(&plx[cNRS[i]], &ply[cNRS[i]], nrsi,
	&plx[cNRS[jj]], &ply[cNRS[jj]], nrsj, Dsnap, crit+1L);
	}
	if (khit > crit) {
	card[i]++;
	card[jj]++;
	is[ii] = i;
	jjs[ii] = jj;
	ii++;
	}
	}
	}
	}

correct that if 1 vertex between the polygons is within [0, snap] for queen contiguity it is considered a neighbor and 2 vertices for rook contiguity?

[rsbivand]

@JosiahParry would it help if I created an instrumented branch, printing out values to the terminal?

There is quite a lot of legacy code - the parts used when useC=FALSE, some of which is untested since useC=TRUE became the default.

IIRC the bounding boxes (envelopes) are expanded by snap and used to check the foundInBox candidates, then the snap distance is used in

spdep/src/polypoly.c

Lines 85 to 110 in 4cb8ba9

	int polypolyC(double *px1, double *py1, int n1, double *px2, double *py2,
	int n2, double sn, int crit) {
	int i, j, k=0;
	double dist;
	double x1, x2, y1, y2, xd, yd;
	for (i=0; (i < n1) && (k < crit); i++) {
	R_CheckUserInterrupt();
	x1 = px1[i];
	y1 = py1[i];
	for (j=0; (j < n2) && (k < crit); j++) {
	x2 = px2[j];
	y2 = py2[j];
	xd = x1-x2;
	if (fabs(xd)>sn) { continue; }
	yd = y1-y2;
	if (fabs(yd)>sn) { continue; }
	dist = hypot(xd, yd);
	if (dist <= sn) {
	k++;
	}
	}
	}
	return(k);
	}

to compare the vertices of the pair of candidate polygons.

In the queen case, crit is 0, so the first i, j distance < snap sets k to 1 and exits the function with i, j as neighbours, if rook then crit is 1, so two "close enough" vertices set k to 2 and find the i. j neighbour.

I'd need to check more to see whether symmetry is used (if i neighbours j, then j must neighbour i).

[JosiahParry]

This is incredibly helpful. Thank you! If you'd like to make a branch with some messages along the way, that could be quite helpful.

I think the logic itself makes sense though. For queen contiguity we just need to be sure that there is a vertex somewhere close enough. And an edge consists to two points, so if we have two points within the tolerance they might be an edge

[rsbivand]

Yes, the foundInBox candidates appear to take account of symmetry,

spdep/src/polypoly.c

Lines 200 to 201 in 4cb8ba9

	is[ii] = i;
	jjs[ii] = jj;

It looks as though the foundInBox object respects snap for probably all but using s2 cases for sf objects - I need to check

spdep/R/poly2nb.R

Lines 146 to 153 in 4cb8ba9

	if (!is.na(st_is_longlat(pl)) &&
	st_is_longlat(pl) && sf_use_s2()) {
	if (packageVersion("sf") < "1.0.4") {
	fB1 <- st_intersects(pl, sparse=TRUE)
	} else {
	fB1 <- st_intersects(pl, sparse=TRUE, model="closed")
	}
	} else {

You are right that the rook criterion is two vertices within snap of each other. I did look at two proximate vertices, but ring direction makes this hard; a line segment with vertices A, B, C may be checked against a line segment with the vertices ordered C, B, A. So just counting two was easier at the time. Do you think matching proximate points matters? It would mean doing a topological analysis, there is something like this in rgrass::vect2neigh https://osgeo.github.io/rgrass/reference/read_VECT.html.

[JosiahParry]

In thinking about this a bit more—I think it would be somewhat straightforward to perform a check to see if two sequential vertices are present within snap to ensure that it is an edge. But making that topological comparison between the two features when separated by some distance is beyond my computational geometry abilities!

For determining if the two sequential vertices are within snap I think the following idea could work.

Check the first vertex. If it is within snap we can increment the counter k. If the first one is a hit, check if the last vertex is within snap. This way we don't have to worry about ring direction. If it is within snap early return, wonderful! Otherwise, continue iterating through the remaining vertices. If the next vertex is not within snap, reset k = 0. If k = 2 then we can return and indicate that there is a symmetric neighbors in

 is[ii] = i; 
 jjs[ii] = jj;