Distancia entre geometrías tipo polígono utilizando PyQGIS

Continuando con el uso de los métodos y situaciones que pueden plantearse con la clase QgsGeometry hoy se va a considerar la distancia entre geometrías de tipo polígono. Con vectoriales tipo punto es fácil intuir que la distancia entre ellos es la recta que los une. Con polígonos cabría esperar que esta fuese la menor distancia de separación.

Para probar este aserto vamos a considerar la situación correspondiente a la imagen siguiente:

distance1

Consideramos un caso sencillo previendo futuras automatizaciones. Por ello, intuitivamente se observa que la menor distancia debería ser el segmento de recta entre el vértice del polígono azul que apunta hacia el verde y que es perpendicular al lado superior (mas cercano) de este último.

El código que permite determinar rápidamente esa distancia es el siguiente:

layers = mapcanvas.layers()

feat_p2 = layers[0].getFeatures().next()

feat_p9 = layers[1].getFeatures().next()

geom2 = feat_p2.geometry()
geom9 = feat_p9.geometry()

distance1 = geom2.distance(geom9)

print distance1

cuyo valor viene siendo 18406.4938087 m.

La siguiente porción de código permite conocer los índices de los puntos que están más cercanos a cada simetría.

polygon2 = geom2.asPolygon()
polygon9 = geom9.asPolygon()

idx_p9 = [ geom9.closestVertexWithContext(point)[1]
           for point in polygon2[0] ]

d = {}

d.fromkeys(idx_p9).keys()

nl1 = [d.setdefault(x,x) for x in idx_p9 if x not in d ]

print nl1

points_p9 = [ polygon9[0][i] for i in nl1 ]

print points_p9

idx_p2 = [ geom2.closestVertexWithContext(point)[1]
           for point in polygon9[0] ]

nl2 = [d.setdefault(x,x) for x in idx_p2 if x not in d ]

print nl2

points_p2 = [ polygon2[0][i] for i in nl2 ]

print points_p2

El resultado de la ejecución es la siguiente:

[5]
[(428798,4.4507e+06)]
[2, 3]
[(409514,4.45616e+06), (411269,4.44282e+06)]

El índice 5 corresponde al vértice más cercano del polígono azul. Para encontrar la distancia más corta hasta el polígono verde solo hay que usar el método ‘closestSegmentWithContext’ para determinar cual es el punto de intersección con ese segmento más cercano. Una vez determinado, se imprime la distancia y se crea la respectiva memory layer.

La porción faltante del código es la siguiente:

p3 = geom2.closestSegmentWithContext(points_p9[0])[1]

d3 = sqrt(points_p9[0].sqrDist(p3))

print d3

crs = layers[0].crs()

epsg = crs.postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&field=distance&index=yes"

dist = QgsVectorLayer(uri, 
                      'dist', 
                      'memory')

QgsMapLayerRegistry.instance().addMapLayer(dist)

prov = dist.dataProvider()

feat = QgsFeature()

feat.setGeometry(QgsGeometry.fromPolyline([points_p9[0], p3]))
feat.setAttributes([1, feat.geometry().length()])

prov.addFeatures([feat])

y permite verificar que el segmento de recta que une ambas geometrías en la imagen siguiente:

distance2

es el esperado y con la misma longitud impresa al comienzo: 18406.4938087 m.

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s