Agregar a todas las tramas binarias en un directorio usando calculadora de trama en un Script en Python independiente

Question

Tengo un directorio con el binario de archivos raster (1 banda) que necesito para agregar una secuencia de comandos de python independiente. El problema es que los rásteres tiene distinto número de filas y columnas y resolución diferente (la lectura como arrays de numpy y la adición de ellos no funciona). Cuando yo uso la calculadora ráster, se las arregla para hacer lo que quiero, pero ahora necesito para usarlo fuera de QGIS y hacerlo de forma automática para todos los rásteres en la carpeta. ¿Cómo puedo agregar a todos los de mi rásteres en mi directorio inputFolder y crear un ráster de salida? Tengo el siguiente código:

inputFolder = "/home/usr/Desktop/"

def findRasters(path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join(root, file)

entries = []

for l in findRasters(inputFolder, '*.tif'):
    fileInfo = QFileInfo(l)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(l,baseName)
    QgsMapLayerRegistry.instance().addMapLayer(rlayer) 
    layer = QgsRasterCalculatorEntry()
    layer.ref = l.name() + '@1'
    layer.raster = l 
    layer.bandNumber = 1
    entries.append(layer) 

expression = '(' + entries[1].ref + ' + ' + entries[0].ref + ')'
calc = QgsRasterCalculator(expression,'/home/usr/Desktop/final.tif',
                           'GTiff',
                           layers[0].extent(),
                           layers[0].width(),
                           layers[0].height(),
                           entries)

calc.processCalculation()

No estoy seguro de cómo debo formular mi expression que tiene en cuenta todos los rásteres en la carpeta

Answer 1

Veo más código del que se han agregado desde mi anterior respuesta :)

Corrió esto desde la consola de python, usted puede necesitar otras de las importaciones si llama desde fuera de QGIS.

Hizo unos pequeños cambios, a saber:

• Agregado de las importaciones
• necesidad de establecer la QgsRasterCalculatorEntry.la trama a rlayer (un QgsLayer ejemplo), no el nombre de la capa
• Se ha añadido una lista, layers, que se llena de las capas en el bucle
• Añade una línea para crear el importe total de la expresión.
• Añade el resultado de la capa de

Este se ejecuta hasta el final. Me temo que no tengo muchos buenos conjuntos de datos para probar esto. Pero parece que funciona. Al igual que con el código que se utilice la extensión y la resolución de la primera capa.

import os
import fnmatch
from PyQt4.QtCore import QFileInfo
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

inputFolder = "/path/to/rasters"
outputfilename = '/tmp/final.tif'

def findRasters(path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join(root, file)

entries = []
layers = []

for l in findRasters(inputFolder, '*.tif'):
    print l
    fileInfo = QFileInfo(l)
    baseName = fileInfo.baseName()
    rlayer = QgsRasterLayer(l,baseName)
    QgsMapLayerRegistry.instance().addMapLayer(rlayer) 
    layer = QgsRasterCalculatorEntry()
    layer.ref = rlayer.name() + '@1'
    layer.raster = rlayer
    layers.append(rlayer)
    layer.bandNumber = 1
    print layer
    entries.append(layer) 

reflist = " + ".join([ent.ref for ent in entries]) 
expression = '(' +  reflist + ')'
print expression

calc = QgsRasterCalculator(expression,outputfilename,
                           'GTiff',
                           layers[0].extent(),
                           layers[0].width(),
                           layers[0].height(),
                           entries)

calc.processCalculation()
print "Finished!"

fileInfo = QFileInfo(outputfilename)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(outputfilename, baseName)
if not rlayer.isValid():
    print "Layer failed to load!"
else:
    QgsMapLayerRegistry.instance().addMapLayer(rlayer) 

Answer 2

Usted podría utilizar existentes gdal herramientas.

hay dos gdal utilidades de línea de comandos que se pueda llamar desde la secuencia de comandos mediante os.system o subprocess.call

gdalwarp - utilice esta opción para crear copias de cada trama, la escala para el mismo tamaño (el uso de la -ts opción para establecer la resolución en píxeles, o -tr en unidades de mapa)

gdal_calc - esto le permite hacer la suma en la normalización de los rásteres. Si usted no tiene más de 26 rásteres usted puede hacer esto en una llamada, de lo contrario tendrá que hacer esto en un bucle (agregar raster1+raster2, a continuación, añadir raster3, a continuación, añadir raster4 etc)