Tengo una función que crea aerogeneradores representados como puntos. Esencialmente, usa el código de los puntos aleatorios dentro de la herramienta de polígonos (fijo) aunque con algunos pequeños cambios.
El objetivo es crear puntos aleatorios dentro de los polígonos teniendo en cuenta la distancia mínima especificada. Esto funciona muy bien, especialmente con polígonos que no están cerca de otro (por ejemplo, un solo polígono):
Sin embargo, si el polígono está cerca o junto a otro polígono (por ejemplo, como se muestra a continuación), los puntos de cada polígono pueden estar dentro de la distancia mínima como se muestra en rojo:
¿Cómo podría cambiar el código para que esos puntos en rojo no estén cerca de otro de un polígono cercano?
Idealmente, me gustaría que varios puntos sean reemplazados por un solo punto:
Aquí está el código que se puede reproducir en la Consola Python , se debe seleccionar una capa de polígono con un CRS relevante antes de ejecutar la función:
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:
Disolver y / o convertir los polígonos en partes individuales no parece ayudar mucho ya que los puntos generados todavía parecen estar dentro de la distancia mínima.
Probado en QGIS 2.18.3 .
Respuestas:
Tienes que cambiar dos cosas para que esto funcione. Sin embargo, no obtendrá el máximo de aerogeneradores por área. Para esto, necesitaría ejecutar algunas iteraciones para cada valor y obtener el recuento máximo de puntos.
Me he movido
index = QgsSpatialIndex()ypoints = dict()fuera del bucle for. El código se vería así:EDITAR:
Joseph tenía razón. mis cambios solo funcionaron para un área realmente pequeña. He probado y encontrado una nueva solución moviendo dos variables fuera del ciclo for y cambié la
pointCountvariable.Lo he probado con 500m y este es el resultado (dos intentos diferentes):
fuente
Un método podría ser crear otra función que haga lo siguiente:
Esta es la función que utilicé:
La función se puede ejecutar inmediatamente al final de la
generate_wind_turbines()función utilizando:Esto da los resultados como se muestra en la imagen de la pregunta. Probablemente no sea la más eficiente de las soluciones, pero parece funcionar.
Algunos ejemplos donde los puntos rojos son los generados inicialmente y los puntos mostrados como turbinas con un límite son el resultado final:
generate_wind_turbines(500)generate_wind_turbines(1000)fuente