# Calculating the cross-nearest-neighbor distance in spatstat

I wish to calculate the cross-nearest-neighbor distance between two point patterns in `spatstat`. According to Pr. Sadahiro, the cross-nearest-neighbor distance is defined as:

My inputs are `sf` datasets, which I convert to ppp format, before entering into a formula for distance.

``````#Create the sf datasets
library(sf)
#10-point dataset
pub <- structure(list(X = c(959207.877070254, 959660.734838225, 951483.685462513,  951527.767554883, 958310.673042469, 950492.05212104, 959660.734838225,  959207.877070254, 960500.020456073, 959660.734838225), Y = c(1944457.42827898,  1955543.76027363, 1939982.16629396, 1940216.55143212, 1954704.68186897,  1951434.68524296, 1955543.76027363, 1944457.42827898, 1955292.64874361,  1955543.76027363), geometry = structure(list(structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(951483.685462513,  1939982.16629396), class = c("XY", "POINT", "sfg")), structure(c(951527.767554883,  1940216.55143212), class = c("XY", "POINT", "sfg")), structure(c(958310.673042469,  1954704.68186897), class = c("XY", "POINT", "sfg")), structure(c(950492.05212104,  1951434.68524296), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg")), structure(c(959207.877070254,  1944457.42827898), class = c("XY", "POINT", "sfg")), structure(c(960500.020456073,  1955292.64874361), class = c("XY", "POINT", "sfg")), structure(c(959660.734838225,  1955543.76027363), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 950492.05212104,  ymin = 1939982.16629396, xmax = 960500.020456073, ymax = 1955543.76027363 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS[\"Korea 2000 / Unified CS\",\n    BASEGEOGCRS[\"Korea 2000\",\n        DATUM[\"Geocentric datum of Korea\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4737]],\n    CONVERSION[\"Korea Unified Belt\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",38,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",127.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",1000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",2000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Korea, Republic of (South Korea)\"],\n        BBOX[28.6,122.71,40.27,134.28]],\n    ID[\"EPSG\",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(4177L,  15721L, 21365L, 21973L, 24836L, 59359L, 66313L, 70379L, 83277L,  90828L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))
#21-point dataset
pat <- structure(list(X = c(950037.869142169, 952809.658320316, 957446.529265191,  957446.529265191, 951548.096552647, 953896.66691363, 906343.850625855,  959380.829011949, 951936.426429872, 950288.722667828, 959109.852180828,  958181.105739686, 959968.373867043, 959380.829011949, 957446.529265191,  959380.829011949, 932679.405187587, 953993.414749656, 960049.499276811,  952735.358695265, 958624.334501739), Y = c(1955782.34759634,  1946446.42376891, 1954520.57051429, 1954520.57051429, 1942270.95246119,  1951222.182335, 1940637.27489082, 1952879.9377109, 1932938.1054862,  1953450.8956205, 1951649.75074957, 1960119.73774404, 1946530.91451559,  1952879.9377109, 1954520.57051429, 1952879.9377109, 1948892.87536131,  1952741.63981237, 1963350.22042542, 1948920.98629227, 1942820.79162891 ), geometry = structure(list(structure(c(950037.869142169, 1955782.34759634 ), class = c("XY", "POINT", "sfg")), structure(c(952809.658320316,  1946446.42376891), class = c("XY", "POINT", "sfg")), structure(c(957446.529265191,  1954520.57051429), class = c("XY", "POINT", "sfg")), structure(c(957446.529265191,  1954520.57051429), class = c("XY", "POINT", "sfg")), structure(c(951548.096552647,  1942270.95246119), class = c("XY", "POINT", "sfg")), structure(c(953896.66691363,  1951222.182335), class = c("XY", "POINT", "sfg")), structure(c(906343.850625855,  1940637.27489082), class = c("XY", "POINT", "sfg")), structure(c(959380.829011949,  1952879.9377109), class = c("XY", "POINT", "sfg")), structure(c(951936.426429872,  1932938.1054862), class = c("XY", "POINT", "sfg")), structure(c(950288.722667828,  1953450.8956205), class = c("XY", "POINT", "sfg")), structure(c(959109.852180828,  1951649.75074957), class = c("XY", "POINT", "sfg")), structure(c(958181.105739686,  1960119.73774404), class = c("XY", "POINT", "sfg")), structure(c(959968.373867043,  1946530.91451559), class = c("XY", "POINT", "sfg")), structure(c(959380.829011949,  1952879.9377109), class = c("XY", "POINT", "sfg")), structure(c(957446.529265191,  1954520.57051429), class = c("XY", "POINT", "sfg")), structure(c(959380.829011949,  1952879.9377109), class = c("XY", "POINT", "sfg")), structure(c(932679.405187587,  1948892.87536131), class = c("XY", "POINT", "sfg")), structure(c(953993.414749656,  1952741.63981237), class = c("XY", "POINT", "sfg")), structure(c(960049.499276811,  1963350.22042542), class = c("XY", "POINT", "sfg")), structure(c(952735.358695265,  1948920.98629227), class = c("XY", "POINT", "sfg")), structure(c(958624.334501739,  1942820.79162891), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",  "sfc"), precision = 0, bbox = structure(c(xmin = 906343.850625855,  ymin = 1932938.1054862, xmax = 960049.499276811, ymax = 1963350.22042542 ), class = "bbox"), crs = structure(list(input = "EPSG:5179",      wkt = "PROJCRS[\"Korea 2000 / Unified CS\",\n    BASEGEOGCRS[\"Korea 2000\",\n        DATUM[\"Geocentric datum of Korea\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4737]],\n    CONVERSION[\"Korea Unified Belt\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",38,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",127.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",1000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",2000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (X)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (Y)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Korea, Republic of (South Korea)\"],\n        BBOX[28.6,122.71,40.27,134.28]],\n    ID[\"EPSG\",5179]]"), class = "crs"), n_empty = 0L)), row.names = c(11459L,  13177L, 17483L, 17484L, 17491L, 17494L, 1099074L, 1099643L, 1100354L,  1100723L, 1100844L, 2427603L, 2427604L, 2427605L, 2427606L, 2427607L,  2427608L, 2427609L, 2427610L, 2427611L, 2427612L), class = c("sf",  "data.frame"), sf_column = "geometry", agr = structure(c(X = NA_integer_,  Y = NA_integer_), .Label = c("constant", "aggregate", "identity" ), class = "factor"))

#Convert to point patterns
library(spatstat)
pub <- unmark(as.ppp(pub))
pat <- unmark(as.ppp(pat))

#Calculate cross-nearest-neighbor distance
(sum(nncross(pub, pat)\$dist) + sum(nncross(pat, pub)\$dist))/(npoints(pub) + npoints(pat))
[1] 4648.36
``````

Is this correct? Is there a better way of doing it?

I feel uncertain about the conversion from sf to ppp and the relative integrity of the coordinates between the two patterns.

• I've added the R tag which should help people find this. Feb 19, 2021 at 9:12

I can get the same result using only `sf` and `FNN` packages:

Starting with the `sf` versions of `pub` and `pat`:

``````> library(FNN)
> nnpat = knnx.dist(st_coordinates(pat), st_coordinates(pub), k=1)[,1]
> nnpub = knnx.dist(st_coordinates(pub), st_coordinates(pat), k=1)[,1]
> (sum(nnpub)+sum(nnpat))/(nrow(pub)+nrow(pat))
[1] 4648.36
``````

I think that confirms both answers are (probably) correct, but note that this is all being done in cartesian coordinates.

As to whether which is "better"? That will depend on your exact problem and you'd have to compare some examples.

Here's two functions that compute V given two `ppp` objects:

``````Vspat = function(pub,pat){
pub = unmark(pub)
pat = unmark(pat)
(sum(nncross(pub, pat)\$dist) + sum(nncross(pat, pub)\$dist))/(npoints(pub) + npoints(pat))
}
``````

or two 2-column matrices (which is what `st_coordinates` gets from spatial points):

``````Vmatrix = function(pub,pat){
nnpat = knnx.dist(pat, pub, k=1, algorithm="kd_tree")[,1]
nnpub = knnx.dist(pub, pat, k=1, algorith="kd_tree")[,1]
(sum(nnpub)+sum(nnpat))/(nrow(pub)+nrow(pat))
}
``````

Generate some test data with random normal density returning both spatstat and plain matrix form:

``````btestdata <- function(n1=100, n2=100){
m1 = cbind(rnorm(n1), rnorm(n1))
m2 = cbind(rnorm(n2), rnorm(n2))

sf1 = st_as_sf(data.frame(m1), coords=1:2)
sf2 = st_as_sf(data.frame(m2), coords=1:2)

p1 = as.ppp(sf1)
p2 = as.ppp(sf2)

return(
list(
spat=list(p1=p1, p2=p2),
matrix = list(m1=m1, m2=m2)
)
)
}
``````

Check our functions produce the same answer on a data set:

``````> d = btestdata(100,100)
> Vmatrix(d\$matrix\$m1, d\$matrix\$m2)
[1] 0.2525389
> Vspat(d\$spat\$p1, d\$spat\$p2)
[1] 0.2525389
``````

Here's our benchmark test using `microbenchmark` package:

``````bench <- function(datalist){
return(microbenchmark(
matrix = Vmatrix(datalist\$matrix\$m1, datalist\$matrix\$m2),
spatstat = Vspat(datalist\$spat\$p1, datalist\$spat\$p2)
)
)
}
``````

for 10 points in each class:

``````> d = btestdata(10,10)
> bench(d)
Unit: microseconds
expr      min       lq     mean    median        uq      max neval
matrix   80.930   91.138  122.512  109.9305  140.3835  290.656   100
spatstat 3359.651 3519.349 4014.001 3727.6450 4223.4865 8256.213   100
``````

the matrix algorithm wins bigly.

For 100 points in each class:

``````> bench(d)
Unit: microseconds
expr      min       lq      mean    median       uq      max neval
matrix  201.669  211.656  305.8227  229.1925  267.937 6684.524   100
spatstat 3307.209 3519.107 3823.0726 3644.3945 3827.220 9534.207   100
``````

The matrix algorithm still wins pretty clearly.

Up to 4000 points in each class:

``````> d = btestdata(4000,4000)
> bench(d)
Unit: milliseconds
expr      min       lq     mean   median       uq      max neval
matrix 6.220665 6.650980 6.932406 6.828636 7.171964 10.83241   100
spatstat 5.937900 6.337493 6.762920 6.591006 6.908898 10.19699   100
``````

spatstat wins! Let's do 100,000 points:

``````> d = btestdata(100000,100000)
> bench(d)
Unit: milliseconds
expr      min       lq     mean   median       uq      max neval
matrix 305.2815 315.6360 346.5156 324.9187 373.6829 504.6390   100
spatstat 201.5483 207.0548 218.3515 211.2062 219.0227 341.9087   100
>
``````

spatstat still winning!

So if you have to do a lot of small point-count cases, use the matrix algorithm from `FNN` package, otherwise if you've got large point-count cases, the overhead in converting to `ppp` might win (note there's no conversion overhead here, I'm benchmarking by feeding the algorithms with the food they can digest easiest).

• The advantage of @Spacedman answer is that, if so inclined, you could easily extend this to a multivariate problem. Oh, and there are the massive improvements in speed. Feb 19, 2021 at 16:16
• `spatstat::nncross` does seem to do some sort of indexing to make the algorithm better than `O(n x m)` efficient, but the `FNN` package has super-optimised algorithms. See now I have to benchmark these things!! Feb 19, 2021 at 16:22
• Hmmm for anything less than 3000 pts vs 3000 pts `FNN` is faster, but then `spatstat` wins. I suspect maybe FNN is optimised for higher dimensional problems... Feb 19, 2021 at 16:41
• I had marked this as a solution but unmarked it because it is described as only "probably" correct.
– syre
Mar 1, 2021 at 5:10