Tuesday, June 29, 2021

GIS Programming with Python and QGIS - Part 3

GIS Programming with Python and QGIS - Part 3 >> The PyQGIS module

If you missed Part 1 and Part 2 read them before you continue reading...

In this part I will cover the following topics:-
1~ Introduction to PyQGIS
2~ Basic GIS Operations using pyqgis

Let's get started...


1~ Introduction to PyQGIS
PyQGIS is the official python module for the QGIS software. It allows you to perform spatial operations and automate spatial processes in the QGIS software.

The official python documentation is available on this website url (the C++ documentation is available here. Some notable resources for learning to use the module includes:-
i) PyQGIS developer cookbook - part of QGIS Documentation
ii) PyQGIS Programmer's Guide by Gary Sherman
iii) PyQGIS 101 by Anita Graser
iv) Open Source Option YouTube Channel on PyQGIS

The resources above are good enough to get you started.



2~ Basic GIS Operations using PyQGIS
Lets take a look at some common operations we usually do in QGIS software and see how we can do the same in PyQGIS.

i) Loading and rendering shapefile in PyQGIS
The two important methods needed to load a vector layer onto the map canvas are: QgsVectorLayer(vector_file, 'DISPLAYNAME', 'ogr') and QgsProject.instance().addMapLayer(layer)

See proper usage below:-

# Path to the shapefile...
vector_file = r"C:\Users\Yusuf_08039508010\Desktop\Working_Files\GIS Data\NGR\Nig Adm\NGA_adm2.shp"

# Create vector layer object...
layer_1 = QgsVectorLayer(vector_file, 'DISPLAYNAME', 'ogr')

# Add shapefile to map canvas...
QgsProject.instance().addMapLayer(layer_1)
If you want to see a more useful scenario of adding/loading vector layers into qgis layer panel using pyqgis see this page.


We can also interact with the project file (.qgs or .qgz) as follow:-
# Create project instance object
project = QgsProject.instance()

# Get project file path name
project.fileName()

# Get project file base name
project.baseName()

# Set project file to new title
project.setTitle('New Title Name...')

# Get lenght/number of layers in a project layer panel
len(project.mapLayers())

# Find more attribute and methods for QgsProject
dir(QgsProject.instance())




ii) Accessing Shapefile Attributes Table
To access the attribute table fields/column use the layer_1.fields() method like this in a forloop expression...

# get attribute columns names
for f in layer_1.fields():
    print(f.name())
This will printout all the field/column names. You can manipulate individual column from the temporary variable f.





iii) Zoom functions
You can perform many zoom types using pyqgis, let see some in action below.

If you check the dir(...) on the map canvas object, you will get the following list of zoom functions that can carryout various types of zooming on the map canvas.
'zoomByFactor', 'zoomIn', 'zoomLastStatusChanged', 'zoomNextStatusChanged', 'zoomOut', 'zoomResolutions', 'zoomScale', 'zoomToFeatureExtent', 'zoomToFeatureIds', 'zoomToFullExtent', 'zoomToNextExtent', 'zoomToPreviousExtent', 'zoomToSelected', 'zoomWithCenter'


canvas = iface.mapCanvas() # Or: canvas = qgis.utils.iface.mapCanvas()
canvas.zoomIn()    

canvas.zoomOut()

canvas.zoomByFactor(10) # zoomout by factor of 10

canvas.zoomToFullExtent()

canvas.zoomScale(1000000) # zoom to scale 1:1000000

# Get cuurent canvas extent in the form of QgsRectangle
e = iface.mapCanvas().extent()
canvas.zoomToFeatureExtent(e) # it takes in QgsRectangle

# Just as the name implies zoomToPreviousExtent
canvas.zoomToPreviousExtent()

canvas.zoomWithCenter(7, 8, True)
canvas.zoomWithCenter(7, 8, False)

canvas.zoomToSelected()
Above are a handful of zoom implementations in pyqgis. You can use the doc together with dir() and help() function to lookup expected parameters for some of the function.




iv) Rename layer
# Get active layer...
layer = iface.activeLayer()

# Display name...
layer.name()

# Set new name...
layer.setName('NewName')


v) Making attribute selection
# Get active layer...
layer = iface.activeLayer()

