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