El código que se presenta a continuación permite determinar la diferencia en elevación entre un punto y las direcciones cardinales a su alrededor. Inicialmente, emplea un método de numpy para generar puntos correspondientes a un polígono de ocho lados alrededor de cada punto original y luego agruparlos por pares para producir cada dirección cardinal como memory layer.
Posteriomente, cada dirección cardinal es interpolada por el tamaño de la resolución del ráster y cada punto es usado para obtener los raster values a lo largo de las referidas direcciones para finalmente seleccionar la del valor máximo con su índice respectivo.Estos valores son usados para computar las diferencias de elevación entre cada punto original y el punto más elevado sobre cada dirección cardinal y sus respectivas distancias; para luego introducirlos en la correspondiente tabla de atributos del vectorial de línea.
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 92 93 94 95 96 97 98 99 100 101 | import numpy as np print "Processing..." bufferLength = 15000 polygonSides = 8 registry = QgsMapLayerRegistry.instance() layer = registry.mapLayersByName('point') raster = registry.mapLayersByName('utah_demUTM2') xsize = raster[0].rasterUnitsPerPixelX() prov_raster = raster[0].dataProvider() feat_points = [ feat for feat in layer[0].getFeatures() ] points = [ feat.geometry().asPoint() for feat in layer[0].getFeatures() ] epsg = layer[0].crs().postgisSrid() uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes" mem_layer = QgsVectorLayer(uri, 'points', 'memory') prov = mem_layer.dataProvider() group_points = [] for i, point in enumerate(points): int_points = [ QgsPoint(point[0] + np.sin(angle)*bufferLength, point[1] + np.cos(angle)*bufferLength) for angle in np.linspace(0, 2*np.pi, polygonSides, endpoint = False) ] group_points.append(int_points) lines = [] idx_lines = [] for i, group in enumerate(group_points): for point in group: lines.append([points[i], point]) idx_lines.append(i) for group in group_points: feats = [ QgsFeature() for i in range(len(group)) ] for i, feat in enumerate(feats): feat.setAttributes([i]) feat.setGeometry(QgsGeometry.fromPoint(group[i])) prov.addFeatures(feats) QgsMapLayerRegistry.instance().addMapLayer(mem_layer) uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer&field=dif_elev:double&field=distance:double""&index=yes" mem_layer = QgsVectorLayer(uri, 'Output_line', 'memory') prov = mem_layer.dataProvider() feats = [ QgsFeature() for i in range(len(lines)) ] dif_elev = [] distances = [] for i, feat in enumerate(feats): geom = QgsGeometry.fromPolyline(lines[i]) int_points = [] for distance in range(0, int(geom.length()), int(xsize)): values = [] point = geom.interpolate(distance) pt = point.asPoint() int_points.append(pt) for p in int_points: value = prov_raster.identify(p, QgsRaster.IdentifyFormatValue).results()[1] values.append(value) init_value = prov_raster.identify(points[idx_lines[i]], QgsRaster.IdentifyFormatValue).results()[1] elev = np.max(values) - init_value distance = values.index(np.max(values))*xsize dif_elev.append(elev) distances.append(distance) feat.setAttributes([i,float(elev), float(distance)]) feat.setGeometry(geom) prov.addFeatures(feats) QgsMapLayerRegistry.instance().addMapLayer(mem_layer) print "Done!" |
El código anterior fue ejecutado con las capas de punto y ráster de la imagen siguiente.
La tabla de atributos que se observa es la del vectorial de línea correspondiente a las 8 direcciones cardinales para cada punto; donde se han colocado como campos la diferencia de elevación máxima entre ellos y la dirección cardinal a lo largo de un trayecto de 15 km y la separación para la cual ésto ocurre.
No hay comentarios:
Publicar un comentario