¿Crear una gran cantidad de puntos aleatorios en el ráster binario?

Quiero crear un conjunto de datos vectoriales de puntos de 10000 puntos (o más grandes) dentro de un ráster binario, donde los puntos deben estar restringidos a áreas donde el valor del ráster es 1.

Intenté los siguientes pasos.

1. Poligonalizar ráster
2. QGIS: Vector -> Herramientas de investigación -> Puntos aleatorios

Esto funciona bien hasta 2000 puntos, pero todo lo anterior solo hace que QGIS se bloquee.

¿Hay alguna manera de crear un conjunto de datos vectoriales con una gran cantidad de entidades de puntos restringidas por un ráster binario (o la versión poligonalizada del mismo)?

Las siguientes herramientas están a mi disposición, clasificadas de más a menos favorables: QGIS, Python, R, ArcGIS

Esto es lo que busco, solo con 10 veces las características de punto.

Aquí hay una manera en R:

Haga un ráster de prueba, 20x30 celdas, haga 1/10 de las celdas establecidas en 1, trace:

> require(raster)
> m = raster(nrow=20, ncol=30)
> m[] = as.numeric(runif(20*30)>.9)
> plot(m)

Para un ráster existente en un archivo, por ejemplo, un geoTIFF, puede hacer lo siguiente:

> m = raster("mydata.tif")

Ahora obtenga una matriz de las coordenadas xy de las celdas 1, trace esos puntos y veremos que tenemos centros de celdas:

> ones = xyFromCell(m,1:prod(dim(m)))[getValues(m)==1,]
> head(ones)
       x    y
[1,] -42 85.5
[2,] 102 85.5
[3,] 162 85.5
[4,]  42 76.5
[5,] -54 67.5
[6,]  30 67.5
> points(ones[,1],ones[,2])

Paso 1. Genera 1000 pares (xo, yo) centrados en 0 en un cuadro del tamaño de una sola celda. Tenga en cuenta el uso de respara obtener el tamaño de celda:

> pts = data.frame(xo=runif(1000,-.5,.5)*res(m)[1], yo=runif(1000,-.5,.5)*res(m)[2])

Paso 2. Determine en qué celda se encuentra cada uno de los puntos anteriores al muestrear aleatoriamente 1000 valores de 1 a la cantidad de 1 celdas:

> pts$cell = sample(nrow(ones), 1000, replace=TRUE)

Finalmente calcule la coordenada agregando el centro de la celda al desplazamiento. Parcela para verificar:

> pts$x = ones[pts$cell,1]+pts$xo
> pts$y = ones[pts$cell,2]+pts$yo
> plot(m)
> points(pts$x, pts$y)

Aquí hay 10,000 puntos (reemplaza los 1000 anteriores con 10000), trazados con pch=".":

Bastante instantáneo para 10,000 puntos en una trama de 200x300 con la mitad de los puntos como unos. Aumentará en el tiempo linealmente con cuántos en la trama, creo.

Para guardar como un archivo shape, conviértalo en un SpatialPointsobjeto, dele la referencia correcta del sistema de coordenadas (igual que su ráster) y guarde:

> coordinates(pts)=~x+y
> proj4string(pts)=CRS("+init=epsg:4326") # WGS84 lat-long here
> shapefile(pts,"/tmp/pts.shp")

Eso creará un archivo de forma que incluye el número de celda y las compensaciones como atributos.

Siempre que trabajo con grandes conjuntos de datos, me gusta ejecutar herramientas / comandos fuera de QGIS, como desde un script independiente o desde OSGeo4W Shell . No tanto porque QGIS falla (incluso si dice "No responde", probablemente todavía está procesando los datos que puede verificar desde el Administrador de tareas ), sino porque hay más recursos de CPU como RAM disponibles para procesar los datos. QGIS en sí consume una buena cantidad de memoria para ejecutarse.

De todos modos, para ejecutar una herramienta fuera de QGIS ( necesitaría haber instalado QGIS a través del instalador OSGeo4W ), siga los primeros 2 pasos descritos por @gcarrillo en esta publicación: Problema con la importación qgis.core al escribir un script PyQGIS independiente (Sugiero descargar y usar su archivo .bat).

Una vez que se configuran los CAMINOS, escriba pythonen la línea de comando. Para mayor comodidad, copie el siguiente código en un editor de texto como el Bloc de notas, edite los parámetros como el nombre de ruta de su archivo de forma, etc. y luego pegue todo en la línea de comando haciendo clic con el botón derecho> Pegar :

import os, sys
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *

from os.path import expanduser
home = expanduser("~")

QgsApplication( [], False, home + "/AppData/Local/Temp" )

QgsApplication.setPrefixPath("C://OSGeo4W64//apps//qgis", True)
QgsApplication.initQgis()
app = QApplication([])

sys.path.append(home + '/.qgis2/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

shape = home + "/Desktop/Polygon.shp"
result = home + "/Desktop/Point.shp"
general.runalg("qgis:randompointsinlayerbounds", shape, 10000, 0, result)

Usando el script, ejecuté la herramienta Puntos aleatorios en límites de capa para un archivo de forma bastante grande y tardé menos de 20 segundos en generar 10k puntos. Ejecutarlo dentro de QGIS tomó casi 2 minutos, así que al menos para mí, hay una diferencia significativa.

También puede usar GRASS GIS directamente para este trabajo - Muestreo aleatorio estratificado: Muestreo aleatorio de un mapa vectorial con restricciones espaciales :

https://grass.osgeo.org/grass72/manuals/v.random.html#stratified-random-sampling:-random-sampling-from-vector-map-with-spatial-constraints

Además, en el comando se implementa un muestreo aleatorio del mapa vectorial por atributo y algunos otros métodos.

Nota: La versión v.random expuesta en QGIS a través del procesamiento no refleja la funcionalidad completa sino solo una vista simplificada.