# Randomly sampling points in R with minimum distance constraint?

I'm trying to randomly select a number of points within my data frame (example data below) with a constraint that the minimum distance between the selected point must be greater than a certain distance. I managed to do the randomly selection bit using the sample function in R, but I can't figure out how to add the constraint bit into my code. I suppose this must involve using some spatial analysis package in R but I haven't got a clue where to start.

I know ArcGIS has a tool called Create Random Points which can specify the minimum distance between points. But my situation requires a larger number of repeated sampling, thus making me feel doing this in R would be much easier because it can be incorporated with a loop.

Example data:

grid_index x        y
grid_168   323012.5 674187.5
grid_169   323012.5 674212.5
grid_292   323037.5 672287.5
grid_293   323037.5 672312.5
grid_368   323037.5 674187.5
grid_369   323037.5 674212.5

If I understand you correctly, you want to draw a distance-constrained random sample from your data for each observation in the data. This is akin to a K nearest neighbor analysis.

Here is an example workflow that will create a kNN random sample, using a minimum distance constraint, and add the corresponding rowname back to your data.

library(sp)
data(meuse)
coordinates(meuse) <- ~x+y

Calculate a distance matrix using spDists

dmat <- spDists(meuse)

Define minimum sample distance and set to NA in distance matrix. Here is where you would create any type of constraint say, a distance range.

min.dist <- 500
dmat[dmat <= min.dist] <- NA

Here we iterate through each row in the distance matrix and select a random sample != NA. The "samples" object is a data.frame where ID is the rownames of the source object and kNN is the rowname of the nearest neighbor. Note; there is some NA handling added just in case no neighbor is found, which could happen with distance constraints.

samples <- data.frame(ID=rownames(meuse@data), kNN=NA)
for(i in 1:nrow(dmat) ) {
x <- as.vector( dmat[,i] )
names(x) <- samples\$ID
x <- x[!is.na(x)]
if(!length(x) == 0) {
samples[i,][2] <- names(x)[sample(1:length(x), 1)]
} else {
samples[i,][2] <- NA
}
}

We can then add the kNN column, containing the rownames of the nearest neighbor, to the original data.

meuse@data <- data.frame(meuse@data, kNN=samples\$kNN)

We could also subset the unique nearest neighbor observations.

meuse.sub <- meuse[which(rownames(meuse@data) %in% unique(samples\$kNN)),]

There are much more elegant ways to perform this analysis but this workflow gets the general idea across. I would recommend taking a hard look at the spdep library and dnearneigh or knearneigh functions for a more advanced solution.

• Very helpful and very well described. Much appreciated.
– Gray
Oct 21, 2020 at 1:49

You can do this either with R or ArcGIS independently.

With ArcGIS, first create a feature class (e.g. shape file) from your grid coordinates. Then use this grid feature class as constraining_extent parameter of "Create Random Points" tool.

The only coding you have to do is to put that tool in a loop which can be achieved via Model Builder or Arcpy.

here is a sample (100 iteration):

import arcpy

for i in range(100):
print i
arcpy.CreateRandomPoints_management("c:/data/project", "samplepoints", "c:/data/studyarea.shp", "", 500, "", "POINT", "")

For R, use genrandompnts from spatialecology website. This tool is similiar to ArcGis "Create Random Points" tool.

How to create random points outside polygons?

• Thanks for the comments. Are you aware of any equivalent functions in R. Because my previous codes were written in R and I'm not very familiar with Python, which could make my further customization of the code difficult. In addition what I really need for the output is just a data frame containing the coordinate rather than a shapefile produced by the ArcGIS. Cheers. Sep 17, 2015 at 11:32
• check the updated answer (use genrandompnts function in R) Sep 17, 2015 at 19:03
• I'm reverting back to doing this in ArcGIS, because the genrandompnts only allows selecting from polygon but doesn't allow selecting from a set of points. May I ask in Python, how do I write the results in different files with names differentiated by the loop time? something like paste("c:/data/", "i") in R? Sep 28, 2015 at 20:43
• "c:/data/Points_"+str(i) Sep 28, 2015 at 20:45

If one is interested in sampling points with a distance constraint for each polygon, in the meanwhile, there are two nice and fast possibilities: (1) using QGIS "random points inside polygons" through RQGIS-package , or (2) using spatstat::rSSI-function.

In the following examples for (1) RQGIS and (2) spatstat::rSSI:

if(!require("pacman")) install.packages("pacman")
pacman::p_load(sf, sp, rgdal, dplyr, mapview, spatstat, maptools, devtools)

## load data and convert to sf
sf::st_as_sf(.)

## start random sampling with distance constraint

# # # # # # # # # # # # # # # # # # # #
# (1) ... using RQGIS
# # # # # # # # # # # # # # # # # # # #

devtools::install_github("jannes-m/RQGIS")
library("RQGIS")

# ... open QGIS tunnel
RQGIS::open_app() # QGIS must be installed

# ... find suitable algorithm
RQGIS::find_algorithms(search_term = "random")
# [6] "Random points inside polygons (fixed)
# ---------------->qgis:randompointsinsidepolygonsfixed"

# [7] "Random points inside polygons (variable)
# ------------->qgis:randompointsinsidepolygonsvariable"

# ... check usage
RQGIS::get_usage(alg = "qgis:randompointsinsidepolygonsfixed")
RQGIS::get_args_man(alg = "qgis:randompointsinsidepolygonsfixed")

# ... process random points with minimum distance (0.25 degree)
#     using a fixed maximum number (10)
rdnmPts.RQGIS <- RQGIS::run_qgis(alg = "qgis:randompointsinsidepolygonsfixed",
show_output_paths = TRUE,
load_output = TRUE, params = list(
VECTOR =  columbus, MIN_DISTANCE = "0.25",
VALUE = "10",OUTPUT = "rndm_pts.shp"))

# ... take a look on the result
mapview::mapview(list(rdnmPts.RQGIS, columbus))

# ... check number of random points
nrow(rdnmPts.RQGIS)
# [1] 187

# # # # # # # # # # # # # # # # # # # #
# (2) ... using spatstat::rSSI ------------------------------------
# # # # # # # # # # # # # # # # # # # #

# spatstat::rSSI uses a special format input. Therefore, a function is created
# to transform the simple feature to owin-format.

# init function
genRandomPtsDist <- function(x, seed = 123, dist = 10, n = Inf,
maxit = 100, quiet = TRUE, ...)
{

# get start time of process
process.time.start <- proc.time()

# get crs
crs <- sf::st_crs(x = x)

# convert simple feature to spatial polygons
x.sp <- x %>% as(., "Spatial") %>%  as(., "SpatialPolygons")

# convert to owin object
x.owin <- x.sp %>%
slot(., "polygons") %>%
lapply(X = ., FUN = function(x){sp::SpatialPolygons(list(x))}) %>%
lapply(X = ., FUN = spatstat::as.owin)

# generate random sampling with distant constraint (can be parallelized)
pts.ppp <- lapply(X = 1:length(x.owin), FUN = function(i, x.owin, r, n, quiet,
seed, maxit, ...)
{
if(quiet == FALSE) cat("Run ", i, " of ", length(x.owin), "\n")
set.seed(seed)
spatstat::rSSI(r = r, n = n, giveup = maxit, win = x.owin[[i]], ...)
}, quiet = quiet, x.owin = x.owin, r = dist, n = n, seed = seed, maxit = maxit, ...)

# back-conversion to simple feature
pts.sf <- pts.ppp %>%
lapply(X = ., FUN = function(x) sf::st_as_sfc(as(x, "SpatialPoints"))) %>%
do.call(c, .) %>%
sf::st_sf(., crs = crs)

# get intersected items
pts.inter.x <- sf::st_intersects(x = pts.sf, y = x) %>% unlist

if(length(pts.inter.x) != nrow(pts.sf))
{
warning("Some sample points are outside a polygon")
} else{
pts.sf\$In <- pts.inter.x
}

# get time of process
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat(paste0("------ Run of genRandomPtsDist: " ,
round(x = process.time.run["elapsed"][[1]]/60, digits = 3), " Minutes ------\n"))

return(pts.sf)
} # end of function

## Using the defined function, now one can generate a random sample.
# ... process random points with minimum distance (0.25 degree)
#     using a fixed maximum number (10)
rdnmPts.rSSI <- genRandomPtsDist(x = columbus, dist = 0.25, n = 10, quiet = FALSE)

# ... take a look on the result