sábado, 11 de noviembre de 2017

Cómo recuperar puntos rasterizados con el orden como value y construir una bounding box orientada

En la imagen a continuación se tienen unos puntos rasterizados con una resolución 10x10 donde como value se ha incorporado el orden de los puntos. El objetivo es vectorizar nuevamente esos puntos, en el orden correcto, para contruir con ellos una bounding box orientada.

Para ello se va a emplear una mezcla de procedimientos de GDAL/OGR y PyQGIS.


El código completo desarrollado es el siguiente:

 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
from osgeo import gdal, osr
import processing

raster_path = '/home/zeito/pyqgis_data/trajectory_L_shape.tif'

dataset = gdal.Open(raster_path)
band = dataset.GetRasterBand(1)
geotransform = dataset.GetGeoTransform()

data = band.ReadAsArray()

nonzero = zip(data.nonzero()[0], data.nonzero()[1])

values = [ data[element[0], element[1]] for element in nonzero ]

points_idx = [ [QgsPoint(geotransform[0] + geotransform[1]/2 + indices[1]*geotransform[1],
               geotransform[3] + geotransform[5]/2 + indices[0]*geotransform[5]), values[i]]
               for i, indices in enumerate(nonzero) ]

points_idx.sort(key=lambda x: x[1]) #to sort by second column

wkt = dataset.GetProjection()

srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
epsg = srs.GetAttrValue("AUTHORITY", 1)

uri = "Point?crs=epsg:" + epsg + "&field=id:integer&field=order:integer""&index=yes"

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

prov = mem_layer.dataProvider()

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

for i, feat in enumerate(feats):
    feat.setAttributes([i,int(points_idx[i][1])])
    feat.setGeometry(QgsGeometry.fromPoint(points_idx[i][0]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

processing.runandload('qgis:orientedminimumboundingbox',
                      mem_layer,
                      False,
                      None)

Cuando se ejecuta en la Python Console de QGIS se obtienen los resultados que se encuentran a continuación:


En la tabla de atributos del vectorial de puntos el orden se incrementa con el id; tal como se espera.

No hay comentarios: