¿Cómo sumar valores múltiples para cada píxel en PyQGIS?

Question

Tengo una capa rasterizada y una capa de punto vectorial. Calculé la distancia entre el centro de cada píxel y cada punto basándome en @YoLecomte y @mgri, este enlace y xunilk, este enlace . Para cada pixel hay multidistancia, quiero sumar esta multidistancia para cada pixel en la capa rasterizada y escribirlos en un nuevo archivo rasterizado. En la última parte de este código para sumar las distancias para cada píxel, ¿cuál es el error?

def pixel2coord(x, y):
    xp = (pixelWidth * x) + originX + (pixelWidth/2)
    yp = (pixelHeight * y) + originY + (pixelHeight /2)
    return QgsPoint(xp, yp)

pntLayer = QgsVectorLayer("/Data/points.shp","pointLayer",'ogr')

feats = [ feat for feat in pntLayer.getFeatures() ]

# Open tif file
ds = QgsRasterLayer("/Data/study.tif","Study")

pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()

# extent of the layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()
drive = gdal.GetDriverByName( 'GTiff' )
src_ds = gdal.Open("/Data/study.tif")
outBand = src_ds.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

pntRstList = []

for i in range(0, src_cols):
    for j in range(0, src_rows):
        rspnt = pixel2coord(i,j)
        pntRstList.append(rspnt)

sumP = []
for rpoint in pntRstList:
    for ft in feats:
        vgeometry = ft.geometry()
        rgeometry = QgsGeometry.fromPoint(rpoint)
        dist=vgeometry.distance(rgeometry) # get distance between cell center to point
    if dist < 200:
        sumP.append(dist)
        print sumP

Answer 1

Usé tu código (ligeramente modificado):

from osgeo import gdal
import numpy

def pixel2coord(x, y):
    xp = (pixelWidth * x) + originX + (pixelWidth/2)
    yp = (pixelHeight * y) + originY + (pixelHeight /2)
    return QgsPoint(xp, yp)

#pntLayer = QgsVectorLayer("/Data/points.shp","pointLayer",'ogr')
pntLayer = iface.activeLayer()

feats = [ feat for feat in pntLayer.getFeatures() ]

# Open tif file
ds = QgsRasterLayer("/home/zeito/pyqgis_data/aleatorio1.tif","Study")

pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()

# extent of the layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()
drive = gdal.GetDriverByName( 'GTiff' )
src_ds = gdal.Open("/home/zeito/pyqgis_data/aleatorio1.tif")
outBand = src_ds.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

pntRstList = []

for i in range(0, src_cols):
    for j in range(0, src_rows):
        rspnt = pixel2coord(i,j)
        pntRstList.append(rspnt)

sumP = []
for rpoint in pntRstList:
    for ft in feats:
        vgeometry = ft.geometry()
        rgeometry = QgsGeometry.fromPoint(rpoint)
        dist=vgeometry.distance(rgeometry) # get distance between cell center to point
    if dist < 200:
        tmp_sum += dist
        sumP.append(dist)
        print sumP, numpy.sum(sumP)

para adaptarse a la situación de la siguiente imagen:

Creo que su requerimiento se satisface sólo con una instrucción "numpy.sum(sumP)".

Después de ejecutar el código, las sumas impresas en la Consola de Python para cada lista se ven correctas.

Answer 2

Si quieres sumar, por cada píxel, todas las distancias de las características del punto inferiores a 200 m, puedes usar este código:

for rpoint in pntRstList:
    tmp_sum = 0
    sumP = []
    for ft in feats:
        vgeometry = ft.geometry()
        rgeometry = QgsGeometry.fromPoint(rpoint)
        dist=vgeometry.distance(rgeometry) # get distance between cell center to point
        if dist < 200:
            tmp_sum += dist
            sumP.append(dist)

El tmp_sum es una variable que almacena la suma de todas las distancias inferiores a 200 m, mientras que sumP es una lista que almacena cada distancia por separado (bot tmp_sum y sumP reiniciar para cada nueva iteración).