Reading value of point (x, y) in raster and changing its value to another raster using PyQGIS?

Question

I have two raster: raster1 with values [1,2,4,8,16,32,64,128] and raster2, same extension, with all values -9999. Starting from the coordinates of a point, I would like to read its value in raster1 and to do something like this:

If in raster1 the value of the point (x, y) = 1 and the value of the point (x-1, y) = 16, in the raster2 at point (x, y) I leave the -9999 value.

If the value of the point (x, y) = 1 in the raster1 and the value of the point (x-1, y) = 4 then in raster2 I'd write the value 1.

As output I should have the raster2 containing the values 1 only in points that interest me. Can someone help me?

This is what I have done so far:

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

To do that is preferable to use a QgsRasterBlock object to get raster values and python GDAL module to write resulting raster values in a new raster. In this case you only need raster1. Complete code is:

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

Above code was run with a aleatory raster produced with your values [1,2,4,8,16,32,64,128]. Resulting raster was explored for row, column index printed at Python Console of QGIS by using Value Tool plugin. Results obtained were as expected; as it can be observed at next image.

Editing Note:

This is the new script (based in your commentary):

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)

After running it, at the Python Console of QGIS you can observe point coordinates for points in each condition. At Map Canvas, it's visualized each memory point layer separately. Points are at the middle of each cell raster.