Hasilkan poin dalam banyak fitur dengan mempertimbangkan jarak minimum

8

Saya memiliki fungsi yang menciptakan Turbin Angin yang direpresentasikan sebagai poin. Pada dasarnya, ia menggunakan kode dari titik acak di dalam alat poligon (tetap) meskipun dengan sedikit perubahan.

Tujuannya adalah untuk membuat titik acak di dalam poligon dengan mempertimbangkan jarak minimum yang ditentukan. Ini bekerja sangat baik terutama dengan poligon yang tidak dekat dengan yang lain (misalnya poligon tunggal):

Contoh

Namun, jika poligon dekat atau berdampingan dengan poligon lain (misalnya seperti yang ditunjukkan di bawah), titik-titik dari setiap poligon dapat berada dalam jarak minimum seperti yang ditunjukkan dengan warna merah:

Contoh masalah

Bagaimana saya bisa mengubah kode sehingga titik-titik merah tidak dekat dengan yang lain dari poligon terdekat?

Idealnya, saya ingin beberapa poin diganti oleh satu titik:

Hasil


Berikut adalah kode yang dapat direproduksi dalam Konsol Python , lapisan poligon harus dipilih dengan CRS yang relevan sebelum menjalankan fungsi:

import random
from PyQt4.QtCore import QVariant

def checkMinDistance(point, index, distance, points):
    if distance == 0:
        return True
    neighbors = index.nearestNeighbor(point, 1)
    if len(neighbors) == 0:
        return True
    if neighbors[0] in points:
        np = points[neighbors[0]]
        if np.sqrDist(point) < (distance * distance):
            return False
    return True

def generate_wind_turbines(spacing):
    layer = iface.activeLayer()
    crs = layer.crs()
    # Memory layer
    memory_lyr = QgsVectorLayer("Point?crs=epsg:" + unicode(crs.postgisSrid()) + "&index=yes", "Wind turbines for " + str(layer.name()), "memory")
    QgsMapLayerRegistry.instance().addMapLayer(memory_lyr)
    memory_lyr.startEditing()
    provider = memory_lyr.dataProvider()
    provider.addAttributes([QgsField("ID", QVariant.Int)])
    # Variables
    point_density = 0.0001
    fid = 1
    distance_area = QgsDistanceArea()
    # List of features
    fts = []
    # Create points
    for f in layer.getFeatures():
        fGeom = QgsGeometry(f.geometry())
        bbox = fGeom.boundingBox()
        pointCount = int(round(point_density * distance_area.measure(fGeom)))
        index = QgsSpatialIndex()
        points = dict()
        nPoints = 0
        fid += 1
        nIterations = 0
        maxIterations = pointCount * 200
        random.seed()
        while nIterations < maxIterations and nPoints < pointCount:
            rx = bbox.xMinimum() + bbox.width() * random.random()
            ry = bbox.yMinimum() + bbox.height() * random.random()
            pnt = QgsPoint(rx, ry)
            geom = QgsGeometry.fromPoint(pnt)
            if geom.within(fGeom) and checkMinDistance(pnt, index, spacing, points):
                f = QgsFeature(nPoints)
                f.setAttributes([fid])
                f.setGeometry(geom)
                fts.append(f)
                index.insertFeature(f)
                points[nPoints] = pnt
                nPoints += 1
            nIterations += 1
    provider.addFeatures(fts)
    memory_lyr.updateFields()
    memory_lyr.commitChanges()

generate_wind_turbines(500)

Edit:

Melarutkan dan / atau mengonversi poligon menjadi bagian tunggal tampaknya tidak banyak membantu karena poin yang dihasilkan tampaknya masih berada dalam jarak minimum.

Diuji pada QGIS 2.18.3 .

Yusuf
sumber
1
jika ini pilihan: Sudahkah Anda mencoba menggunakan poligon multipart sebagai poligon input?
LaughU
@ LaughU - Saya baru saja mengujinya dengan poligon multipart tetapi tidak ada poin yang dihasilkan. Ini akan menjadi opsi untuk mengubah fitur singlepart menjadi multi-bagian jika kodenya dapat di-tweak agar berfungsi dengan multi-poligon.
Joseph
Apakah Anda berpikir untuk melarutkan poligon, diikuti oleh multipart-ke-singlepart ke lapisan sementara yang Anda gunakan untuk pembuatan titik? Anda juga dapat menggunakan buffer negatif untuk lapisan sementara untuk menghindari tumpang tindih simbol di bagian tepinya.
Matte
@Matte - Terima kasih, saya mencoba melarutkan semua poligon sebelumnya. Saya mencoba lagi dan mengubahnya menjadi singlepart (tidak yakin apakah ini melakukan apa-apa karena sudah merupakan fitur tunggal) tetapi beberapa titik di dekat tepi berada dalam titik di poligon lain. Saya ingin menghindari penggunaan buffer negatif karena saya ingin mengizinkan poin berada di dekat tepi :)
Joseph

Jawaban:

5

Anda harus mengubah dua hal agar ini berfungsi. Namun, Anda tidak akan mendapatkan turbin angin maksimum per area. Untuk ini, Anda perlu menjalankan beberapa iterasi untuk setiap nilai dan mendapatkan jumlah poin maksimal.

Saya telah memindahkan index = QgsSpatialIndex()dan di points = dict()luar for loop. Kode akan terlihat seperti ini:

import random
def generate_wind_turbines(spacing):
    layer = self.iface.activeLayer()
    crs = layer.crs()
    # Memory layer
    memory_lyr = QgsVectorLayer("Point?crs=epsg:" + unicode(crs.postgisSrid()) + "&index=yes", "Wind turbines for " + str(layer.name()), "memory")
    QgsMapLayerRegistry.instance().addMapLayer(memory_lyr)
    memory_lyr.startEditing()
    provider = memory_lyr.dataProvider()
    provider.addAttributes([QgsField("ID", QVariant.Int)])
    # Variables
    point_density = 0.0001
    fid = 1
    distance_area = QgsDistanceArea()
    # List of features
    fts = []
    # Create points
    points = dict() # changed from me 
    index = QgsSpatialIndex()# changend from me 
    nPoints = 0 # changed in the edit 
    pointCount = 0 # changed in the edit 

    for f in layer.getFeatures():
        fGeom = QgsGeometry(f.geometry())
        bbox = fGeom.boundingBox()
        # changed here as well 
        pointCount = int(round(point_density * distance_area.measure(fGeom))) + int(pointCount)
        fid += 1
        nIterations = 0
        maxIterations = pointCount * 200
        random.seed()
        while nIterations < maxIterations and nPoints < pointCount:
            rx = bbox.xMinimum() + bbox.width() * random.random()
            ry = bbox.yMinimum() + bbox.height() * random.random()
            pnt = QgsPoint(rx, ry)
            geom = QgsGeometry.fromPoint(pnt)
            if geom.within(fGeom) and checkMinDistance(pnt, index, spacing, points):
                f = QgsFeature(nPoints)
                f.setAttributes([fid])
                f.setGeometry(geom)
                fts.append(f)
                index.insertFeature(f)
                points[nPoints] = pnt
                nPoints += 1
            nIterations += 1
    provider.addFeatures(fts)
    memory_lyr.updateFields()
    memory_lyr.commitChanges()

def checkMinDistance( point, index, distance, points):
    if distance == 0:
        return True
    neighbors = index.nearestNeighbor(point, 1)
    if len(neighbors) == 0:
        return True
    if neighbors[0] in points:
        np = points[neighbors[0]]
        if np.sqrDist(point) < (distance * distance):
            return False
    return True

EDIT:

Joseph benar. perubahan saya hanya bekerja untuk area yang benar-benar kecil. Saya telah menguji sekitar dan menemukan solusi baru dengan memindahkan dua variabel dari for for dan mengubah pointCountvariabel.

Saya sudah mengujinya dengan 500m dan ini hasilnya (dua percobaan berbeda):

masukkan deskripsi gambar di sini

