1

Relatively new to GIS. I've created a vector layer of points based on row/col position of high population density pixels in a raster layer. The layers do not align even though they have the same CRS, presumably because they do not have the same extent.

Is it possible to edit the extent of the raster to match that of the vector layer or when creating the new vector layer can I define extent based on my raster layer?

enter image description here

  • If my answer solved your issue you should mark it as accepted. – xunilk Apr 28 '17 at 21:17
2

To do that you need to get simultaneously x,y coordinates of each point by using some raster parameters (as xmin, ymax, and x, y resolution). Next code uses these parameters for obtaining a point memory layer where each point is situated at the middle of each raster cell.

layer = iface.activeLayer()

provider = layer.dataProvider()

extent = layer.extent()
xmin,ymin,xmax,ymax = extent.toRectF().getCoords()
rows = layer.height()
cols = layer.width()
xSize = layer.rasterUnitsPerPixelX()
ySize = layer.rasterUnitsPerPixelY()

block = provider.block(1, extent, cols, rows)

xinit = xmin + xSize/2
yinit = ymax - ySize/2

x = xinit
y = yinit

points = []
values = []

for i in range(rows):
    for j in range(cols):
        value = block.value(i,j)
        if value > 150 and value < 255:
            points.append(QgsPoint(x,y))
            values.append(value)
        x += xSize
    y -= ySize
    x = xinit

epsg = layer.crs().postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=value:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i, values[i]])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

I used this population raster; where the criterion to get points was they represent values greater than 150 (254 = 10,000 persons/km2 and 255 = No data). After running the code at Python Console of QGIS I got:

enter image description here

Points over raster match as expected.

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.