2

I would like to create point layer for storing values retrieved from a raster file.

My idea was to iterate through all raster cells. If cell value is different than NULL, I wish to create a new point which represent mentioned raster cell.

Unfortunately, result point layer, do not match raster data - points are far away from the raster file. After long time of investigation, I have still no idea what might be wrong. Could you help?

from osgeo import gdal
from PyQt4.QtCore import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsVectorLayer, QgsField, QgsMapLayerRegistry, QgsFeature, QgsGeometry, QgsPoint

#read CRS as EPSG from the raster file:
myCRS = MyRasters['first'].crs().authid()
#read Raster path:
RasterfilePath=MyRasters['first'].dataProvider().dataSourceUri()
#Create temporary vector layer and add to map
vl = QgsVectorLayer("Point?crs=" + myCRS, "temporary_points", "memory")
pr = vl.dataProvider()
QgsMapLayerRegistry.instance().addMapLayer(vl)
# Add points:

pr = vl.dataProvider()

# Enter editing mode
vl.startEditing()

# add fields
pr.addAttributes( [ QgsField("first", QVariant.Double),
    QgsField("second",  QVariant.Int),
    QgsField("third", QVariant.Double) ] )

raster_or = gdal.Open(RasterfilePath)

nodata = raster_or.GetRasterBand(1).GetNoDataValue()
(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = raster_or.GetGeoTransform()

numpy_array = raster_or.ReadAsArray()

width,height = numpy_array.shape
#width=raster_or.GetRasterBand(1).XSize
#height=raster_or.GetRasterBand(1).YSize

for row in range(0,width):
    for col in range(0,height):
        if  ~(nodata==numpy_array[row,col]):
            #Add new point:
            x = row * x_size + upper_left_x + (x_size / 2) #add half the cell size
            y = col * y_size + upper_left_y + (y_size / 2) #to centre the point
            # add a feature
            fet = QgsFeature()
            fet.setGeometry( QgsGeometry.fromPoint(QgsPoint(x,y )
            fet.initAttributes(3)
            fet.setAttribute(0, float(numpy_array[row,col]))
            pr.addFeatures( [ fet ] )
# Commit changes
vl.commitChanges()
  • Width and height are inverted (and there are other errors produced by bad copy/paste). – xunilk Apr 10 '15 at 9:34
2

They are inverted (width, height). Try this:

.
.
.

for row in range(height):
    for col in range(width):
.
.
.     

Editing Note:

I edited your code:

from osgeo import gdal
from PyQt4.QtCore import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsVectorLayer, QgsField, QgsMapLayerRegistry, QgsFeature, QgsGeometry, QgsPoint

#read CRS as EPSG from the raster file:
myCRS = iface.activeLayer().crs().authid()
#read Raster path:
RasterfilePath=iface.activeLayer().dataProvider().dataSourceUri()
#Create temporary vector layer and add to map
vl = QgsVectorLayer("Point?crs=" + myCRS, "temporary_points", "memory")
pr = vl.dataProvider()
QgsMapLayerRegistry.instance().addMapLayer(vl)
# Add points:

pr = vl.dataProvider()

# Enter editing mode
vl.startEditing()

# add fields
pr.addAttributes( [ QgsField("first", QVariant.Double),
    QgsField("second",  QVariant.Int),
    QgsField("third", QVariant.Double) ] )

raster_or = gdal.Open(RasterfilePath)

nodata = raster_or.GetRasterBand(1).GetNoDataValue()
(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = raster_or.GetGeoTransform()

numpy_array = raster_or.ReadAsArray()

#width,height = numpy_array.shape
width=raster_or.GetRasterBand(1).XSize
height=raster_or.GetRasterBand(1).YSize

# add a feature
fet = QgsFeature()

for row in range(height):
    for col in range(width):
        if  (nodata == numpy_array[row,col]):
#            #Add new point:
            x = col * x_size + upper_left_x + (x_size / 2) #add half the cell size
            y = row * y_size + upper_left_y + (y_size / 2) #to centre the point
            fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(x,y)))
            fet.initAttributes(3)
            fet.setAttribute(0, float(numpy_array[row,col]))
            fet.setAttribute(1, row)
            fet.setAttribute(2, col)
            pr.addFeatures( [ fet] )


# Commit changes
vl.commitChanges()

and I tested it by using a raster with nodata values:

enter image description here

It works! Points are in nodata values places.

  • Very nice! Always awesome to see people testing and reporting back, already gave you +1 =) – Joseph Apr 10 '15 at 11:12
  • Thanks. The original code had some errors but now it works nicely. – xunilk Apr 10 '15 at 11:39

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.