Qgis: lee el valor de un punto (x, y) en un ráster y cambia su valor a otro ráster

Question

Tengo dos raster: raster1 con los valores de [1,2,4,8,16,32,64,128] y raster2, extensión de la misma, con todos los valores -9999. A partir de las coordenadas de un punto, me gustaría leer su valor en raster1 y hacer algo como esto:

Si en raster1 el valor del punto (x, y) = 1 y el valor del punto (x-1, y) = 16, en el raster2 en el punto (x, y) les dejo el -9999 valor.

Si el valor del punto (x, y) = 1 en el raster1 y el valor del punto (x-1, y) = 4, a continuación, en raster2 me gustaría escribir el valor 1.

Como resultado, tendría yo la raster2 que contiene los valores 1 sólo en los puntos que me interesan. Alguien me puede ayudar?

Esto es lo que he hecho hasta ahora:

import processing
import gdal

layer_1 = iface.addRasterLayer("path/to/raster1.tif")
if layer.isValid():
    QgsMapLayerRegistry.instance().addMapLayer(layer)

processing.runalg('gdalogr:rastercalculator', 'path/to/raster1.tif','1',None,'1',None,'1',None,'1',None,'1',None,'1','A*0-9999',None,5,None,'path/to/raster2.tif')
layer_2 = iface.addRasterLayer('path/to/raster2')
QgsMapLayerRegistry.instance().addMapLayer(layer_2)

x = 603608
y = 4359398
point = QgsPoint(x,y)
ident_1 = layer_1.dataProvider().identify(QgsPoint(x,y),QgsRaster.IdentifyFormatValue)
ident_2 = layer_2.dataProvider().identify(QgsPoint(x,y),QgsRaster.IdentifyFormatValue)
if ident.isValid():
    a = ident_1.results()
    b = ident_2.results()
    if a == 1 &&...

Answer 1

Para hacer eso es preferible usar un QgsRasterBlock objeto de obtener ráster de valores y python GDAL módulo para escribir resultante ráster de valores en una nueva trama. En este caso sólo se necesita raster1. Código completo es:

from osgeo import gdal, osr
import numpy as np

layer = iface.activeLayer()

provider = layer.dataProvider()

extent = provider.extent()

rows = layer.height()
cols = layer.width()

xmin = extent.xMinimum()
ymax = extent.yMaximum()
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

print rows, cols

block = provider.block(1, extent, cols, rows)

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

for i in range(rows):
    for j in range(cols):
        if block.value(i,j) == 1 and block.value(i-1,j) == 16:
            print "yes1", i, j
            block.setValue(i,j,-9999)
            values[i].append(block.value(i,j))
        elif block.value(i,j) == 1 and block.value(i-1,j) == 4:
            print "yes2", i, j
            block.setValue(i,j,1)
            values[i].append(block.value(i,j))
        else:
            values[i].append(block.value(i,j))

raster = np.array(values)

geotransform = [xmin, xsize, 0, ymax, 0, -ysize]

# Create gtif file with rows and columns from parent raster 
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/aleatorio_block.tif"

dst_ds = driver.Create(output_file, 
                       cols, 
                       rows, 
                       1, 
                       gdal.GDT_Int16)

##writting output raster
band = dst_ds.GetRasterBand(1)
band.WriteArray( raster )
band.SetNoDataValue(-9999)

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)

# setting spatial reference of output raster 
epsg = layer.crs().postgisSrid()
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset 
dst_ds = None

Código anterior se ejecutó con un aleatorio de trama producido con sus valores [1,2,4,8,16,32,64,128]. Resultante de trama fue explorado por fila, columna de índice impreso en la Consola de Python de QGIS usando el Valor de la Herramienta de plugin. Los resultados obtenidos fueron como se esperaba; como se puede observar en la siguiente imagen.

La Edición De La Nota:

Este es el nuevo guión (basado en el comentario):

from osgeo import gdal, osr
import numpy as np

layer = iface.activeLayer()

provider = layer.dataProvider()

extent = provider.extent()

rows = layer.height()
cols = layer.width()

xmin = extent.xMinimum()
ymax = extent.yMaximum()
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

print rows, cols

block = provider.block(1, extent, cols, rows)

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

x = xmin + xsize/2
y = ymax - ysize/2

first_cond_points = [] 
second_cond_points = [] 

for i in range(rows):
    for j in range(cols):
        if block.value(i,j) == 1 and block.value(i-1,j) == 16:
            print "yes1", i, j, x, y
            first_cond_points.append(QgsPoint(x,y))
            block.setValue(i,j,-9999)
            values[i].append(block.value(i,j))
        elif block.value(i,j) == 1 and block.value(i-1,j) == 4:
            print "yes2", i, j, x, y
            second_cond_points.append(QgsPoint(x,y))
            block.setValue(i,j,1)
            values[i].append(block.value(i,j))
        else:
            values[i].append(block.value(i,j))

        x += xsize

    y -= ysize

    x = xmin + xsize/2

raster = np.array(values)

geotransform = [xmin, xsize, 0, ymax, 0, -ysize]

# Create gtif file with rows and columns from parent raster 
driver = gdal.GetDriverByName("GTiff")

output_file = "/home/zeito/pyqgis_data/aleatorio_block.tif"

dst_ds = driver.Create(output_file, 
                       cols, 
                       rows, 
                       1, 
                       gdal.GDT_Int16)

##writting output raster
band = dst_ds.GetRasterBand(1)
band.WriteArray( raster )
band.SetNoDataValue(-9999)

#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)

# setting spatial reference of output raster 
epsg = layer.crs().postgisSrid()
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset 
dst_ds = None

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

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

prov = mem_layer.dataProvider()

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

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

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

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

prov = mem_layer.dataProvider()

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

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

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Después de ejecutarlo, en la Consola de Python de QGIS se puede observar el punto de las coordenadas de los puntos en cada condición. En el Mapa de Lona, es visualizada cada punto de memoria capa por separado. Puntos que se encuentran a la mitad de cada celda del raster.