Friday, 9 October 2015

Measuring distances from each pixel center (in raster) to each point (in vector) in PyQGIS?


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

Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...