viernes, 26 de mayo de 2017

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:


No hay comentarios: