5

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
3

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.

In addition, there is another thread similiar to your quesiton.

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. – moho wu Sep 17 '15 at 11:32
  • check the updated answer (use genrandompnts function in R) – Farid Cheraghi Sep 17 '15 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? – moho wu Sep 28 '15 at 20:43
  • 1
    "c:/data/Points_"+str(i) – Farid Cheraghi Sep 28 '15 at 20:45
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.

Add libraries and example 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)
  head(meuse@data)

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.

0

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:

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


## load data and convert to sf
columbus <- readOGR(system.file("shapes/columbus.shp", package="spData")[1]) %>%
              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
mapview::mapview(list(rdnmPts.rSSI, columbus))


# ... check number of random points
nrow(rdnmPts.rSSI)
# [1] 171

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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