Erstellen Sie eine große Anzahl zufälliger Punkte in einem binären Raster?

Ich möchte einen Punktvektordatensatz mit 10000 Punkten (oder mehr) in einem binären Raster erstellen, wobei die Punkte auf Bereiche beschränkt werden sollten, in denen der Rasterwert 1 ist.

Ich habe die folgenden Schritte versucht.

1. Raster polygonisieren
2. QGIS: Vektor -> Recherchetools -> Zufällige Punkte

Dies funktioniert bis zu 2000 Punkten, aber alles, was darüber liegt, führt nur zum Absturz von QGIS.

Gibt es eine Möglichkeit, ein Vektordatensatz mit einer großen Anzahl von Punktmerkmalen zu erstellen, die durch ein binäres Raster (oder die polygonisierte Version davon) eingeschränkt werden?

Die folgenden Tools stehen mir zur Verfügung, von den günstigsten bis zu den ungünstigsten: QGIS, Python, R, ArcGIS

Dies ist, was ich will, nur mit 10x der Punktfunktionen.

Hier ist ein Weg in R:

Machen Sie ein Testraster, 20x30 Zellen, machen Sie 1/10 der Zellen auf 1 gesetzt, zeichnen Sie:

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

Für ein vorhandenes Raster in einer Datei, z. B. ein geoTIFF, können Sie einfach Folgendes tun:

> m = raster("mydata.tif")

Holen Sie sich nun eine Matrix der xy-Koordinaten der 1 Zellen, zeichnen Sie diese Punkte auf und wir sehen, dass wir Zellzentren haben:

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

Schritt 1. Generieren Sie 1000 (xo, yo) Paare, die in einem Feld von der Größe einer einzelnen Zelle auf 0 zentriert sind. Beachten Sie die Verwendung von res, um die Zellengröße zu erhalten:

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

Schritt 2. Ermitteln Sie, in welche Zelle sich jeder der oben genannten Punkte befindet, indem Sie 1000 Werte von 1 bis zur Anzahl von 1 Zellen zufällig auswählen:

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

Berechnen Sie abschließend die Koordinate, indem Sie die Zellmitte zum Versatz hinzufügen. Grundstück zu überprüfen:

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

Hier sind 10.000 Punkte (ersetzen Sie die 1000 oben durch 10000), dargestellt mit pch=".":

Ziemlich augenblicklich für 10.000 Punkte auf einem 200x300- Raster mit der Hälfte der Punkte als Eins. Wird mit der Zeit linear mit der Anzahl der Raster zunehmen, denke ich.

Um als Shapefile zu speichern, konvertieren Sie es in ein SpatialPointsObjekt, geben Sie ihm die richtige Koordinatensystemreferenz (die gleiche wie Ihr Raster) und speichern Sie:

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

Dadurch wird ein Shapefile erstellt, das die Zellennummer und die Offsets als Attribute enthält.

Wenn ich mit großen Datenmengen arbeite, führe ich gerne Tools / Befehle außerhalb von QGIS aus, z. B. über ein eigenständiges Skript oder über die OSGeo4W-Shell . Nicht so sehr, weil QGIS abstürzt (selbst wenn "Nicht reagiert" angezeigt wird, werden wahrscheinlich immer noch die Daten verarbeitet, die Sie im Task-Manager überprüfen können ), sondern weil mehr CPU-Ressourcen wie RAM zur Verarbeitung der Daten verfügbar sind. QGIS selbst verbraucht ziemlich viel Speicher, um ausgeführt zu werden.

Um ein Tool außerhalb von QGIS auszuführen ( Sie müssten QGIS über das OSGeo4W-Installationsprogramm installiert haben ), befolgen Sie die ersten beiden Schritte, wie von @gcarrillo in diesem Beitrag beschrieben: Problem beim Importieren von qgis.core beim Schreiben eines eigenständigen PyQGIS-Skripts (Ich schlage vor, seine .bat-Datei herunterzuladen und zu verwenden).

Sobald die Pfade festgelegt sind, geben Sie pythonin die Befehlszeile ein. Kopieren Sie der Einfachheit halber den folgenden Code in einen Texteditor wie Notepad, bearbeiten Sie die Parameter wie den Pfadnamen Ihres Shapefiles usw. und fügen Sie das Ganze mit der rechten Maustaste> Einfügen in die Befehlszeile ein :

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)

Mit dem Skript habe ich das Werkzeug Zufällige Punkte in Ebenengrenzen für ein ziemlich großes Shapefile ausgeführt und es dauerte weniger als 20 Sekunden, um 10.000 Punkte zu generieren. Das Ausführen in QGIS dauerte fast 2 Minuten. Zumindest für mich gibt es einen signifikanten Unterschied.

Sie können GRASS GIS auch direkt für diesen Job verwenden - Geschichtete Zufallsstichprobe: Zufallsstichprobe aus einer Vektorkarte mit räumlichen Einschränkungen :

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

Zusätzlich werden im Befehl Zufallsstichproben aus der Vektorkarte nach Attributen und einige andere Methoden implementiert.

Hinweis: Die in QGIS durch Verarbeitung bereitgestellte v.random-Version spiegelt nicht die volle Funktionalität wider, sondern nur eine vereinfachte Ansicht.