Leer el valor de un ráster y cambiarlo según condiciones dadas mediante PyQGIS

El código que se presenta a continuación lee todos los valores de un ráster mediante un objeto de la clase QgsRasterBlock y los cambia según condiciones establecidas para imprimirlos en un nuevo ráster. Para facilitar la ubicación de tales valores, complementariamente, se crean memory layers de puntos, según cada condición, y se colocan en el medio de las celdas del ráster.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116

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)

El código anterior se ejecutó con un ráster aleatorio cuyos valores son [1,2,4,8,16,32,64,128] y las dos condiciones establecidas se cumplieron sólo para 8 celdas; generando diferencialmente dos tipos de puntos (6 para la primera condición y 2 para la segunda). Esto se ve reflejado en la imagen siguiente: