¿Cómo acelerar este script en QGIS?

Question

Soy nuevo en PyQGIS, quiero crear una capa de puntos a partir de una capa raster. En la tabla de la capa de puntos creado, hay 3 campos por el nombre de "X_coord", "Y_coord" y "Azimut", cada punto en el Punto de la capa de muestra x ("X_coord de campo"), y ("Y_coord") centro de coordenadas de cada píxel de la capa raster y "Azimut" campo muestra el azimut desde otro punto único de la capa para cada centro de píxeles. Al ejecutar esta secuencia de comandos para la capa de trama con 67*52 píxeles de 30 metros de resolución tiempo de ejecución es de 5 segundos, para 133*105 píxeles de tiempo es de 15 segundos y para 291*220 píxeles tiempo de 22 minutos!!!! Creo que esto es inusual.¿Cómo puedo reducir el tiempo para ejecutar este scrip de datos de gran tamaño? Gracias.

Mi script es:

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) Si su pnt shapefile tiene sólo un elemento, utilice el iterador método next():

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

Si el archivo de forma tiene muchos elementos, utilice la lista de comprensión

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

2) en lugar de GDAL, ¿por qué no usar directamente 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) capa de memoria en la consola de Python (no hay necesidad de 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

Intente esta versión (modificada de la respuesta del gene):

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

Esta versión sólo llama a addFeatures una vez, en lugar de una vez por píxel en el bucle.

Actualización: Aquí hay una versión para tratar de evitar los errores de memoria de almacenar todo en una sola lista. Ahora, las características se agregan una vez por fila. No es tan rápido como el método anterior, pero menos hambre de memoria:

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