3

Given a series of rasters, and the raster image of the minimum values of the series, I am looking for a python script capable of selecting, for each raster, the pixels whose values are greater than the minimum, i.e. "A/(A>minimum)", being A a multiple raster.

I made the following script, it seems like is calculating something, but at the end I get the message: 'NoneType' object is not iterable. Any idea for solution?

##NAME=name
##Raster=group
##A=multiple raster
##minimum=raster
##X=output raster

import qgis, os, processing 
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator


# Split rasters
layers = A.split(';')
output = X + "/"
suffix = "_suffix.tif"
# Get layer object
layer1 = processing.getObject(minimum)

def x(A,minimum,output):
    for ras in layers:
        fileInfo = QFileInfo(ras)
        baseName = fileInfo.baseName()
        lyr = QgsRasterLayer(ras, baseName)

        # Get layer object
        lyr = processing.getObjectFromUri(ras)
        entries = []
        ras = QgsRasterCalculatorEntry()
        ras.ref = 'lyr@1'
        ras.raster = lyr
        ras.bandNumber = 1
        entries.append(ras)  

        #define raster 1 ("minimum")
        raster1=QgsRasterCalculatorEntry()
        raster1.ref='minimum@1'
        raster1.raster=minimum
        raster1.bandNumber=1
        entries.append(raster1)

        calc = QgsRasterCalculator( 'lyr@1/(lyr@1>minimum@1)', output + lyr.name() + suffix, 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
        calc.processCalculation()

for results in x(layers, layer1,output + "*" + suffix): 
    fileInfo = QFileInfo(results)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

In the case of solution proposed by Ali, I don't understand if the solution is suitable only for a fixed number of input rasters (7). Actually I need a flexible tool which work independently of the number of input, being them 3 or 10 for instance. Anyway I was not able to make a useful script due to my basic knowledge about python.

On the other hand I also tried the solution as proposed by xunilk but I get the following messages:

xmin = layers[0].extent().xMinimum() IndexError: list index out of range

In both cases I would need to use the script within the "processing toolbox" as a standalone tool as in the image:

enter image description here

Actually in the case of using just two layers I have a solution which is the following:

##NAME=name
##A=raster
##Minimum=raster
##X=output raster

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

# Get layer object
layer1 = processing.getObject(A)
# Get layer object
layer2 = processing.getObject(Minimum)

def x(A,Minimum,output):
    entries=[]

    #define raster 1 ("A")
    raster1=QgsRasterCalculatorEntry()
    raster1.ref='A@1'
    raster1.raster=A
    raster1.bandNumber=1
    entries.append(raster1)

    #define raster 2 ("Minimum")
    raster2=QgsRasterCalculatorEntry()
    raster2.ref='Minimum@1'
    raster2.raster=Minimum
    raster2.bandNumber=1
    entries.append(raster2)

    #X Processing
    calc=QgsRasterCalculator('A@1/(A@1>Minimum@1)',output,'GTiff',Minimum.extent(),Minimum.width(),Minimum.height(),entries)
    calc.processCalculation()

x(layer1,layer2,X)

The problem here is that I need the factor A as a multiple raster instead of a single raster in order to iterate the same operation, but I don't know how to solve it.

Is there the possibility of using the script for making a tool from the "processing toolbox" by entering the variables directly from the QGIS workspace by using the following way:

##NAME=name
##Raster=group
##A=multiple raster
##minimum=raster
##X=output raster

Instead of referring the script to folder path?

  • This: output + "*" + suffix is, probably, part of your issue. Caracter "*" is in the name of output raster. By the way, its easier to solve this by using the 'identify' method of QgsRasterDataProvider. – xunilk Jun 12 '16 at 16:04
  • @Nico_77 I have updated my answer to provide an alternative solution – Ali Jun 18 '16 at 12:30
3

It looks like your issue here is with your function x. Since nothing is being returned from the function, you cannot iterate over it.

I would suggest making use of the yield keyword. This will return a generator which you can iterate over. Essentially it returns each result in the function in turn. For example:

def a_function(x):
    for i in range(x):
        yield i

You can then use your for results... loop to loop through each yielded result. Based on the function above:

for result in a_function(7):
    # Optionally do something with the result
    print result

This will simply print each value in range(7) in turn, i.e.,

0
1
2
3
4
5
6

If I understand your script correctly, I think you should yield your new raster path in your x function. This will then be passed into your final loop as results and therefore your QFileInfo call will point to your new raster :

def x(A,minimum,output):
    for ras in layers:
        fileInfo = QFileInfo(ras)
        baseName = fileInfo.baseName()
        lyr = QgsRasterLayer(ras, baseName)

        # Get layer object
        lyr = processing.getObjectFromUri(ras)
        entries = []
        ras = QgsRasterCalculatorEntry()
        ras.ref = 'lyr@1'
        ras.raster = lyr
        ras.bandNumber = 1
        entries.append(ras)  

        #define raster 1 ("minimum")
        raster1=QgsRasterCalculatorEntry()
        raster1.ref='minimum@1'
        raster1.raster=minimum
        raster1.bandNumber=1
        entries.append(raster1)

        calc = QgsRasterCalculator( 'lyr@1/(lyr@1>minimum@1)', output + lyr.name() + suffix, 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
        calc.processCalculation()

        # Now yield your output raster
        yield output + lyr.name() + suffix

UPDATE

My answer above does indeed work for any number of rasters, however, having spent a bit more time with this, I have noticed some more fundamental issues with your initial script e.g. unused variables, redundant objects. I have therefore written my own interpretation of your script.

The following is my Python script which I successfully used to remove the black borders from Landsat imagery and have them loaded into the QGIS map canvas. As you mentioned limited python knowledge I have tried to provide useful comments throughout the script.

Essentially the process_raster function will, given an input raster, minimum raster, and output path, carry out the processing you require and return the path to the new file.

This script will loop through a list of rasters (this list can be any length!) and, for each one, create the new raster and store its path. It then loads the new raster into the map canvas before continuing to the next raster in the list.

import qgis, os, processing 
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator

# Define list of input raster paths
rasters_list = [r"PATH_TO_IN_RASTER1.tif", r"PATH_TO_IN_RASTER2.tif", r"ETC..."]

# Define path to minimum raster
# Do not use minimum or min as variable names as they are reserved words.
min_rast = r"PATH_TO_MIN_RASTER.tif"

# Define output (braces allow insertion of variable into string)
output = r"X:\{}_suffix.tif"

def process_raster(in_rast,min_rast,out_rast):
    '''
    in_rast(string): Path to input raster
    min_rast(string): Path to minimum raster
    out_rast(string): Path to which output raster will be saved
    '''

    entries = []

    # Get layer object from input raster
    in_lyr = processing.getObjectFromUri(in_rast)
    # Create first raster calculator entry
    rast_ent1 = QgsRasterCalculatorEntry()
    rast_ent1.ref = 'lyr@1'
    rast_ent1.raster = in_lyr
    rast_ent1.bandNumber = 1
    entries.append(rast_ent1)  

    # Get layer object from minimum raster
    min_lyr = processing.getObjectFromUri(min_rast)
    # Create second raster calculator entry
    rast_ent2=QgsRasterCalculatorEntry()
    rast_ent2.ref='minimum@1'
    rast_ent2.raster=min_lyr
    rast_ent2.bandNumber=1
    entries.append(rast_ent2)

    out_rast = out_rast.format(in_lyr.name()) # layer name will be inserted into the output string, replacing the braces
    calc = QgsRasterCalculator( 'lyr@1/(lyr@1>minimum@1)', out_rast, 'GTiff', in_lyr.extent(), in_lyr.width(), in_lyr.height(), entries )
    calc.processCalculation()

    return out_rast

for raster in rasters_list:
    new_rast = process_raster(raster,min_rast,output)
    # since the process_raster function returns the output path to the
    # new tiff, the new_rast variable now stores this path and can be used
    # to add the raster to the map canvas.
    fileInfo = QFileInfo(new_rast)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

UPDATE 2

In order to get the script functioning as a processing tool I defined the parameters and parameter types at the top of the script and removed the hard coded file paths.

##Rasters Minimum=name
##rasters_list=multiple raster
##min_rast= raster
##output=output raster

import qgis, os, processing 
from PyQt4.QtCore import QFileInfo
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator

# Define list of input raster paths
# (as input is now a semi-colon separated string it must be split)
rasters_list = rasters_list.split(";")

def process_raster(in_rast,min_rast,out_rast):
    '''
    in_rast(string): Path to input raster
    min_rast(string): Path to minimum raster
    out_rast(string): Path to which output raster will be saved
    '''

    entries = []

    # Get layer object from input raster
    in_lyr = processing.getObjectFromUri(in_rast)
    # Create first raster calculator entry
    rast_ent1 = QgsRasterCalculatorEntry()
    rast_ent1.ref = 'lyr@1'
    rast_ent1.raster = in_lyr
    rast_ent1.bandNumber = 1
    entries.append(rast_ent1)  

    # Get layer object from minimum raster
    min_lyr = processing.getObjectFromUri(min_rast)
    # Create second raster calculator entry
    rast_ent2=QgsRasterCalculatorEntry()
    rast_ent2.ref='minimum@1'
    rast_ent2.raster=min_lyr
    rast_ent2.bandNumber=1
    entries.append(rast_ent2)

    # The specified output raster is split into its name and extension components
    # then reconstructed to include the name of the raster currently being processed.
    # Therefore each output raster path is autogenerated from the single input parameter.
    out_rast = os.path.splitext(out_rast)[0] + "_" + in_lyr.name() + os.path.splitext(out_rast)[1]
    calc = QgsRasterCalculator( 'lyr@1/(lyr@1>minimum@1)', out_rast, 'GTiff', in_lyr.extent(), in_lyr.width(), in_lyr.height(), entries )
    calc.processCalculation()

    return out_rast

for raster in rasters_list:
    new_rast = process_raster(raster,min_rast,output)
    # since the process_raster function returns the output path to the
    # new tiff, the new_rast variable now stores this path and can be used
    # to add the raster to the map canvas.
    fileInfo = QFileInfo(new_rast)
    path = fileInfo.filePath()
    baseName = fileInfo.baseName()
    layer = QgsRasterLayer(path, baseName)
    QgsMapLayerRegistry.instance().addMapLayer(layer)

Note: This requires only a single tif specified as the output; the path of each individual output raster will be derived from this and therefore the specified output wont actually exist. Therefore, be sure to uncheck Open output file after running algorithm. The process of adding the outputs to the map canvas is taken care of in the script.

  • Is there the possibility of using the script for making a tool from the "processing toolbox" by entering the variables directly from the QGIS workspace by using the following way: ##NAME=name ##Raster=group ##A=multiple raster ##minimum=raster ##X=output raster Instead of referring the script to folder path? – Nico_77 Jun 18 '16 at 20:00
  • @Nico_77 I've added a further update which I've successfully used as a tool in QGIS. – Ali Jun 18 '16 at 21:33
  • It was exactly what I was looking for! it works to me as well, perfectly. Thanks a lot! – Nico_77 Jun 18 '16 at 22:28
2

Its easier to solve this issue by using the 'identify' method of QgsRasterDataProvider and a little bit of GDAL Python. I used a very simple rasters (20 x 20) to test my approach. The first one was an aleatory raster with values between 1 and 100. The another one was the 'minimum' raster; whose values were all 50. It looks like this:

enter image description here

The code is:

from osgeo import gdal, osr
import numpy as np

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()

xmin = layers[0].extent().xMinimum()
ymax = layers[0].extent().yMaximum()

xsize = layers[0].rasterUnitsPerPixelX()
ysize = layers[0].rasterUnitsPerPixelY()

rows = layers[0].height()
columns = layers[0].width()

minimums = [ [] for i in range(rows) ] 
points = [ [] for i in range(rows) ]

y = ymax - ysize/2

#conforming a QgsPoints matrix for all rasters
for i in range(rows):

    x = xmin + xsize/2

    for j in range(columns):
        points[i].append(QgsPoint(x,y))
        x += xsize
    y -= ysize

provider_lay = layers[0].dataProvider()
provider_min = layers[1].dataProvider()

values = [ [] for i in range(rows) ]

for i in range(rows):
    for j in range(columns):
        value1 = provider_lay.identify(points[i][j], QgsRaster.IdentifyFormatValue).results()[1]
        value2 = provider_min.identify(points[i][j], QgsRaster.IdentifyFormatValue).results()[1]
        if value1 > value2:
            values[i].append(value1)
        else:
            values[i].append(-1)

#transforming matrix in array
raster = np.asarray(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/values_gtmin.tif"

dst_ds = driver.Create(output_file, 
                       columns, 
                       rows, 
                       1, 
                       gdal.GDT_Float32)

#writting output raster
dst_ds.GetRasterBand(1).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)

crs = layers[0].crs()

epsg = crs.postgisSrid()

# setting spatial reference of output raster 
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dst_ds.SetProjection( srs.ExportToWkt() )

#Close output raster dataset 
dst_ds = None

After running the above code at the Python Console of QGIS, I got a raster (values_gtmin.tif) where pixels greater than minimum image were adequately selected and pixels lesser (or equal) than minimum image were replaced by -1. At the next image, it can be observed the corroboration of obtained results with the help of 'Value Tool' plugin. Black areas are pixels with -1 values.

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.