Créer une grande quantité de points aléatoires dans un raster binaire?

Je veux créer un jeu de données vectorielles ponctuelles de 10000 points (ou plus) dans un raster binaire, où les points doivent être contraints aux zones où la valeur du raster est 1.

J'ai essayé les étapes suivantes.

1. Polygoniser un raster
2. QGIS: Vector -> Outils de recherche -> Points aléatoires

Cela fonctionne très bien jusqu'à 2000 points, mais tout ce qui précède provoque simplement le blocage de QGIS.

Existe-t-il un moyen de créer un jeu de données vectorielles avec un grand nombre d'entités ponctuelles contraintes par un raster binaire (ou sa version polygonisée)?

Les outils suivants sont à ma disposition, classés du plus favorable au moins favorable: QGIS, Python, R, ArcGIS

C'est ce que je veux, seulement avec 10x les fonctionnalités ponctuelles.

Voici un moyen dans R:

Faire un raster de test, 20x30 cellules, faire 1/10 des cellules définies sur 1, tracer:

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

Pour un raster existant dans un fichier, par exemple un geoTIFF, vous pouvez simplement faire:

> m = raster("mydata.tif")

Maintenant, obtenez une matrice des coordonnées xy des 1 cellules, tracez ces points, et nous voyons que nous avons des centres de cellules:

> 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])

Étape 1. Générez 1000 paires (xo, yo) centrées sur 0 dans une boîte de la taille d'une seule cellule. Notez l'utilisation de respour obtenir la taille des cellules:

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

Étape 2. Déterminez dans quelle cellule chacun des points ci-dessus va entrer en échantillonnant au hasard 1000 valeurs de 1 au nombre de 1 cellules:

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

Enfin, calculez les coordonnées en ajoutant le centre de la cellule au décalage. Terrain à vérifier:

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

Voici 10000 points (remplacez les 1000 ci-dessus par 10000), tracé avec pch=".":

Quasiment instantané pour 10 000 points sur un raster 200x300 avec la moitié des points comme des uns. Augmentera le temps de façon linéaire avec le nombre de pixels dans le raster, je pense.

Pour enregistrer en tant que fichier de formes, convertissez-le en SpatialPointsobjet, donnez-lui la bonne référence de système de coordonnées (la même que votre raster) et enregistrez:

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

Cela créera un fichier de formes qui comprend le numéro de cellule et les décalages en tant qu'attributs.

Chaque fois que je travaille avec de grands ensembles de données, j'aime exécuter des outils / commandes en dehors de QGIS, comme à partir d'un script autonome ou d' OSGeo4W Shell . Pas tant parce que QGIS plante (même s'il dit "Ne répond pas", il traite probablement toujours les données que vous pouvez vérifier à partir du Gestionnaire des tâches ), mais parce que davantage de ressources CPU telles que la RAM sont disponibles pour traiter les données. QGIS lui-même consomme une bonne partie de la mémoire pour fonctionner.

Quoi qu'il en soit, pour exécuter un outil en dehors de QGIS ( vous devez avoir installé QGIS via le programme d'installation OSGeo4W ), suivez les 2 premières étapes décrites par @gcarrillo dans cet article: Problème avec l'importation qgis.core lors de l'écriture d'un script PyQGIS autonome (Je suggère de télécharger et d'utiliser son fichier .bat).

Une fois les CHEMINS définis, tapez pythondans la ligne de commande. Pour plus de commodité, copiez le code suivant dans un éditeur de texte tel que le Bloc-notes, modifiez les paramètres tels que le chemin d'accès de votre fichier de formes, etc., puis collez le tout dans la ligne de commande en cliquant avec le bouton droit sur> Coller :

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)

À l'aide du script, j'ai exécuté l' outil Points aléatoires dans les limites des couches pour un fichier de formes assez volumineux et il a fallu moins de 20 secondes pour générer 10 000 points. L'exécuter dans QGIS a pris près de 2 minutes, donc au moins pour moi, il y a une différence significative.

Vous pouvez également utiliser GRASS GIS directement pour ce travail - Échantillonnage aléatoire stratifié: Échantillonnage aléatoire à partir d'une carte vectorielle avec des contraintes spatiales :

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

De plus, l'échantillonnage aléatoire de la carte vectorielle par attribut et quelques autres méthodes sont implémentés dans la commande.

Remarque: La version v.random exposée dans QGIS via le traitement ne reflète pas la fonctionnalité complète mais juste une vue simplifiée.