Membuat sejumlah besar titik acak dalam raster biner?

9

Saya ingin membuat dataset vektor titik 10.000 poin (atau lebih besar) dalam raster biner, di mana titik harus dibatasi ke area di mana nilai raster adalah 1.

Saya mencoba langkah-langkah berikut.

  1. Polygonize raster
  2. QGIS: Vektor -> Alat Penelitian -> Poin Acak

Ini berfungsi dengan baik hingga 2000 poin tetapi apa pun di atas hanya menyebabkan QGIS lumpuh.

Apakah ada cara untuk membuat dataset vektor dengan sejumlah besar fitur titik yang dikekang oleh biner raster (atau versi poligonnya)?

Alat-alat berikut ini siap membantu saya, peringkat dari yang paling tidak menguntungkan: QGIS, Python, R, ArcGIS

Ini yang saya maksudkan, hanya dengan fitur 10x poin.

1k poin acak

Kersten
sumber
Seberapa besar raster Anda, biasanya?
Spacedman
Yang dalam contoh di atas adalah 19200 x 9600. Raster khas adalah sekitar 10.000 x 10.000 piksel.
Kersten
Oke, semakin banyak RAM yang dimiliki mesin Anda, semakin baik. Saya tidak berani menguji pada raster 10.000x10.000 pada PC kecil saya di sini, meskipun Anda selalu dapat membagi raster, sampel dalam beberapa bagian, dan bergabung ...
Spacedman
mengapa mem-poligon raster? apakah menurut Anda jawaban ini bermanfaat bagi Anda? gis.stackexchange.com/questions/22601/…
Luigi Pirelli
Karena dengan begitu saya dapat menggunakan fungsi "Poin Acak dalam Poligon", sedangkan QGIS tidak memiliki fungsi "Poin Acak di dalam nilai-nilai spesifik fungsi Raster".
Kersten

Jawaban:

7

Berikut cara di R:

Buat raster pengujian, 20x30 sel, buat 1/10 dari sel-sel diatur ke 1, plot:

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

Untuk raster yang ada dalam file, misalnya geoTIFF, Anda bisa melakukan:

> m = raster("mydata.tif")

Sekarang dapatkan matriks koordinat xy dari 1 sel, plot titik-titik itu, dan kami melihat kami memiliki pusat sel:

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

Langkah 1. Hasilkan 1000 (xo, yo) pasangan yang berpusat pada 0 dalam kotak ukuran sel tunggal. Perhatikan penggunaan resuntuk mendapatkan ukuran sel:

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

Langkah 2. Tentukan sel mana dari masing-masing poin di atas yang masuk secara acak dengan mengambil 1000 nilai dari 1 hingga jumlah 1 sel:

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

Akhirnya hitung koordinat dengan menambahkan pusat sel ke offset. Plot untuk memeriksa:

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

Berikut 10.000 poin (ganti 1000 di atas dengan 10.000), diplot dengan pch=".":

poin dalam satu

Cukup instan untuk 10.000 poin pada raster 200x300 dengan setengah poin sebagai raster. Akan bertambah waktu secara linear dengan berapa banyak yang ada di raster, saya pikir.

Untuk menyimpan sebagai shapefile, konversikan ke SpatialPointsobjek, berikan referensi sistem koordinat yang tepat (sama seperti raster Anda) dan simpan:

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

Itu akan membuat shapefile yang menyertakan nomor sel dan offset sebagai atribut.

Spacedman
sumber
Ini terlihat sangat menjanjikan. R saya menjadi sedikit berkarat: Bagaimana saya bisa mengekspor poin ke format vektor (Shapefile, geojson, gml, ... apapun benar-benar) - Saya perlu menyimpan lokasi titik sampel untuk digunakan nanti.
Kersten
Suntingan menunjukkan cara membaca raster dan mengonversi Poin ke
Shapefile
3

Setiap kali saya bekerja dengan dataset besar, saya suka menjalankan alat / perintah di luar QGIS seperti dari skrip mandiri atau dari OSGeo4W Shell . Bukan karena QGIS lumpuh (walaupun dikatakan "Tidak merespons", mungkin masih memproses data yang dapat Anda periksa dari Task Manager ), tetapi karena lebih banyak sumber daya CPU seperti RAM tersedia untuk memproses data. QGIS sendiri menghabiskan banyak memori untuk dijalankan.

Bagaimanapun, untuk menjalankan alat di luar QGIS ( Anda perlu menginstal QGIS melalui installer OSGeo4W ), ikuti 2 langkah pertama seperti yang dijelaskan oleh @ gcarrillo di posting ini: Masalah dengan impor qgis.core saat menulis skrip PyQGIS yang berdiri sendiri (Saya sarankan untuk mengunduh dan menggunakan file .batnya).

Setelah PATHS diatur, ketikkan pythonke dalam baris perintah. Untuk kenyamanan, salin kode berikut ke editor teks seperti Notepad, edit parameter seperti pathname dari shapefile Anda dll. Dan kemudian rekatkan semuanya ke baris perintah dengan Klik kanan> Tempel :

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)

Dengan menggunakan skrip, saya menjalankan titik Acak di alat batas lapisan untuk shapefile yang cukup besar dan butuh waktu di bawah 20 detik untuk menghasilkan 10k poin. Menjalankannya di dalam QGIS membutuhkan waktu hampir 2 menit sehingga setidaknya bagi saya, ada perbedaan yang signifikan.

Yusuf
sumber
1
Alternatif luar biasa, +1. Baru saja menguji ini untuk aplikasi saya dan sementara itu sedikit lebih lambat daripada pendekatan R itu menciptakan hasil yang diinginkan.
Kersten
@Kersten - Luar biasa, senang itu bekerja :)
Joseph
1

Anda juga dapat menggunakan GRASS GIS secara langsung untuk pekerjaan ini - Pengambilan sampel acak bertingkat: Pengambilan sampel acak dari peta vektor dengan batasan spasial :

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

Selain itu, pengambilan sampel acak dari peta vektor berdasarkan atribut dan beberapa metode lain diimplementasikan dalam perintah.

Catatan: Versi v.random yang diekspos di QGIS melalui pemrosesan tidak mencerminkan fungsionalitas penuh tetapi hanya tampilan yang disederhanakan.

markN
sumber