# Make selection...
layer.selectByExpression( " state_name = 'Kogi' " ) # Or ( " \"state_name\" = 'Kaduna' " )

You can use this selection concept to make a game. I have done it! I named it KnowYourState (KYS), see the code on this page.



vi) Change Canvas Color, Title and Overall Look & Feel
from PyQt5 import QtTest

for _ in range(5):
    app = QApplication.instance()
    app.setStyleSheet("QWidget {color: blue; background-color: yellow;}")
    iface.mainWindow().setWindowTitle("My QGIS")
    
    print("Change to: Black color")
    iface.mapCanvas().setCanvasColor(Qt.black)
    iface.mapCanvas().refresh()
    
#    time.sleep(3)
    QtTest.QTest.qWait(1000)
    
    
    app = QApplication.instance()
    app.setStyleSheet("QWidget {color: blue; background-color: green;}")
    iface.mainWindow().setWindowTitle("I love QGIS")
    
    print("Change to: Red color")
    iface.mapCanvas().setCanvasColor(Qt.red)
    iface.mapCanvas().refresh()
    print("-"*10)
    
    QtTest.QTest.qWait(1000)


Note that time.sleep() causes QGIS canvas to freeze, so I used PyQt5 built-in timer for the delay sequence.



vii) Changing layer visibility to On or Off
# Uncheck or Hide layer visibility...

qgis_prjt_lyrs = QgsProject.instance().layerTreeRoot().findLayers()
qgis_prjt_lyr_ids = QgsProject.instance().layerTreeRoot().findLayerIds()


for ly_name, ly_id in zip(qgis_prjt_lyrs, qgis_prjt_lyr_ids):
    print(ly_name.name())
    print(ly_id)

    # By layer name...
    ly_name.setItemVisibilityChecked(True) # True=On, False=Off

    # By layer Id...
    # QgsProject.instance().layerTreeRoot().findLayer(ly_id).setItemVisibilityChecked(False)



viii) Renaming attribute table field/column
# https://gis.stackexchange.com/questions/241149/renaming-field-in-attribute-table-via-pyqgis/359929
# Delete attribute column
vlayer = qgis.utils.iface.activeLayer()

# Open editing session
vlayer.startEditing()

# Rename field
for field in vlayer.fields():
    if field.name() == 'Old_Field_Name':
        idx = vlayer.fields().indexFromName(field.name())
        vlayer.renameAttribute(idx, 'New_Field_Name')

# Close editing session and save changes
vlayer.commitChanges()



ix) Load raster layer and read its properties

Loading a raster layer is similar to loading vector layer.

# Load raster layer...
path = r"C:\Users\Yusuf_08039508010\Desktop\IMG\n08_e007_3arc_v2.tif"

rlayer = QgsRasterLayer(path, 'Raster_Layer_Name', 'gdal')
QgsProject.instance().addMapLayer(rlayer)

Or
# Load raster layer...
path = r"C:\Users\Yusuf_08039508010\Desktop\IMG\n08_e007_3arc_v2.tif"
iface.addRasterLayer(path,"Raster_Layer_Name","gdal")
Note that layers such as OSM and other web map services (WMS) can be added via the raster layer. See more raster driver layers on this page.

With the rsater layer loaded, you can access many raster methods and properties as available when you do a dir(...) on the raster layer object.

>>> my_raster = iface.activeLayer()

>>> my_raster
<QgsMapLayer: 'Raster_Layer_Name' (gdal)>

