0

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)
  • I think that the problem could be here: originX ,originY = (ext.xMinimum(),ext.yMinimum()) because raster origin must be: originX ,originY = (ext.xMinimum(),ext.yMaximum()). – xunilk Jul 14 '17 at 22:22
  • Hi @xunilk, Sorry I changed code based on your suggestion but points are same before but moved above the raster. You can show it this link:dropbox.com/s/d1llohdi4qixzqb/… – nickan Jul 15 '17 at 5:47
  • Just now, I remembered your commentary. After seeing the image, I think that your problem it could be here: pixelHeight = ds.rasterUnitsPerPixelY(). It should be: pixelHeight = -ds.rasterUnitsPerPixelY() #minus sign. – xunilk Jul 19 '17 at 2:09
  • Thanks for your attention. Sorry, The point layer result completely moved to below of raster layer for your suggestion. I also tested it for "pixelWidth = -ds.rasterUnitsPerPixelX()" and for it the point layer result moved to left side of raster layer. Did you test the code for your data? I am worry that my question is fake!! – nickan Jul 19 '17 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.