2

I have one raster layer and point vector layer. I want to measure distance between center of each pixel in raster to each point in vector layer. I extracted x,y center of each cell and added them to a list, but I can't make a relation between this list and vector points to calculate distances. How can I measure distances?

from qgis.core import *
from osgeo import gdal, ogr, osr

pntLayer = QgsVectorLayer("/Data/Villages.shp","pointLayer")

feat = [ 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()

outDs = drive.Create("/Data/pop_Rst.tif", src_cols, src_rows, 1, gdal.GDT_Float32)
outBand = outDs.GetRasterBand(1)

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

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

pntRstList = []

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

print pntRstList 
4

According to this question, How to calculate distance in meter between QgsPoint objects in PyQgis?, it's easy.

You just need to return a QgsPoint with your function pixel2coord() :

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

Then you will store a list of QgsPoint than you can compare to your vector Point (store in your feat list) like that:

for vpoint in feat :
    for rpoint in pntRstList:
        vgeometry = QgsGeometry().fromPoint(vpoint) 
        rgeometry = QgsGeometry().fromPoi‌​nt(rpoint)
        dist=vgeometry.distance(rgeometry)

I'm not sure it's the most efficient way to do this but it should works.

  • Hi @YoLecomte, Thanks for answer, but I receive "invalid syntax" in last line "distance=QgsGeometry().fromPoint(vpoint).distance(QgsGeometry().fromPoi‌​nt(rpoint))". – nickan Jul 6 '17 at 8:20
  • 1
    @nikan I would suggest not using a method as a variable name because may confuses QGIS. Try using dist as a variable, for example. – mgri Jul 6 '17 at 8:27
  • 1
    @nikan I've edited my answer, I test the new code in QGIS python console and it works well. – YoLecomte Jul 6 '17 at 8:32
  • 1
    @nikan my comment was not to be intended as a solution, but as a suggestion. However, does it change something if you use pntRstList.append(QgsPoint(rspnt)) instead of what you provided? – mgri Jul 6 '17 at 8:43
  • 1
    @nikan ok that's probably why it doesn't works... You need a list of QgsPoint. I have edited my answer to store the point in a variable and return this variable. It works for me in python console. – YoLecomte Jul 6 '17 at 8:59
3

There were some missing pieces of code. I used the solution by @YoLecomte for obtaining a working solution:

from qgis.core import *
from osgeo import gdal, ogr, osr
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/Villages.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' )
outDs = drive.Create("/Data/pop_Rst.tif", src_cols, src_rows, 1, gdal.GDT_Float32) # ?
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)

#print pntRstList

for ft in feats:
    vgeometry = ft.geometry()
    for rpoint in pntRstList:
        rgeometry = QgsGeometry.fromPoint(rpoint)
        dist=vgeometry.distance(rgeometry)
        #print dist

If it works also for you, please accept the solution by @YoLecomte since I simply merged his code with yours.

  • 1
    It's work, That's great job. – nickan Jul 6 '17 at 9:48
  • 1
    Great! I understand now... I miss to get the geometry of the vector point before calculating the distance... – YoLecomte Jul 6 '17 at 10:05

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.