Asignación de valores ráster más cercanos a puntos sobre ‘None’ values con PyQGIS

En el post anterior se consideró la conformación de una memory layer de puntos en las celdas del perímetro de una máscara ráster. La máscara es equivalente a “recortar” un ráster por una geometría tipo polígono y asigna como valores ‘nodata’ a los píxeles que se encuentran más allá del perímetro de la geometría.

La máscara representa un caso complejo a la hora de determinar distancias a puntos más allá de su extensión porque puede asignar un contorno irregular a las celdas de su perímetro si el polígono de “corte” también lo es; en contraste con su aspecto de cuadrilátero antes de la extracción de la máscara.

En el post anterior se incluyó el código para obtener la memory layer de puntos que representan el centro de las celdas del perímetro. Sin embargo, su visualización no es necesaria. Sólo se requieren sus rasgos para determinar las distancias a los puntos más allá de la extensión y escoger aquella que sea mínima. El value a asignar, con el método ‘identify’ de QgsRasterDataProvider, corresponderá entonces a la celda del perímetro ráster más cercana al punto en cuestión.

El código siguiente realiza el proceso y ha sido adaptado para ser ejecutado con procedimientos de PyQGIS.

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()
provider = layers[1].dataProvider()

cols = layers[1].width()
rows = layers[1].height()

extent = layers[1].extent()

print rows, cols

xOrigin = extent.xMinimum()
yOrigin = extent.yMaximum()
pixelWidth = layers[1].rasterUnitsPerPixelX()
pixelHeight = layers[1].rasterUnitsPerPixelY()

x0 = xOrigin + pixelWidth/2
y0 = yOrigin - pixelHeight/2

y = y0

values = [ [] for i in range(rows) ]

for i in range(rows):
    x = x0
    for j in range(cols):
        value = provider.identify(QgsPoint(x, y), QgsRaster.IdentifyFormatValue).results()[1]
        values[i].append(value)
        x += pixelWidth
    y -= pixelHeight

idxs = []

try:
    for i, row in enumerate(values):
        for j, element in enumerate(row):
            if values[i-1][j] == None and element != None or \
               values[i+1][j] == None and element != None or \
               values[i][j-1] == None and element != None or \
               element != None and values[i][j+1] == None:
               idxs.append([i, j])

except IndexError:
    print "Done!"

points = []

for item in idxs:
    x = x0 + item[1]*pixelWidth
    y = y0 - item[0]* pixelHeight
    points.append(QgsPoint(x,y))

crs = layers[1].crs()

epsg = crs.postgisSrid()

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

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

prov = mem_layer.dataProvider()

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

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

prov.addFeatures(feats)

#QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

#for memory layer
feats_m = [ feat for feat in mem_layer.getFeatures()]

#for point layer
feats = [ feat for feat in layers[0].getFeatures()]

for feat in feats:
    value = provider.identify(feat.geometry().asPoint(), 
                              QgsRaster.IdentifyFormatValue).results()[1]
    
    if value != None:
        print value
    
    if value == None:
        distances = []
        for feat_m in feats_m:
            distance = feat.geometry().distance(feat_m.geometry())
            distances.append(distance)

        min_dist = min(distances)
        idx = distances.index(min_dist)
        new_point = points[idx]

        new_value = provider.identify(new_point, 
                                      QgsRaster.IdentifyFormatValue).results()[1]
        print new_value

Cuando se ejecuta el código anterior, considerando la capa Points que contiene 10 features (mezclando la posibilidad de estar sobre values del ráster, nodata values del ráster y más allá de su extensión) de la imagen siguiente:

identify1

se obtienen en la Python Console los resultados a continuación:

62.0
62.0
37.0
54.0
34.0
9.0
37.0
97.0
51.0
22.0

donde todos fueron verificados uno por uno, mediante el plugin ‘Value Tools’, para corroborar que son los esperados.

Esta entrada fue publicada en PyQGIS, QGIS, SIG, 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