3

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?

  • Good question. Related post, using osgeo: How to Create a raster layer from a numpy array using pyqgis? – ArMoraer May 4 '16 at 14:33
  • @ArMoraer yeah I saw that, but what I really want is to change the raster values of a currently available layer instead of creating a new one if I can help it. It is my backup option though – Semicolons and Duct Tape May 4 '16 at 14:35
  • Yes, I realize that. Creating a new raster is already a rather tedious task with PyQGIS, let alone updating one! Worst-case scenario: you could consider deleting you raster before re-creating it... Sorry I can't help more. – ArMoraer May 4 '16 at 14:42
  • I think that it is preferable to produce a new raster with the changed pixels for preserving original raster. You can see my answer. – xunilk May 27 '16 at 22:55
3

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.