5

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
3

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:

enter image description here

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.

enter image description here

  • 1
    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
  • 1
    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
  • 1
    You need to put original question for reopening. I voted for that. – xunilk Jul 10 '17 at 18:33
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).

  • 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

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.