domingo, 2 de agosto de 2020

Cortando una rejilla (grid) hexagonal creada con PyQGIS 3

En el post anterior se consideró la creación de una rejilla hexagonal. Esta contenía unos “salientes” laterales y en las partes superior e inferior de la rejilla. Además, el punto de referencia usado para su confección luce desplazado por la necesidad de encajar las filas de la manera adecuada.

Las modificaciones al código, además de aquellas que requieren necesariamente las nuevas clases implementadas en QGIS 3, se necesita ubicar el punto de inicio y final del corte para crear una geometría asociada a la clase QgsRectangle. Esta geometría se intersectará con la correspondiente a cada feature en el mismo loop de creación de la memory layer.

El código completo se incluye a continuación:

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import numpy as np
 
bufferLength = 1000
polygonSides = 6
 
inc_x = 2 * bufferLength * np.cos(np.pi*30/180)
inc_y = inc_x * np.cos(np.pi*30/180)
 
p_init = QgsPointXY(354972.451507, 4473926.04508)
 
p1 = QgsPointXY(p_init.x(), p_init.y() - bufferLength/2)
 
points = []
 
x = p1.x() + inc_x/2
y = p1.y()
 
rows = 8
cols = 10
 
p_final = QgsPointXY(p_init.x() + cols*inc_x - inc_x/2, p_init.y() - (rows-1)*inc_y - bufferLength)
     
for i in range(rows):
    for j in range(cols):
        point = QgsPointXY(x, y)
        points.append(point)
        x += inc_x
    y -= inc_y
    ver = i%2
 
    if ver ==0:
        h = 0
    else:
        h = inc_x/2
    x = p1.x() + h
 
epsg = 32612
 
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
 
mem_layer = QgsVectorLayer(uri,
                           'hexgrid',
                           'memory')
 
prov = mem_layer.dataProvider()
 
extent = QgsRectangle(p_init, p_final)
 
geom_rect = QgsGeometry.fromRect(extent)
 
for i, point in enumerate(points):
    outFeat = QgsFeature()
    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) ]])
 
    new_geom = geom.intersection(geom_rect)
 
    outFeat.setGeometry(new_geom)
     
    outFeat.setAttributes([i])
    prov.addFeatures([outFeat])
 
QgsProject.instance().addMapLayer(mem_layer)

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


En la imagen anterior, sólo como referencia, se incluyen los puntos inicial y final de corte. Por otra parte, la lista ‘points’ contiene todos los centroides de los hexágonos antes del corte. Si se desea, se puede conformar una memory layer de puntos con ellos.

No hay comentarios: