Gere pontos em vários recursos, levando em consideração a distância mínima

8

Eu tenho uma função que cria turbinas eólicas representadas como pontos. Essencialmente, ele usa o código dos pontos Aleatórios dentro da ferramenta polígonos (fixos) , embora com algumas pequenas alterações.

O objetivo é criar pontos aleatórios dentro dos polígonos, levando em consideração a distância mínima especificada. Isso funciona muito bem, especialmente com polígonos que não estão próximos de outro (por exemplo, um único polígono):

Exemplo

No entanto, se o polígono estiver próximo ou adjacente a outro polígono (por exemplo, como mostrado abaixo), os pontos de cada polígono podem estar dentro da distância mínima, como mostrado em vermelho:

Problema de exemplo

Como eu poderia mudar o código para que esses pontos em vermelho não fiquem próximos de outros de um polígono próximo?

Idealmente, gostaria que vários pontos fossem substituídos por um único ponto:

Resultado


Aqui está o código que pode ser reproduzido no Python Console , uma camada de polígono deve ser selecionada com um CRS relevante antes de executar a função:

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)

Editar:

Dissolver e / ou converter os polígonos em peças únicas não parece ajudar muito, pois os pontos gerados ainda parecem estar dentro da distância mínima.

Testado no QGIS 2.18.3 .

Joseph
fonte
1
se for uma opção: você tentou usar um polígono de várias partes como o polígono de entrada?
precisa saber é o seguinte
@ LaughU - Acabei de testar com um polígono de várias partes, mas nenhum ponto é gerado. Seria uma opção para converter os recursos de peça única em multipartes se o código puder ser ajustado para funcionar com polígonos multipartes.
Joseph Joseph
Você pensou em dissolver os polígonos, seguidos de peças múltiplas para peça única em uma camada temporária usada para a geração de pontos? Você também pode usar um buffer negativo para a camada temporária para evitar sobreposições de símbolos nas bordas.
Matte
@ Matte - Obrigado, eu tentei dissolver todos os polígonos antes. Tentei novamente e converti-lo em peça única (não tenho certeza se isso faz alguma coisa, pois já era um recurso único), mas alguns pontos próximos às arestas estão dentro dos pontos dos outros polígonos. Eu gostaria de evitar o uso de buffers negativos, pois quero permitir que os pontos fiquem próximos das bordas :) #
Joseph

Respostas:

5

Você precisa mudar duas coisas para que isso funcione. No entanto, você não obterá o máximo de turbinas eólicas por área. Para isso, você precisaria executar algumas iterações para cada valor e obter a contagem máxima de pontos.

Mudei o index = QgsSpatialIndex()e points = dict()fora do loop for. O código ficaria assim:

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

EDITAR:

Joseph estava certo. minhas alterações funcionaram apenas em uma área muito pequena. Eu testei e encontrei uma nova solução movendo duas variáveis ​​para fora do loop for e alterei a pointCountvariável.

Eu testei com 500m e este é o resultado (duas tentativas diferentes):

insira a descrição da imagem aqui

LaughU
fonte
1
Bom, obrigado por perceber! Eu tenho certeza que isso melhora muito a eficiência =) #
Joseph
1
@ Joseph, você estava certo. Eu testei com pequenas áreas que funcionaram. Eu adicionei algumas pequenas alterações e funciona para mim agora. Deixe-me saber se ele ainda precisa melhorar :)
LaughU
1
Interessante! Eu testarei isso amanhã e apresentarei um relatório, mas o resultado parece realmente bom;) #
Joseph
1
Sim, isso é bom! Seu script também é mais rápido e ainda oferece os resultados que eu estava procurando. Eventualmente, você receberá a recompensa, espero que você não se importe, mas eu editei seu código levemente;) #
Joseph
1
@Joseph Não, eu não me importo;) eu estava testando ao redor com o código e se esqueceu de limpar a bagunça aka whitespaces
LaughU
4

Um método pode ser criar outra função que faça o seguinte:

  1. Gere buffers com o mesmo raio do espaçamento em torno de todos os pontos e os armazene em uma camada de polígono.
  2. Crie uma lista contendo o ID de todos os buffers que se cruzam.
  3. Exclui os buffers usando esta lista de IDs.
  4. Gere centróide dos buffers restantes e armazene-os em uma camada de pontos.

Esta é a função que eu usei:

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)

A função pode ser executada imediatamente no final da generate_wind_turbines()função usando:

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

Isso fornece os resultados, como mostrado na imagem na pergunta. Provavelmente não é a mais eficiente das soluções, mas parece funcionar.


Alguns exemplos em que os pontos vermelhos são os gerados inicialmente e os pontos exibidos como turbinas com um limite são o resultado final:

  • generate_wind_turbines(500)

    Cenário 1


  • generate_wind_turbines(1000)

    Cenário 2

Joseph
fonte