Cómo obtener un vectorial de puntos con valores escogidos de un ráster mediante PyQGIS

Hace unos meses ayudé a un usuario a resolver un problema relacionado con la obtención de un vectorial de puntos para todos los nodata values de un ráster. A continuación presento el código generalizado, en forma de función, para muestrear cualquier tipo de valor y producir el vectorial de puntos correspondiente.

from osgeo import gdal
from PyQt4.QtCore import *
from qgis.core import *
from qgis.gui import QgsMessageBar

def pointsData(value = -999):

    layer = iface.activeLayer()
    
    if layer is None or layer.type() != QgsMapLayer.RasterLayer:
        iface.messageBar().pushMessage("Warning: ",
                                       "There is not an Active Raster Layer",
                                       QgsMessageBar.WARNING, 5)
        return

    #read CRS as EPSG from the raster file:
    myCRS = layer.crs().authid()

    #read Raster path:
    path = layer.dataProvider().dataSourceUri()

    #Create temporary vector layer and add to map
    memory_layer = QgsVectorLayer("Point?crs=" + myCRS, "memory_points", "memory")
    pr = memory_layer.dataProvider()
    QgsMapLayerRegistry.instance().addMapLayer(memory_layer)

    # Add points:

    # Enter editing mode
    memory_layer.startEditing()

    # add fields
    pr.addAttributes( [ QgsField("id", QVariant.Int),
                        QgsField("value", QVariant.Double),
                        QgsField("row",  QVariant.Int),
                        QgsField("column", QVariant.Int) ] )

    raster_or = gdal.Open(path)

    if value == -999:
        value = raster_or.GetRasterBand(1).GetNoDataValue()
    
    (upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = raster_or.GetGeoTransform()

    numpy_array = raster_or.ReadAsArray()

    #columns,rows = numpy_array.shape
    columns = raster_or.GetRasterBand(1).XSize
    rows = raster_or.GetRasterBand(1).YSize

    # add a feature
    fet = QgsFeature()

    i = 0
    
    for row in range(rows):
        for col in range(columns):
            if  (value == numpy_array[row,col]):
                #Add new point:
                x = col * x_size + upper_left_x + (x_size / 2) #add half the cell size
                y = row * y_size + upper_left_y + (y_size / 2) #to centre the point
                fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(x,y)))
                fet.initAttributes(4)
                fet.setAttribute(0, i)
                fet.setAttribute(1, float(numpy_array[row,col]))
                fet.setAttribute(2, row)
                fet.setAttribute(3, col)
                pr.addFeatures( [ fet] )
                
                i += 1

    # Commit changes
    memory_layer.commitChanges()
    
    if i == 0:
        iface.messageBar().pushMessage("Warning: ",
                                       "No Points. Memory layer was not created",
                                       QgsMessageBar.WARNING, 5)
        
        QgsMapLayerRegistry.instance().removeMapLayer(memory_layer.id())
        
        return

Si no se le suministran parámetros a la función se asume el valor por defecto -999 para los nodata values y determina, en consecuencia, el vectorial; una memory layer de puntos. La función maneja las excepciones en el sentido de que si no existe capa activa o ésta no es un ráster se advierte al usuario mediante una QgsMessageBar. Advertencias similares se indican si no existe el valor a muestrear; con lo que el número de puntos es cero y, por tanto, se elimina la memory layer del registro.

Para probar la función al principio utilicé rásteres sintéticos. En el caso siguiente el ráster test.tif tiene todas sus celdas como nodata values y, por tanto, el vectorial de puntos obtenido abarca toda su extensión.

points1

Observen que la tabla atributiva del vectorial captura el valor del ráster para cada punto así como su posición relativa en términos de fila y columna.

En el caso siguiente se utiliza un ráster con valores aleatorios entre 1 y 3. La función se prueba para producir el vectorial de puntos con values iguales a 3.

points2

Probando con rásteres reales (un DEM) para el valor arbitrario de 1572. Se obtienen prácticamente unas “curvas de nivel”.

points3

Para el caso siguiente, recorté el dem un poco para agilizar el proceso de obtención del vectorial. Sé que para la cota de 1368 existe un lago. La imagen a continuación lo confirma:

points4

Un detalle con mayor acercamiento:

points5

El vectorial obtenido tiene 46.532 puntos.

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Cómo obtener un vectorial de puntos con valores escogidos de un ráster mediante PyQGIS

  1. Pingback: Cómo usar una QProgressBar en un Plugin de QGIS | El Blog de José Guerrero

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