Cómo producir rejillas (grid) triangulares mediante PyQGIS

En posts pasados (1,2) se consideró la manera de producir rejillas hexagonales mediante PyQGIS y en este le toca el turno a las triangulares. Aunque pueda parecer extraño, a menor número de lados mayor complejidad.

En las hexagonales, la “distancia buffer” era igual al lado del polígono; pero esto no se cumple en las triangulares. Es lo primero que hay que resolver. Por otra parte, en la secuencia en una fila, la alternacia de triangulos exige que unos esten rotados 60 grados para que encajen; por lo que hay que implementar ésto en el código. Con relación al ajuste entre las diferentes filas, éste puede hacerse por desplazamiento. Sin embargo, esto produce una zona algo compleja en la cual las intersecciones con el área de corte produce puntos en lugar de polígonos. Esto hay que preverlo en el código para evitar errores al seleccionar los features.

El código completo se presenta a continuación para producir un rejilla triangular de 8 filas por 20 columnas que comienza en el punto (355472.451507, 4474003.3953491896) y con un EPSG:32612. El punto final de corte se determina automáticamente con base al número de filas y columnas.

import numpy as np

sideLength = 1000 
bufferLength = sideLength/(2*np.sin(np.pi*60/180))
polygonSides = 3

p_init = QgsPoint(355472.451507, 4474003.3953491896)
p1 = QgsPoint(p_init.x()-sideLength/2, p_init.y()-bufferLength)

points = []

inc_x = sideLength/2
inc_y = (sideLength/2)*np.tan(np.pi*30/180)

x = p1.x() + inc_x
y = p1.y()

rows = 8
cols = 20

if cols%2 == 0:
    coef = 1
else:
    coef = 1.5
    cols += 1

p_final = QgsPoint(p_init.x() + (cols+2)*sideLength/2-coef*sideLength, p_init.y() - rows*sideLength*np.cos(np.pi*30/180) )

for i in range(rows):
    ver2 = i%2
    for j in range(cols+2):
        ver1 = j%2
        point = QgsPoint(x, y)
        points.append(point)
        x += inc_x
        if ver1 == 0:
            y += inc_y
        else:
            y -= inc_y
    
    y -= sideLength*np.cos(np.pi*30/180)
    
    if ver2 == 0:
        h = 0
    else:
        h = inc_x

    x = p1.x() + h

epsg = 32612
 
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
 
mem_layer = QgsVectorLayer(uri,
                           'grid',
                           'memory')
 
prov = mem_layer.dataProvider()

extent = QgsRectangle(p_init, p_final)

geom_rect = QgsGeometry.fromRect(extent)

k=0

for i, point in enumerate(points):
    ver = i%2
    
    geom = QgsGeometry.fromPolygon([[ QgsPoint(point[0] + np.sin(angle)*bufferLength, point[1] + np.cos(angle)*bufferLength)
                                   for angle in np.linspace(0, 2*np.pi, polygonSides + 1, endpoint = True) ]])

    if ver != 0:
        geom.rotate(60,point)

    if geom.intersects(geom_rect):
        outFeat = QgsFeature()
        new_geom = geom.intersection(geom_rect)
        
        if new_geom.type() == QGis.Polygon:
 
            outFeat.setGeometry(new_geom)
     
            outFeat.setAttributes([k])
    
            prov.addFeatures([outFeat])

            k += 1

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Después de ejecutado en la Python Console de QGIS, se obtiene:

triangular1

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