How to sum multi values for each pixel in PyQGIS?

Question

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 08 '17 at 13:16

score 3 · Accepted Answer · answered Jul 08 '17 at 13:13

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. — HMadadi, Jul 10 '17 at 18:20
I am very interesting to see your answer in https://gis.stackexchange.com/questions/246996/how-to-convert-a-dictionary-to-a-raster-layer-in-pyqgis. — HMadadi, 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
You need to do the question visible again and ask to administrator for opening if it was edited adequately. — xunilk, Jul 10 '17 at 18:52

score 2 · Answer 2 · answered Jul 08 '17 at 11:27

2

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).

answered Jul 08 '17 at 11:27

mgri

16,159
6
47
80

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. – HMadadi Jul 08 '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 08 '17 at 16:40

How to sum multi values for each pixel in PyQGIS?

2 Answers2