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.