5

I am a new in PyQGIS, I want to create a Point layer from a raster layer. In table of point layer created, there is 3 fields by name of "X_coord", "Y_coord" and "Azimuth", each point in Point layer shows x ("X_coord" field), y ("Y_coord") center coordinate of each pixel from raster layer and "Azimuth" field shows azimuth from another single point layer to each center pixels. When I run this script for raster layer with 67*52 pixels in 30 meter resolution time for running is 5 seconds, for 133*105 pixels time is 15 seconds and for 291*220 pixels time is 22 minutes! I think this is unusual.How can I decrease time for run this script for large data?

My script is:

from qgis.core import *
from osgeo import gdal
from PyQt4.QtCore import *
import time

startTime = time.clock()

pnt = QgsVectorLayer(r'E:/Sample_Dataset_Lorestan/Point.shp', 'POINT','ogr')
for feat in pnt.getFeatures():
   geom = feat.geometry()
   pntxy=geom.asPoint()

# Open tif file
ds = gdal.Open(r'E:\Sample_Dataset_Lorestan\subset_a1.tif')

numpy_array = ds.ReadAsArray()
# GDAL affine transform parameters

geotransform = ds.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
rotationX = geotransform[2]
rotationY = geotransform[4]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]

cols=ds.RasterXSize
rows=ds.RasterYSize

def pixel2coord(x, y):
   xp =  (pixelWidth * x) + originX + (pixelWidth/2)
   yp =  (pixelHeight * y) + originY + (pixelHeight /2)
   return(xp, yp)

#Create temporary vector layer and add to map
vl = QgsVectorLayer("Point?crs=", "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("X_coord", QVariant.Int), QgsField("Y_coord", QVariant.Int), QgsField("Azimuth", QVariant.Int) ] )

# add a feature
fet = QgsFeature()

# get columns and rows of your image from gdalinfo

for row in range(0,rows):
   for col in range(0,cols):
    rspnt = pixel2coord(col,row)
    rspnt2 = QgsPoint(rspnt[0], rspnt[1])
    fet.setGeometry(QgsGeometry.fromPoint(rspnt2))
    fet.initAttributes(3)
    fet.setAttribute(0, rspnt[0] )
    fet.setAttribute(1, rspnt[1])

    az = pntxy.azimuth(rspnt2)
    if az < 0.0:
         azim = 360 + az
    else:
        azim = az

    fet.setAttribute(2, azim)
    pr.addFeatures( [fet] )

# Commit changes
vl.commitChanges()
print "completed within:", time.clock() - startTime, "seconds", "\n"
1

Try this version (modified from gene's response):

vl = QgsVectorLayer("Point?crs=", "temporary_points2", "memory")
vl.dataProvider().addAttributes([ QgsField("X_coord", QVariant.Int), QgsField("Y_coord", QVariant.Int), QgsField("Azimuth", QVariant.Int) ])

new_features = []
for row in range(0,rows):
   for col in range(0,cols):
      rspnt = pixel2coord(col,row)
      fet = QgsFeature()
      fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(rspnt[0], rspnt[1])))
      fet.initAttributes(3)
      fet.setAttribute(0, rspnt[0] )
      fet.setAttribute(1, rspnt[1])
      az = pntxy.azimuth(QgsPoint(rspnt[0], rspnt[1]))
      if az < 0.0:
         azim = 360 + az
      else:
         azim = az
      fet.setAttribute(2, azim)
      new_features.append(fet)

vl.dataProvider().addFeatures( new_features )
vl.updateExtents()
vl.updateFields()    
QgsMapLayerRegistry.instance().addMapLayers([vl])

This version only calls addFeatures once, rather than once per pixel in the loop.

Update: Here's a version to try to avoid the memory errors from storing everything in a single list. Now, the features are added once per row. It's not as fast as the above method, but less memory hungry:

vl = QgsVectorLayer("Point?crs=", "temporary_points2", "memory")
vl.dataProvider().addAttributes([ QgsField("X_coord", QVariant.Int), QgsField("Y_coord", QVariant.Int), QgsField("Azimuth", QVariant.Int) ])

for row in range(0,rows):
   new_features = []
   for col in range(0,cols):
      rspnt = pixel2coord(col,row)
      fet = QgsFeature()
      fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(rspnt[0], rspnt[1])))
      fet.initAttributes(3)
      fet.setAttribute(0, rspnt[0] )
      fet.setAttribute(1, rspnt[1])
      az = pntxy.azimuth(QgsPoint(rspnt[0], rspnt[1]))
      if az < 0.0:
         azim = 360 + az
      else:
         azim = az
      fet.setAttribute(2, azim)
      new_features.append(fet)

   vl.dataProvider().addFeatures( new_features )
vl.updateExtents()
vl.updateFields()    
QgsMapLayerRegistry.instance().addMapLayers([vl])
  • Hi ndawson, you are a professional programmer, It's work very well, for 291*220 pixels, time change from 5 minutes to 8 seconds, for 418*326 pixels time is 14 seconds, It's very good,but when I want to run for 10393* 7170 pixels layer below error appear!!!! Traceback (most recent call last):........... new_features.append(fet) MemoryError. what's happend!!? – nickan Jan 27 '16 at 21:26
  • 1
    Try the updated version – ndawson Jan 27 '16 at 21:43
  • I used updated version, first time without any error, closed QGIS program, in second time, after 1 minute, appeared an "Exception window" with "bad allocation" message!! Is this relation to my Computer? My system is win 8.1 64bit, RAM 4 GB, Core i7 2.2 GHz. I installed QGIS 2.12.3 Lyon from OSGeo4W. – nickan Jan 27 '16 at 22:15
  • Ok, revision posted – ndawson Jan 27 '16 at 22:32
  • Hi ndawson, I run two times your revision posted, but after 3 minutes, without any error closed QGIS program! – nickan Jan 28 '16 at 5:58
4

1) If your pnt shapefile has only one element, use the iterator method next():

feat = pnt.getFeatures().next()
pntxy = feat.geometry().asPoint()

If the shapefile has many elements, use list comprehensions

pntxy = [feat.geometry().asPoint() for feat in pnt.getFeatures()]
# select pntxy[0] or pntxy[1] or...

2) instead of GDAL, why don't you use directly PyQGIS ?

ds = QgsRasterLayer('E:\Sample_Dataset_Lorestan\subset_a1.tif','subset_a1')
pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()
# extent of the layer
ext = ds.extent()
originX ,originY  = (ext.xMinimum(),ext.yMinimum())
cols = ds.width()
rows = ds.height()

3) memory layer in the Python console (no need of vl.startEditing())

vl = QgsVectorLayer("Point?crs=", "temporary_points2", "memory")
vl.dataProvider().addAttributes([ QgsField("X_coord", QVariant.Int), QgsField("Y_coord", QVariant.Int), QgsField("Azimuth", QVariant.Int) ])
fet = QgsFeature()
for row in range(0,rows):
   for col in range(0,cols):
      rspnt = pixel2coord(col,row)
      fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(rspnt[0], rspnt[1])))
      fet.initAttributes(3)
      fet.setAttribute(0, rspnt[0] )
      fet.setAttribute(1, rspnt[1])
      az = pntxy.azimuth(QgsPoint(rspnt[0], rspnt[1]))
      if az < 0.0:
         azim = 360 + az
      else:
         azim = az
      fet.setAttribute(2, azim)
      vl.dataProvider().addFeatures( [fet] )
      vl.updateExtents()
      vl.updateFields()

QgsMapLayerRegistry.instance().addMapLayers([vl])
  • That's great job, very thanks gene, new code change time of run for 291*220 pixels from 22 minutes to 5 minutes!!! But my main raster layer is 10393* 7170 pixels, It is important to reduce this time again. can I do another changes for receive lower time of run? Thanks. – nickan Jan 24 '16 at 17:50
  • You can use Fiona and Rasterio in Python only without PyQGIS (faster) – gene Jan 24 '16 at 18:04
  • Try moving the vl.updateExtents() and vl.updateFields() out of the loop. They don't need to be called multiple times. – ndawson Jan 25 '16 at 5:14
  • Another good optimisation would be not calling vl.dataProvider().addFeatures() for every single feature. Instead, add them all to a list and only call addFeatures once for the whole list. – ndawson Jan 25 '16 at 5:17
  • 1
    Well, in your current code you are adding a list of 1 feature for every pixel. Instead, add ALL the pixel features to a list and call addFeatures ONCE with the entire list of features. So 1 call to addFeatures, rather than 1000s. – ndawson Jan 25 '16 at 23:26

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.