How to speed up this script in PyQGIS?

Question

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"

Answer 1

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])

Answer 2

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])