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 6752 pixels in 30 meter resolution time for running is 5 seconds, for 133105 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 scrip for large data? Thanks.

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 6752 pixels in 30 meter resolution time for running is 5 seconds, for 133105 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?

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"

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"