2

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: enter image description here

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.

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

1 Answer 1

Reset to default
5

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

Improve this answer
4
  • 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.
    – Jeffrey Evans
    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!!
    – Spacedman
    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...
    – Spacedman
    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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.