Saltar al contenido

¿Muestreo aleatorio de puntos en R con restricción de distancia mínima?

Solución:

Si lo entiendo correctamente, desea extraer una muestra aleatoria restringida por la distancia de sus datos para cada observación en los datos. Esto es similar a un análisis de K vecino más cercano.

Aquí hay un flujo de trabajo de ejemplo que creará una muestra aleatoria kNN, usando una restricción de distancia mínima, y ​​agregará el nombre de fila correspondiente a sus datos.

Agregar bibliotecas y datos de ejemplo

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

Calcular una matriz de distancia usando spDists

 dmat <- spDists(meuse)

Defina la distancia de muestra mínima y ajústela a NA en la matriz de distancias. Aquí es donde crearía cualquier tipo de restricción, digamos, un rango de distancia.

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

Aquí iteramos a través de cada fila en la matriz de distancias y seleccionamos una muestra aleatoria! = NA. El objeto “samples” es un data.frame donde ID son los nombres de fila del objeto de origen y kNN es el nombre de fila del vecino más cercano. Nota; se agrega algo de manejo de NA en caso de que no se encuentre un vecino, lo que podría suceder con restricciones de distancia.

samples <- data.frame(ID=rownames([email protected]), 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
    }   
  }

Luego podemos agregar la columna kNN, que contiene los nombres de fila del vecino más cercano, a los datos originales.

[email protected] <- data.frame([email protected], kNN=samples$kNN)
  head([email protected])

También podríamos crear un subconjunto de las observaciones únicas del vecino más cercano.

meuse.sub <- meuse[which(rownames([email protected]) %in% unique(samples$kNN)),]

Hay formas mucho más elegantes de realizar este análisis, pero este flujo de trabajo transmite la idea general. Recomendaría echar un vistazo a la biblioteca spdep y las funciones dnearneigh o knearneigh para una solución más avanzada.

Puede hacer esto con R o ArcGIS de forma independiente.

Con ArcGIS, primero cree una clase de entidad (por ejemplo, un archivo de forma) a partir de las coordenadas de su cuadrícula. Luego, use esta clase de entidad de cuadrícula como parámetro constraining_extent de la herramienta “Crear puntos aleatorios”.

La única codificación que tiene que hacer es poner esa herramienta en un bucle que se puede lograr a través de Model Builder o Arcpy.

aquí hay una muestra (100 iteraciones):

import arcpy

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

Para R, use genrandompnts del sitio web de Spatialecology. Esta herramienta es similar a la herramienta “Crear puntos aleatorios” de ArcGis.

Además, hay otro hilo similar a su pregunta.

¿Cómo crear puntos aleatorios fuera de polígonos?

Si uno está interesado en puntos de muestreo con una restricción de distancia para cada polígono, mientras tanto, hay dos posibilidades agradables y rápidas: (1) usar QGIS “puntos aleatorios dentro de polígonos” a través del paquete RQGIS, o (2) usar Spatstat: : función rSSI.

En los siguientes ejemplos para (1) RQGIS y (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
¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)


Tags : / /

Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *