I am trying to calculate a NNA out of two Shapefiles with location data (points). The points are trees in a study area. Unfortuneatly, the result, sum, is always 0. But it should be at least 20. There are no occuring errors.
There are no identical points and all points are within the study area.
library(spatstat)
# Load the data
d_trees <- sf::st_read("digitised_trees.shp")
s_trees <- sf::st_read("simulated_trees.shp")
area <- sf::st_read("study_area.shp")
area_subset <- area
dput(area_subset)
structure(list(Id = 0L, geometry = structure(list(structure(list(
structure(c(673750, 673750, 672550.000000001, 672550.000000001,
673750, 5189650, 5188850, 5188850, 5189650, 5189650), dim = c(5L,
2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
input = "WGS 84 / UTM zone 32N", wkt = "PROJCRS["WGS 84 / UTM zone 32N",n BASEGEOGCRS["WGS 84",n DATUM["World Geodetic System 1984",n ELLIPSOID["WGS 84",6378137,298.257223563,n LENGTHUNIT["metre",1]]],n PRIMEM["Greenwich",0,n ANGLEUNIT["degree",0.0174532925199433]],n ID["EPSG",4326]],n CONVERSION["UTM zone 32N",n METHOD["Transverse Mercator",n ID["EPSG",9807]],n PARAMETER["Latitude of natural origin",0,n ANGLEUNIT["Degree",0.0174532925199433],n ID["EPSG",8801]],n PARAMETER["Longitude of natural origin",9,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",500000,n LENGTHUNIT["metre",1],n ID["EPSG",8806]],n PARAMETER["False northing",0,n LENGTHUNIT["metre",1],n ID["EPSG",8807]]],n CS[Cartesian,2],n AXIS["(E)",east,n ORDER[1],n LENGTHUNIT["metre",1]],n AXIS["(N)",north,n ORDER[2],n LENGTHUNIT["metre",1]],n ID["EPSG",32632]]"), class = "crs"), class = c("sfc_POLYGON",
"sfc"), precision = 0, bbox = structure(c(xmin = 672550.000000001,
ymin = 5188850, xmax = 673750, ymax = 5189650), class = "bbox"))), row.names = c(NA,
-1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(Id = NA_integer_), class = "factor", levels = c("constant",
"aggregate", "identity")))
# Transform the data to the projected CRS
d_trees <- sf::st_transform(d_trees, 32632)
s_trees <- sf::st_transform(s_trees, 32632)
area_subset <- sf::st_transform(area, 32632)
# Extract the coordinates from the sf objects
d_coords <- sf::st_coordinates(d_trees)
s_coords <- sf::st_coordinates(s_trees)
# Output the coordinates using dput
dput(d_coords)
structure(c(672791.910571324, 672799.480667902, 672794.43393685,
672680.146506565, 672694.445577879, 672651.548363936, 672670.315895036,
672662.745798458, 672682.254178773, 672670.946736418, 672695.549550297,
672678.516832996, 672711.320584835, 672704.381329638, 672732.611481461,
672723.148860738, 672730.718957317, 672741.443260803, 672692.868474426,
672847.010291809, 5188901.84124404, 5188865.25244391, 5188880.39263707,
5188881.75946006, 5188860.73141401, 5188869.14263243, 5188914.77349236,
5188954.5164994, 5188954.91661639, 5188948.83892696, 5188964.6099615,
5188931.80620966, 5188967.29103737, 5188951.99313387, 5188961.1403339,
5188967.44874772, 5188978.1730512, 5188980.69641673, 5188981.95809949,
5188861.64647692), dim = c(20L, 2L), dimnames = list(NULL, c("X",
"Y")))
dput(s_coords)
structure(c(673656.645769361, 673025.118431594, 673453.162168716,
673440.619588101, 673589.648795813, 673448.33074397, 673558.477587886,
673140.127848215, 673191.602010456, 673079.423488968, 673225.014156644,
673502.646707603, 673514.304039728, 673589.559533666, 673535.402585275,
673423.8894745, 673665.262841519, 673164.143936007, 673196.893453529,
673404.373844275, 5188907.16103378, 5188931.72983531, 5188928.38909054,
5189028.38824619, 5189023.50717526, 5188884.22643243, 5188944.2819144,
5188976.35201315, 5189012.64219437, 5188931.033843, 5188878.58445754,
5189029.38962503, 5188851.62215933, 5188876.38281282, 5188943.02218533,
5188914.60512405, 5188933.84008721, 5188852.54412112, 5189024.96523525,
5189075.34562712), dim = c(20L, 2L), dimnames = list(NULL, c("X",
"Y")))
# Convert to ppp objects and calculate the mean nearest-neighbor distances (NND)
d_trees_ppp <- spatstat.geom::as.ppp(d_coords, area)
d_nnd <- spatstat.geom::nndist(d_trees_ppp)
d_mnnd_emp <- mean(d_nnd)
s_trees_ppp <- spatstat.geom::as.ppp(s_coords, area)
s_nnd <- spatstat.geom::nndist(s_trees_ppp)
s_mnnd_emp <- mean(s_nnd)
# Calculate theoretical MNNDs for d_trees
d_mnnd_theo <- numeric(100)
for (i in 1:100) {
d_area_sample <- sf::st_sample(area_subset, 4145)
d_area_ppp <- spatstat.geom::as.ppp(d_area_sample)
d_nnd_sample <- spatstat.geom::nndist(d_area_ppp)
d_mnnd_theo[i] <- mean(d_nnd_sample)
}
sum(d_mnnd_theo < d_mnnd_emp)
# Calculate theoretical MNNDs for s_trees
s_mnnd_theo <- numeric(100)
for (i in 1:100) {
s_area_sample <- sf::st_sample(area_subset, 3879)
s_area_ppp <- spatstat.geom::as.ppp(s_area_sample)
s_nnd_sample <- spatstat.geom::nndist(s_area_ppp)
s_mnnd_theo[i] <- mean(s_nnd_sample)
}
sum(s_mnnd_theo < s_mnnd_emp)
New contributor
Regggina is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
3