Tertawa
sumber
1
Bagus, terima kasih sudah memperhatikan! Saya yakin ini sangat meningkatkan efisiensi =)
Joseph
1
@ Joseph kamu benar. Saya sudah mengujinya dengan area kecil yang berfungsi. Saya telah menambahkan beberapa perubahan kecil dan itu berfungsi untuk saya sekarang. Beri tahu saya jika masih perlu perbaikan :)
LaughU
1
Menarik! Saya akan menguji ini besok dan melaporkan kembali tetapi hasilnya terlihat sangat bagus;)
Joseph
1
Ya, ini bagus! Script Anda juga lebih cepat dan tetap memberikan hasil yang saya cari. Akan memberi Anda hadiah pada akhirnya, harap Anda tidak keberatan, tetapi saya sedikit mengedit kode Anda;)
Joseph
1
@ Joseph Tidak, saya tidak keberatan;) Saya menguji sekitar dengan kode dan lupa untuk membersihkan kekacauan alias spasi putih
LaughU
4

Salah satu metode dapat membuat fungsi lain yang melakukan hal berikut:

  1. Hasilkan buffer dengan jari-jari yang sama dengan jarak di semua titik dan simpan ini dalam lapisan poligon.
  2. Buat daftar yang berisi id semua buffer yang saling bersilangan.
  3. Menghapus buffer menggunakan daftar id ini.
  4. Hasilkan centroid dari buffer yang tersisa dan simpan di layer titik.

Ini adalah fungsi yang saya gunakan:

def wind_turbine_spacing_checker(layer, spacing):
    # Create buffers for points
    poly_layer =  QgsVectorLayer("Polygon?crs=epsg:27700", 'Buffers' , "memory")
    pr = poly_layer.dataProvider() 
    pr.addAttributes([QgsField("ID", QVariant.Int)])
    feat_list = []
    for f in layer.getFeatures():
        poly = QgsFeature()
        f_buffer = f.geometry().buffer((spacing / 2), 99)
        f_poly = poly.setGeometry(QgsGeometry.fromPolygon(f_buffer.asPolygon()))
        poly.setAttributes([1])
        poly.setGeometry(f_buffer)
        feat_list.append(poly)

    pr.addFeatures(feat_list)
    poly_layer.updateExtents()
    poly_layer.updateFields()
    QgsMapLayerRegistry.instance().addMapLayers([poly_layer])

    # Get pairs of intersecting buffer features
    features = [feat for feat in poly_layer.getFeatures()]
    ids = []
    for feat in poly_layer.getFeatures():
        for geat in features:
            if feat.id() != geat.id():
                if geat.geometry().intersects(feat.geometry()):
                    ids.append([feat.id(), geat.id()])

    # Set/sort list and get id of intersecting feature(s)
    for x in ids:
        x.sort()

    ids_sort = set(tuple(x) for x in ids)
    ids_list = [list(x) for x in ids_sort]
    ids_firstItem = [item[0] for item in ids_list]
    final_list = list(set(ids_firstItem))

    # Use ids from final_list to remove intersecting buffer features
    with edit(poly_layer):
        poly_layer.deleteFeatures(final_list)

    # Create new point layer and get centroids from buffers
    # (using final_list to delete the original features may not delete those points where the buffer interesects
    # so best to obtain the centroid of the buffers and output this as a new file)
    result_layer = QgsVectorLayer('Point?crs=epsg:27700&field=id:string', 'Result' , 'memory')
    result_layer.startEditing()
    for feat in poly_layer.getFeatures():
        centroid = feat.geometry().centroid()
        name = feat.attribute("ID")
        centroid_feature = QgsFeature(poly_layer.fields())
        centroid_feature.setGeometry(centroid)
        centroid_feature['ID'] = name
        result_layer.addFeature(centroid_feature)

    result_layer.commitChanges()
    QgsMapLayerRegistry.instance().addMapLayer(result_layer)

Fungsi dapat dieksekusi segera di akhir generate_wind_turbines()fungsi dengan menggunakan:

...
memory_lyr.commitChanges()
wind_turbine_spacing_checker(memory_lyr, spacing)

Ini memberikan hasil seperti yang ditunjukkan pada gambar dalam pertanyaan. Mungkin bukan solusi yang paling efisien tetapi tampaknya berhasil.


Beberapa contoh di mana titik merah adalah yang dihasilkan pada awalnya dan titik yang ditampilkan sebagai turbin dengan batas adalah hasil akhir:

  • generate_wind_turbines(500)

    skenario 1


  • generate_wind_turbines(1000)

    Skenario 2

Yusuf
sumber