domingo, 12 de noviembre de 2017

Generando buffer "cartesiano" con PyQGIS

La circunferencia es el lugar geométrico de los puntos de un plano que equidistan de otro punto fijo llamado centro. Su ecuación general, en coordenadas cartesianas, puede ser encontrada aquí.

A partir de la referida ecuación puede despejarse "y" como función de "x" para luego generar diferentes valores de "x", igualmente espaciados, y determinar cuales serían sus correspondientes valores de "y".

En la generación automática llevada a cabo de esta manera se producen puntos repetidos que deben ser luego eliminados antes de visualizarlos. Otro aspecto que también se debe cuidar es el de generarlos en orden secuencial a fin de que también se puedan obtener las correspondientes figuras cerradas (línea, polígono).

El código siguiente genera 80 puntos, igualmente espaciados, para construir una circunferencia y el círculo correspondiente alrededor del centro (392127.908689, 4449521.71257) y con radio de 100 metros. Los puntos, la circunferencia y el círculo se despliegan todos como memory layers en la Map View de QGIS.

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
from math import sqrt

pt = (392127.908689, 4449521.71257)
r = 100

x0 = pt[0]
y0 = pt[1]

xf = pt[0] + r

f0 = r/20

xnew = x0 - r

my_x = [ xnew ]

while(xf - xnew > 0):
    xnew += f0
    my_x.append(xnew)

points1 = [ QgsPoint(x, sqrt(r**2 - (x - x0)**2) + y0) for x in my_x ]

points2 = [ QgsPoint(x, -sqrt(r**2 - (x - x0)**2) + y0) for x in my_x ]

points2.reverse()

points = points1 + points2

epsg = 32612

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

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

prov = mem_layer.dataProvider()

new_points = [ point for i, point in enumerate(points) if point not in points[:i] ]

feats = [ QgsFeature() for i in range(len(new_points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(new_points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

new_points.append(new_points[0])

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

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

prov = mem_layer.dataProvider()

line = [ new_points ]

feats = [ QgsFeature() for i in range(len(line)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(line[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

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

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

prov = mem_layer.dataProvider()

polygon = [ [ new_points ] ]

feats = [ QgsFeature() for i in range(len(polygon)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolygon(polygon[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Después de ejecutado el código anterior, para mejor visualización, en la imagen siguiente sólo se presentan los puntos y la circunferencia. En los extremos derecho e izquierdo se observa la estructura discreta que caracteriza a los puntos usados en su generación. Esto se debe a los elevados cambios de pendiente que se tienen en esas regiones para pequeños cambios en los valores de "x".


El mismo efecto puede también observarse en la imagen a continuación para el círculo generado.


Obviamente, los buffer circulares no se generan en los SIGs con las ecuaciones cartesianas sino con las formas paramétricas de la circunferencia. No obstante, es un buen ejercicio para poner de manifiesto estas limitaciones y para generar puntos de la circunferencia igualmente espaciados en el eje x en el orden secuencial correcto.

No hay comentarios: