- 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...
- 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()
- 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(rspnt2QgsPoint(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])