I have a point vector layer that have a field data. I want to disaggregate this field data to a raster layer by distance. It means when we away from a point, value in each pixel have smaller value and pixels closer to each point have larger values (such as IDW). I also have a raster layer that pixel values = 1 shows study area and pixel values = 0 shows out of study area. Calculation is performed in study area. After computations I must get a point layer (memory layer) that each point is pixel center from raster layer of study area. After run this code a point layer (memory layer) created, that points arranged in it rotated 180 degree related to raster layer.It must not rotated. What is error(s)?
from qgis.core import *
from osgeo import gdal
import numpy
from collections import defaultdict
import itertools
def pixel2coord(x, y):
xp = (pixelWidth * x) + originX + (pixelWidth/2)
yp = (pixelHeight * y) + originY + (pixelHeight /2)
return QgsPoint(xp, yp)
pntLayer = iface.activeLayer()
feats = [ feat for feat in pntLayer.getFeatures() ]
# Open tif file
ds = QgsRasterLayer("/Data/Sample_Polygon.tif","Sample_Polygon")
pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()
# extent of the Raster 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/Sample_Polygon.tif")
outBand = src_ds.GetRasterBand(1)
outData = numpy.zeros((src_cols, src_rows), numpy.float32)
pntRstList = []
provider = ds.dataProvider()
block = provider.block(1, ext, src_cols, src_rows)
for i in range(0, src_cols):
for j in range(0, src_rows):
if block.value(i,j) == 1:
rspnt = pixel2coord(i,j)
pntRstList.append(rspnt)
idxMeanDist = pntLayer.fieldNameIndex("MeanDist")
idxPopulaton = pntLayer.fieldNameIndex("Population")
# dictionary for store distance from each pixel to points in its search radius
dictDist = {}
dictDist00 = {}
# search radius
dictPop = {}
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 < 400:
meanDist = ft.attributes()[idxMeanDist] # get mean dist values in point layer
weightPop = abs((meanDist**2 - dist**2)/(meanDist**2 + dist**2))
pntPop = ft.attributes()[idxPopulaton]
keyR = rpoint
dictDist.setdefault(keyR, [])
dictDist[keyR].append(weightPop)
dictPop.setdefault(keyR, [])
dictPop[keyR].append(pntPop)
for rpoint in pntRstList:
if rpoint in dictDist:
pass
else:
dictDist00.setdefault(rpoint, 0)
dictDistOne = {k: v for k, v in dictDist.iteritems() if len(v) == 1}
dictDistMore = {k: v for k, v in dictDist.iteritems() if len(v) > 1}
newDictMore = {}
for key, value in (dictDistMore.items()):
newDictMore.setdefault(key, [])
for item in value:
newDictMore[key].append(item/(sum(value)))
newDict = newDictMore.copy()
newDict.update(dictDistOne)
C = {k: [x * y for x, y in zip(v, dictPop[k])] for k, v in newDict.items() if k in dictPop}
C_sum = {}
for key, value in (C.items()):
C_sum.setdefault(key, [])
C_sum[key].append(sum(value))
mergeList = [ [ point, mergeDict[point] ] for point in mergeDict ]
points = [ QgsPoint(item[0][0], item[0][1]) for item in mergeList ]
values = [ item[1] for item in mergeList ]
new_values = []
for item in values:
if type(item) == list:
new_values.append(item[0])
else:
new_values.append(item)
epsg = 32639
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=value:real""&index=yes"
mem_layer = QgsVectorLayer(uri,
'points',
'memory')
prov = mem_layer.dataProvider()
feats = [ QgsFeature() for i in range(len(points)) ]
for i, feat in enumerate(feats):
feat.setAttributes([i, new_values[i]])
feat.setGeometry(QgsGeometry.fromPoint(points[i]))
prov.addFeatures(feats)
QgsMapLayerRegistry.instance().addMapLayer(mem_layer)
originX ,originY = (ext.xMinimum(),ext.yMinimum())
because raster origin must be:originX ,originY = (ext.xMinimum(),ext.yMaximum())
.pixelHeight = ds.rasterUnitsPerPixelY()
. It should be:pixelHeight = -ds.rasterUnitsPerPixelY() #minus sign
.