El blog de José Guerrero: Cómo producir rejillas (grid) triangulares mediante PyQGIS 3

En posts pasados (1,2) se consideró la manera de producir rejillas hexagonales mediante PyQGIS 3 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 = QgsPointXY(355472.451507, 4474003.3953491896)
p1 = QgsPointXY(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 = QgsPointXY(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 = QgsPointXY(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.fromPolygonXY([[ QgsPointXY(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.wkbType() == QgsWkbTypes.Polygon:
  
            outFeat.setGeometry(new_geom)
      
            outFeat.setAttributes([k])
     
            prov.addFeatures([outFeat])
 
            k += 1
 
QgsProject.instance().addMapLayer(mem_layer)

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

El blog de José Guerrero

lunes, 3 de agosto de 2020

Cómo producir rejillas (grid) triangulares mediante PyQGIS 3

No hay comentarios: