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
mgri
  • 16,159
  • 6
  • 47
  • 80
HMadadi
  • 1,046
  • 11
  • 24

2 Answers2

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.

YoLecomte
  • 2,895
  • 10
  • 23
  • Hi @YoLecomte, Thanks for answer, but I receive "invalid syntax" in last line "distance=QgsGeometry().fromPoint(vpoint).distance(QgsGeometry().fromPoi‌​nt(rpoint))". – HMadadi Jul 06 '17 at 08: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 06 '17 at 08:27
  • Hi @mgri, I used dist as a variable but It didn't solved. – HMadadi Jul 06 '17 at 08:31
  • 1
    @nikan I've edited my answer, I test the new code in QGIS python console and it works well. – YoLecomte Jul 06 '17 at 08:32
  • @YoLecomte, I think there is a problem in pntRstList because in your edited code I again receive invalid syntax in "rgeometry = QgsGeometry().fromPoi‌​nt(rpoint)" in QGIS python. – HMadadi Jul 06 '17 at 08:38
  • 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 06 '17 at 08:43
  • @nikan What is stored in pntRstList? Can you print it? – YoLecomte Jul 06 '17 at 08:44
  • @mgri, thanks for your comment. Sorry it is not worked for: pntRstList.append(QgsPoint(rspnt)) – HMadadi Jul 06 '17 at 08:49
  • @YoLecomte, yes I can print it. It is a list of pair x, y center of each cell in raster layer. – HMadadi Jul 06 '17 at 08:51
  • 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 06 '17 at 08: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.

mgri
  • 16,159
  • 6
  • 47
  • 80