>>> dir(my_raster)
['Actions', 'AllStyleCategories', 'AttributeTable', 'ColorLayer', 'ColorRampShader', 'ColorShadingAlgorithm', 'CustomProperties', 'Diagrams', 'Fields', 'FlagDontResolveLayers', 'Forms', 'FreakOutShader', 'GeometryOptions', 'GrayOrUndefined', 'Identifiable', 'Labeling', 'LayerConfiguration', 'LayerFlag', 'LayerFlags', 'LayerOptions', 'LayerType', 'MULTIPLE_BAND_MULTI_BYTE_ENHANCEMENT_ALGORITHM', 'MULTIPLE_BAND_MULTI_BYTE_MIN_MAX_LIMITS', 'MULTIPLE_BAND_SINGLE_BYTE_ENHANCEMENT_ALGORITHM', 'MULTIPLE_BAND_SINGLE_BYTE_MIN_MAX_LIMITS', 'MapTips', 'MeshLayer', 'Metadata', 'Multiband', 'Palette', 'PluginLayer', 'PropertyType', 'PseudoColorShader', 'RasterLayer', 'ReadFlag', 'ReadFlags', 'Relations', 'Removable', 'Rendering', 'SAMPLE_SIZE', 'SINGLE_BAND_ENHANCEMENT_ALGORITHM', 'SINGLE_BAND_MIN_MAX_LIMITS', 'Searchable', 'Style', 'StyleCategories', 'StyleCategory', 'Symbology', 'Symbology3D', 'UndefinedShader', 'UserDefinedShader', 'VectorLayer', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 'abstract', 'accept', 'appendError', 'attribution', 'attributionUrl', 'autoRefreshInterval', 'autoRefreshIntervalChanged', 'bandCount', 'bandName', 'beforeResolveReferences', 'blendMode', 'blendModeChanged', 'blockSignals', 'brightnessFilter', 'childEvent', 'children', 'clone', 'configChanged', 'connectNotify', 'constDataProvider', 'createMapRenderer', 'crs', 'crsChanged', 'customEvent', 'customProperty', 'customPropertyKeys', 'dataChanged', 'dataProvider', 'dataSourceChanged', 'dataUrl', 'dataUrlFormat', 'decodedSource', 'deleteLater', 'dependencies', 'dependenciesChanged', 'destroyed', 'disconnect', 'disconnectNotify', 'draw', 'dumpObjectInfo', 'dumpObjectTree', 'dynamicPropertyNames', 'emitStyleChanged', 'encodedSource', 'error', 'event', 'eventFilter', 'exportNamedMetadata', 'exportNamedStyle', 'exportSldStyle', 'extensionPropertyType', 'extent', 'findChild', 'findChildren', 'flags', 'flagsChanged', 'formatLayerName', 'generateId', 'hasAutoRefreshEnabled', 'hasDependencyCycle', 'hasScaleBasedVisibility', 'height', 'htmlMetadata', 'hueSaturationFilter', 'id', 'ignoreExtents', 'importNamedMetadata', 'importNamedStyle', 'inherits', 'installEventFilter', 'isEditable', 'isInScaleRange', 'isRefreshOnNotifyEnabled', 'isSignalConnected', 'isSpatial', 'isTemporary', 'isValid', 'isValidRasterFileName', 'isWidgetType', 'isWindowType', 'keywordList', 'killTimer', 'lastModified', 'legend', 'legendChanged', 'legendSymbologyItems', 'legendUrl', 'legendUrlFormat', 'loadDefaultMetadata', 'loadDefaultStyle', 'loadNamedMetadata', 'loadNamedMetadataFromDatabase', 'loadNamedStyle', 'loadNamedStyleFromDatabase', 'loadSldStyle', 'maximumScale', 'metaObject', 'metadata', 'metadataChanged', 'metadataUri', 'metadataUrl', 'metadataUrlFormat', 'metadataUrlType', 'minimumScale', 'moveToThread', 'name', 'nameChanged', 'objectName', 'objectNameChanged', 'originalXmlProperties', 'paletteAsPixmap', 'parent', 'pipe', 'previewAsImage', 'property', 'providerType', 'publicSource', 'pyqtConfigure', 'rasterType', 'rasterUnitsPerPixelX', 'rasterUnitsPerPixelY', 'readCommonStyle', 'readCustomProperties', 'readLayerXml', 'readOnly', 'readSld', 'readStyle', 'readStyleManager', 'readSymbology', 'readXml', 'recalculateExtents', 'receivers', 'refreshOnNotifyMessage', 'reload', 'removeCustomProperty', 'removeEventFilter', 'renderer', 'renderer3D', 'renderer3DChanged', 'rendererChanged', 'repaintRequested', 'resampleFilter', 'resolveReferences', 'saveDefaultMetadata', 'saveDefaultStyle', 'saveNamedMetadata', 'saveNamedStyle', 'saveSldStyle', 'sender', 'senderSignalIndex', 'setAbstract', 'setAttribution', 'setAttributionUrl', 'setAutoRefreshEnabled', 'setAutoRefreshInterval', 'setBlendMode', 'setContrastEnhancement', 'setCrs', 'setCustomProperties', 'setCustomProperty', 'setDataProvider', 'setDataSource', 'setDataUrl', 'setDataUrlFormat', 'setDefaultContrastEnhancement', 'setDependencies', 'setError', 'setExtent', 'setFlags', 'setKeywordList', 'setLayerOrder', 'setLegend', 'setLegendUrl', 'setLegendUrlFormat', 'setMaximumScale', 'setMetadata', 'setMetadataUrl', 'setMetadataUrlFormat', 'setMetadataUrlType', 'setMinimumScale', 'setName', 'setObjectName', 'setOriginalXmlProperties', 'setParent', 'setProperty', 'setProviderType', 'setRefreshOnNofifyMessage', 'setRefreshOnNotifyEnabled', 'setRenderer', 'setRenderer3D', 'setScaleBasedVisibility', 'setShortName', 'setSubLayerVisibility', 'setTitle', 'setTransformContext', 'setValid', 'shortName', 'showStatusMessage', 'signalsBlocked', 'source', 'startTimer', 'staticMetaObject', 'statusChanged', 'styleChanged', 'styleLoaded', 'styleManager', 'styleURI', 'subLayers', 'thread', 'timerEvent', 'timestamp', 'title', 'tr', 'transformContext', 'triggerRepaint', 'type', 'undoStack', 'undoStackStyles', 'width', 'willBeDeleted', 'writeCommonStyle', 'writeCustomProperties', 'writeLayerXml', 'writeSld', 'writeStyle', 'writeStyleManager', 'writeSymbology', 'writeXml']


Some example are:-
1) Read raster width and height
my_raster.width()
my_raster.height()

2) Read raster Extent
my_raster.extent().toString()

3) Read raster statistics
stats = my_raster.dataProvider().bandStatistics(1)

stats.Range
stats.elementCount
stats.maximumValue
stats.minimumValue
stats.Max
stats.Min
stats.stdDev
stats.statsGathered
stats.sum
stats.sumOfSquares
stats.Mean
stats.mean





x) Delete columns from attribute table

The code below will delete multiple columns from multiple layers.
from qgis.core import QgsProject
layersDict = QgsProject.instance().mapLayers()

# Get attribute column index...
for (id, map) in layersDict.items():
    if (map.type() == QgsMapLayer.VectorLayer):
        col_index_1 = map.fields().indexFromName('source')
        col_index_2 = map.fields().indexFromName('global_id')
        col_index_3 = map.fields().indexFromName('layer')
        col_index_4 = map.fields().indexFromName('path')
        
        print(map.name(), col_index_1, col_index_2, col_index_3, col_index_4)
    
        # Delete by index of the columns...
        map.dataProvider().deleteAttributes([col_index_1, col_index_2, col_index_3, col_index_4])
        map.updateFields() # Save changes made to the attribute table



xi) Create New Vector shapefile layer
file_path = r'C:\Users\Yusuf_08039508010\Desktop\new shp\A_NEW_SHP.shp'

# Create a Attribute fields 
field = QgsFields()
field.append( QgsField('id', QVariant.Int) ) # Integer field
field.append( QgsField('name', QVariant.String) ) # Str field
field.append( QgsField('age', QVariant.Double) ) # Double

