[GIS] Changing pixel values of raster data from QGIS python plugin

pyqgisqgis-pluginsraster

I have been trying my hand at a QGIS plugin. The plug in reads raster layer data into an array, pipes the data to a C++ program, does processing over there, and pipes my results back. After a couple of weeks of fiddling with it, all that works fine.

I would now like to take these results and use them to change the pixel values on the processed raster layers. I have found numerous questions on the site that tell me how to read the raster data from the data provider, but I have not found one that explains to me how I can change the pixel values for the raster layer using Python w/ QGIS.

As a dummy example lets say I have an image of a forest and I want to change all the pixels in a band raster layer with values > 100 to -1 (for example).

How would I do that?

Best Answer

Next function can change raster values greater than 98 for 0.

def changeRasterValues(band):

    fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}

    data_type = band.DataType

    BandType = gdal.GetDataTypeName(band.DataType)

    raster = []

    for y in range(band.YSize):

        scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, data_type)
        values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
        raster.append(values)

    raster = [ list(item) for item in raster ]

    #changing raster values
    for i, item in enumerate(raster):
        for j, value in enumerate(item):
            if value > 98:
                print i, j
                raster[i][j] = 0

    #transforming list in array
    raster = np.asarray(raster)

    return raster

The function is called in the next code and it produces an output raster with the changed pixels. It was tested with one raster tif of 20 x 20 and whose raster values were randomly generated between 1 and 100.

layer = iface.activeLayer()

provider = layer.dataProvider()

path = provider.dataSourceUri()

(raiz, filename) = os.path.split(path)

dataset = gdal.Open(path)

#Get projection
prj = dataset.GetProjection()

#setting band
number_band = 1

band = dataset.GetRasterBand(number_band)

#Get raster metadata
geotransform = dataset.GetGeoTransform()

# Set name of output raster
output_file = "/home/zeito/pyqgis_data/raster_output.tif"

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

raster = changeRasterValues(band)

dst_ds = driver.Create(output_file, 
                       band.XSize, 
                       band.YSize, 
                       number_band, 
                       band.DataType)

#writting output raster
dst_ds.GetRasterBand(number_band).WriteArray( raster )

#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 
srs = osr.SpatialReference(wkt = prj)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset 
dst_ds = None

#Close main raster dataset
dataset = None

After running the code, row and column indexes of changed pixels were printed at the Python Console. When the output raster was loaded to the Map View of QGIS, the pixels changed were compared with the original raster by using the "Value Tool" plugin; as it can be observed at the next image. The changes were as expected.

enter image description here

Related Question