Extraer valores de un ráster bajo una línea (GDAL/OGR/PyQGIS)

Extraer los valores de un ráster bajo puntos es un tema que ya se ha tratado en este Blog. Por otra parte, el muestreo de valores ráster con otro tipo de geometrías (línea, polígono) sugiere, primero, una rasterización del vectorial y luego álgebra de mapas para la incorporación de los values al vectorial rasterizado. Al igual que con el corte (clipper) y la reproyección de rásteres, en QGIS también se usa directamente el comando de GDAL para este procedimiento: gdal_rasterize.

Voy a probar mi método en un script de Python utilizando las herramientas de PyQGIS pensando en un posible plugin. El ráster y el vectorial que me sirvieron como base se encuentran en la imagen siguiente:

rasterize1

La extensión que sirve como base para la rasterización debe ser la del ráster; para que la línea rasterizada quede perfectamente alineada con este. El código propuesto se tiene a continuación:

import os

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

extent = layers[1].extent() #raster layer

xmin = extent.xMinimum()
ymax = extent.yMaximum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()

xsize = layers[1].rasterUnitsPerPixelX()
ysize = layers[1].rasterUnitsPerPixelY()

vlayer = layers[0]  #vector layer
vprovider = vlayer.dataProvider()
uri = vprovider.dataSourceUri()

uri = uri.split("|")

(root, filename) = os.path.split(uri[0])

extension = filename[-4:]

if extension != ".shp":
    
    print "there is not a shapefile"
    
else:
    
    filename = filename[:-4]

    #extent of base raster
    extent = "-te " + str(xmin) + " " + str(ymin) + " " + str(xmax) + " " + str(ymax) + " " 

    #resolution of base raster
    resolution = " -tr " + str(xsize) + " " +  str(ysize) + " "

    attribute_name = " -a id "

    input_vector_line = " -l " + filename + " " + uri[0]
    output_raster = " /home/zeito/pyqgis_data/test_line_raster3.tif"

    cmd = "gdal_rasterize -at " + \
                         extent + \
                     resolution + \
              input_vector_line + \
                 attribute_name + \
                  output_raster

    print cmd

    os.system(cmd)

Cuando se ejecuta en la Python Console de QGIS se obtiene una recta adecuadamente rasterizada; tal como se observa en la imagen siguiente:

rasterize2

Los valores en rojo son cero por lo que el álgebra de mapas con el Raster Calculator de QGIS es un procedimiento bastante sencillo.

Esta entrada fue publicada en GDAL/OGR, PyQGIS, QGIS, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s