# Make a point layer writer obj
writer = QgsVectorFileWriter(file_path, 'utf-8', field, QgsWkbTypes.Point, QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')

# Create features
feat = QgsFeature()
feat.setGeometry( QgsGeometry.fromPointXY( QgsPointXY(5.9, 9.2) ) )
feat.setAttributes( [0, 'Umar', 23.4] )

# Commit/Add the feature...
writer.addFeature(feat)
del(writer)


iface.addVectorLayer(file_path, '', 'ogr')

The above code is for a point shapefile. To create a line and polygon shapefile, you need to change the write object and create features object.

For Line:
file_path = r'C:\Users\Yusuf_08039508010\Desktop\Working_Files\Fiverr\2021\07-July\Shittu Map\new shp\A_NEW_SHP-2.shp'

# Create a Attribute fields 
field = QgsFields()
field.append( QgsField('id', QVariant.Int) ) # Integer field
field.append( QgsField('name', QVariant.String) ) # Str field
field.append( QgsField('age', QVariant.Double) ) # Double

# Make a line layer writer obj
writer = QgsVectorFileWriter(file_path, 'utf-8', field, QgsWkbTypes.LineString, QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')


# Create features
feat = QgsFeature()
p_points=[QgsPointXY(6.68, 8.91), QgsPointXY(7.30, 8.83), QgsPointXY(7.28, 8.56)]
feat.setGeometry( QgsGeometry.fromPolyline(p_points))
feat.setAttributes( [0, 'Umar', 23.4] )
feat..setWidth(3)


# Commit/Add the feature...
writer.addFeature(feat)
del(writer)


iface.addVectorLayer(file_path, '', 'ogr')

For Polygon:
file_path = r'C:\Users\Yusuf_08039508010\Desktop\Working_Files\Fiverr\2021\07-July\Shittu Map\new shp\A_NEW_SHP-3.shp'

# Create a Attribute fields 
field = QgsFields()
field.append( QgsField('id', QVariant.Int) ) # Integer field
field.append( QgsField('name', QVariant.String) ) # Str field
field.append( QgsField('age', QVariant.Double) ) # Double

# Make a polygon layer writer obj
writer = QgsVectorFileWriter(file_path, 'utf-8', field, QgsWkbTypes.Polygon, QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')


# Create features
feat = QgsFeature()
p_points=[QgsPointXY(6.68, 8.91), QgsPointXY(7.30, 8.83), QgsPointXY(7.28, 8.56), QgsPointXY(7.10, 8.45), QgsPointXY(6.61, 8.49)]
feat.setGeometry( QgsGeometry.fromPolygonXY(p_points))
feat.setAttributes( [0, 'Umar', 23.4] )
feat.setColor(QColor(0, 0, 255))
feat.setFillColor(QColor(255,255,0))

# Commit/Add the feature...
writer.addFeature(feat)
del(writer)


iface.addVectorLayer(file_path, '', 'ogr')





xii) Add attribute to Vector shapefile layer
# Adding new attribute FIELD to SHAPEFILES...
from PyQt5.QtCore import QVariant

# Get layers added on the layers panel...
# Get path to layer location on disc...
layers_on_panel = QgsProject.instance().mapLayers()
layer_paths = [l.source() for l in layers_on_panel.values()]


for l in layer_paths:
    displayName = l.split('\\')[-1].replace('.shp', '')
    layer = QgsVectorLayer(l, f"{displayName}", "ogr")
    
    # Define dataProvider for layer
    layer_provider = layer.dataProvider()
    
    # Add string attribute field and update fields...
    layer_provider.addAttributes([QgsField("NewName", QVariant.String)])
    layer.updateFields()
    
    # Print fields in layer...
    print (layer.fields().names())

print('Done...')



xiii) Updating attribute field of a shapefile layer
# Updating Field's Attributes

# Get layers added on the layers panel...
layers_on_panel = QgsProject.instance().mapLayers()

# for k, v in layers_on_panel.items():
for k, v in zip(layers_on_panel.keys(), layers_on_panel.values()):
    layer_name = v.name()
    print('Processing....', layer_name)
    
    # Read the map layer...
    map_file = layers_on_panel[k]
    layer = QgsVectorLayer(map_file.source(), "ogr")

    # Get fields in layer...
    fieldnames = layer.fields().names()
    print (fieldnames)
    
    # get the index of 'NewName' field
    index = fieldnames.index('NewName')    
    print(index)
    
    
    # Define dataProvider for layer
    layer_provider = layer.dataProvider()

    features = layer.getFeatures()
    
    layer.startEditing()
    
    for f in features:
        id = f.id()
        # length = f.geometry().length()
        attr_value = {index:f'{layer_name}'}
        layer_provider.changeAttributeValues({id:attr_value})
    layer.commitChanges()

    # break

print('Finished....')



xiv) Save map canvas to image
# Save map canvas as image, a quick export of the current map canvas
iface.mapCanvas().resize(QSize(1280, 720)) # set canvas size
iface.mapCanvas().saveAsImage('Test123.png', None, 'PNG') # save img


Conclusion

This is not by anyway the only operations you can do with pyqgis API. As you can see, we can do almost anything in QGIS with the PyQGIS API. You can extend the capabilities of QGIS beyond imaginations.

Keep imaging, keep coding.

No comments:

Post a Comment