How to sum multi values for each pixel in PyQGIS?

I have a raster layer and a vector point layer. I calculated distance between center each pixel and each point based on @YoLecomte and @mgri, this link and xunilk, this link. For each pixel there are multi distance,I want to sum this multi distance for each pixel in raster layer and write them into a new raster file. In the last portion of this code to sum distances for each pixel, what is the 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
• I think that your requirement is satisfied only with a 'numpy.sum(sumP)' instruction (please, see my answer). – xunilk Jul 8 '17 at 13:16

I used your code (slightly modified):

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)

for adapting to situation of next image: I think that your requirement is satisfied only with a 'numpy.sum(sumP)' instruction.

After running the code, sums printed at Python Console for each list looks correct. • By the way, why did you retire question about raster from dictionary? It's totally feasible. – xunilk Jul 10 '17 at 16:37
• Sorry @xunilk, I thought it was not a good question, and unfortunately I still could not find the perfect solution for it. So I try it again, thanks. – nickan Jul 10 '17 at 18:20
• I have an answer. – xunilk Jul 10 '17 at 18:21
• I am very interesting to see your answer in gis.stackexchange.com/questions/246996/…. – nickan Jul 10 '17 at 18:28
• You need to put original question for reopening. I voted for that. – xunilk Jul 10 '17 at 18:33

If you want to sum, for each pixel, all the distances from the point features lower than 200 m, you may use this code:

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)

The tmp_sum is a variable that stores the sum of all the distances lower than 200 m, while sumP is a list that stores each separate distance (bot tmp_sum and sumP restart for each new iteration).

• Hi @mgri, Sorry I think I couldn't explain my question clearly. For each pixel there is one or more distance value that sumP in your code can list them, but sum for each pixel is incorrect for tmp_sum. I need result for tmp_sum. – nickan Jul 8 '17 at 12:18
• @nikan I looked the other answer and I confirm I didn't completely understood your issue. However, I'm glad you solved it. – mgri Jul 8 '17 at 16:40