Tworzysz dużą liczbę losowych punktów w binarnym rastrze?

Chcę utworzyć zestaw danych wektora punktowego o wartości 10000 punktów (lub więcej) w binarnym rastrze, w którym punkty powinny być ograniczone do obszarów, w których wartość rastra wynosi 1.

Próbowałem następujące kroki.

1. Poligonizuj raster
2. QGIS: Wektor -> Narzędzia badawcze -> Losowe punkty

Działa to dobrze do 2000 punktów, ale wszystko powyżej powoduje tylko awarię QGIS.

Czy istnieje sposób na utworzenie zbioru danych wektorowych z dużą liczbą cech punktowych ograniczonych przez binarny raster (lub jego poligonizowaną wersję)?

Do mojej dyspozycji są następujące narzędzia, od najbardziej do najmniej korzystnych: QGIS, Python, R, ArcGIS

To jest to, o co zamierzam, tylko z 10-krotnymi cechami punktowymi.

Oto sposób na R:

Zrób testowy raster, komórki 20x30, zrób 1/10 komórek ustawionych na 1, wykreśl:

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

W przypadku istniejącego rastra w pliku, na przykład geoTIFF, możesz po prostu:

> m = raster("mydata.tif")

Teraz uzyskaj macierz współrzędnych xy 1 komórek, wykreśl te punkty i zobaczymy, że mamy centra komórkowe:

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

Krok 1. Wygeneruj 1000 par (xo, yo), które są wyśrodkowane na 0 w polu wielkości pojedynczej komórki. Zwróć uwagę na użycie, resaby uzyskać rozmiar komórki:

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

Krok 2. Sprawdź, do której komórki wchodzi każdy z powyższych punktów, losowo próbkując 1000 wartości od 1 do liczby 1 komórek:

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

Na koniec oblicz współrzędną, dodając środek komórki do odsunięcia. Działka do sprawdzenia:

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

Oto 10 000 punktów (zamień 1000 powyżej na 10000), z wykreśleniem pch=".":

Niemal natychmiastowe dla 10 000 punktów na rastrze 200 x 300 z połową punktów jako jedynymi. Wydaje mi się, że wraz z upływem czasu będzie rósł liniowo.

Aby zapisać jako plik kształtu, przekonwertuj na SpatialPointsobiekt, nadaj mu właściwe odniesienie do układu współrzędnych (taki sam jak twój raster) i zapisz:

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

Spowoduje to utworzenie pliku kształtu zawierającego numer komórki i przesunięcia jako atrybuty.

Ilekroć pracuję z dużymi zestawami danych, lubię uruchamiać narzędzia / polecenia poza QGIS, na przykład z samodzielnego skryptu lub z powłoki OSGeo4W . Nie tyle dlatego, że QGIS ulega awarii (nawet jeśli mówi „Nie odpowiada”, prawdopodobnie nadal przetwarza dane, które można sprawdzić w Menedżerze zadań ), ale ponieważ więcej zasobów procesora, takich jak RAM, jest dostępnych do przetwarzania danych. Sam QGIS zużywa sporą część pamięci.

W każdym razie, aby uruchomić narzędzie poza QGIS ( musiałbyś zainstalować QGIS za pomocą instalatora OSGeo4W ), wykonaj pierwsze 2 kroki opisane w @gcarrillo w tym poście: Problem z importem qgis.core podczas pisania samodzielnego skryptu PyQGIS (Sugeruję pobranie i użycie jego pliku .bat).

Po ustawieniu PATHS wpisz pythonw wierszu polecenia. Dla wygody skopiuj poniższy kod do edytora tekstów, takiego jak Notatnik, edytuj parametry, takie jak nazwa ścieżki pliku shapefile itp., A następnie wklej całość do wiersza poleceń, klikając prawym przyciskiem myszy> Wklej :

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)

Za pomocą skryptu uruchomiłem narzędzie Losowe punkty w obwiedniach warstw dla dość dużego pliku kształtu i wygenerowanie 10 000 punktów zajęło mniej niż 20 sekund. Uruchomienie go w QGIS zajęło prawie 2 minuty, więc dla mnie jest znacząca różnica.

Możesz także użyć GRASS GIS bezpośrednio do tego zadania - Stratified random lospling: Losowe próbkowanie z mapy wektorowej z ograniczeniami przestrzennymi :

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

Ponadto w poleceniu zaimplementowano losowe próbkowanie z mapy wektorów według atrybutu i kilka innych metod.

Uwaga: Wersja v.losowa ujawniona w QGIS podczas przetwarzania nie odzwierciedla pełnej funkcjonalności, a jedynie uproszczony widok.