Cómo modificar un cuadrilátero (trapezoidal o romboidal) para que luzca rectangular mediante geometría analítica

Cuando cursamos cálculo y álgebra lineal algunos cuestionan su utilidad porque la mayoría de los que los dictan los orientan hacia las demostraciones generales excesivamente teóricas. No obstante, cuando trabajamos con SIGs, nos damos cuenta de la importancia de porqué están en el pensum.

Por ejemplo, en este caso, se alega que una digitalización descuidada hizo posible que la mayoría de las edificaciones se observaran con un aspecto distante de su forma verdaderamente rectangular. Tratar de resolver ésto a través de las bounding box se tropieza con el problema de que estas no se orientan con relación al eje longitudinal de los cuadriláteros y, por tanto, el resultado final puede ser exageradamente alto y ancho.

Una forma de resolver esto, sin necesidad de determinar el ángulo de rotación (lo cual luce más complicado), es mediante geometría analítica. La idea es tomar uno de los dos lados del cuadrilátero orientado hacia arriba y obtener la ecuación de la recta paralela que pasa por el “tercer punto” en el sentido horario. Una vez obtenida esta recta sólo queda obtener las ecuaciones de las perpendiculares que pasan por el primer y segundo punto e intersectan la paralela. Esto produce los cuatro puntos necesarios para conformar un rectángulo con las herramientas de QgsGeometry de PyQGIS. El código completo se indica a continuación:

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

n = len(feats)

crs = layer.crs()
epsg = crs.postgisSrid()

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

mem_layer = QgsVectorLayer(uri,
                           'rectangle',
                           'memory')

prov = mem_layer.dataProvider()

for feature in feats:

    geom = feature.geometry()

    xmin, ymin, xmax, ymax = geom.boundingBox().toRectF().getCoords()

    points = feature.geometry().asPolygon()[0]

    for i in range(len(points)-1):
        if points[i][1] == ymax and points[i+1][1] < points[i][1]:
            idx = i
        if points[i][1] == ymax and points[i-1][1] < points[i][1]:
            idx = i-1

    rectangle = []

    #x,y coordinates of first point
    x1 = points[idx][0] 
    y1 = points[idx][1]

    rectangle.append(QgsPoint(x1,y1))

    #x,y coordinates of second point
    x2 = points[idx+1][0] 
    y2 = points[idx+1][1]

    rectangle.append(QgsPoint(x2,y2))

    #slope for first line
    m1 = (y2 - y1) / (x2 - x1)

    #intercept at origin for first line
    int1 = y1 - m1 * x1

    #slope for second line
    m2 = m1

    #x,y coordinates of third point
    x3 = points[idx+2][0] 
    y3 = points[idx+2][1]

    #intercept at origin for second line
    int2 = y3 - m2 * x3

    #first perpendicular
    m3 = -1/m1

    #intercept at origin for second line
    int3 = y2 - m3 * x2

    #intersect point
    x4 = (int3 - int2)/(m2 - m3)
    y4 = m3*x4 + int3

    rectangle.append(QgsPoint(x4, y4))

    #second perpendicular
    m4 = -1/m1

    #intercept at origin for second perpendicular
    int4 = y1 - m4 * x1

    #intersect point
    x5 = (int4 - int2)/(m2 - m4)
    y5 = m4*x5 + int4

    rectangle.extend([QgsPoint(x5, y5),QgsPoint(x1, y1)])

    polygon = []

    polygon.append(rectangle)

    geom = QgsGeometry.fromPolygon(polygon)

    feat = QgsFeature()

    feat.setAttributes([i])
    feat.setGeometry(geom)
    prov.addFeatures( [feat] )

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

El código anterior se probó con la siguiente situación hipotética:

no_rectangular1

y produjo el resultado siguiente después de ser ejecutado en la Python Console de QGIS:

no_rectangular2

Esta entrada fue publicada en PyQGIS, 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