PyQGIS introduction

Let’s start with simple tasks.

First steps

Layers access

Task

• list of layers determined from current project (even unsaved) on line 1

• print layer name and feature count for vector layers (3)

• for other layer types print only name (5)

for layer in QgsProject.instance().mapLayers().values():
    if layer.type() == QgsMapLayer.VectorLayer:
        print (layer.name(), layer.featureCount())
    else:
        print (layer.name())

Script to download.

Fig. 3 Example of running simple script from QGIS built-in editor.

Task

• skip disabled layers on line 2

• process only vector layers with multipolygon geometry type (5)

• loop over features on line 9

• get feature geometry on line 10

• get attribute value nazev (name in Czech) if exists (12)

• print feature info (fid, name, area in hectares - assuming map units meters), see line 16

for layer in QgsProject.instance().mapLayers().values():
    if not QgsProject.instance().layerTreeRoot().findLayer(layer.id()).itemVisibilityChecked():
        continue

    if layer.type() != QgsMapLayer.VectorLayer or layer.wkbType() != QgsWkbTypes.MultiPolygon:
        continue
  
    print (layer.name(), layer.featureCount())
    for feat in layer.getFeatures():
        geom = feat.geometry()
        try:
            nazev = feat['nazev']
        except KeyError:
            nazev = 'no name defined'
            
        print ('{0}: {1} {2:.1f} ha'.format(feat.id(), nazev, geom.area()/1e4))

Script to download.

Vector data, features access and editing

Sample data (ice cream shops in Czech Republic downloaded from OpenStreetMap as points only): ice_cream.csv.

Adding CSV-based layers

Task

• load layer without geometry (as table) (see lines 2 and 5)

• add table layer into layer tree (6)

filename = '/home/martin/Downloads/ice_cream.csv'
uri = 'file:///{}?delimiter=,'.format(filename)

name = os.path.splitext(os.path.basename(filename))[0]
layer = QgsVectorLayer(uri, name, 'delimitedtext')

QgsProject.instance().addMapLayer(layer)

Script to download.

Task

• geometry built from longitute and latitude

• set spatial reference system to EPSG:4326

• add map layer into layer tree

filename = '/home/martin/Downloads/ice_cream.csv'
uri = 'file:///{}?delimiter=,&xField=lon&yField=lat&crs=epsg:4326'.format(filename)

name = os.path.splitext(os.path.basename(filename))[0]
layer = QgsVectorLayer(uri, name, 'delimitedtext')

QgsProject.instance().addMapLayer(layer)

Script to download.

Save layer into file-based GIS format

Task

• reproject to S-JTSK (EPSG:5514) (target CRS set on line 4)

• save layer into Esri Shapefile format (this task performed on lines 5-7)

layer = QgsProject.instance().mapLayersByName('ice_cream')[0]

shp_file = "/tmp/ice_cream.shp"
crs = QgsCoordinateReferenceSystem("EPSG:5514")
QgsVectorFileWriter.writeAsVectorFormat(
    layer, shp_file,
    "UTF-8", driverName="ESRI Shapefile", destCRS=crs)
                                        
layer_shp = QgsVectorLayer(shp_file, "test", "ogr")
print(layer_shp.isValid())

Script to download.

Task

Improvements

• check if input layer found (see line 1)

• check for writter’s errors (see cookbook example)

Editing vector features

Let’s delete all features (ice cream shows) located more than 1km from specified point.

Task

• select features in distance of 1km from current position (14.3881100E, 50.1041200N), see lines 8 and 21

• point of interest must be transformed into layer’s CRS (10-13)

• area of interest is defined by rectangle on lines 15-19

• invert selection on line 24

• delete selected features on line 27

• note that layer must be switched to edinging mode, it can be done eg. by triggers, see lines 26,28

shp_file = "/tmp/ice_cream.shp"
layer = QgsVectorLayer(shp_file, "test", "ogr")

QgsProject.instance().addMapLayer(layer)

x = 14.3881100
y = 50.1041200
geom = QgsGeometry.fromPointXY(QgsPointXY(x, y))

p_crs = QgsCoordinateReferenceSystem("EPSG:4326")
trans = QgsCoordinateTransform(
    p_crs, layer.crs(), QgsProject.instance())
geom.transform(trans)

p = geom.asPoint()
offset = 1500
p1 = QgsPointXY(p.x() - offset, p.y() - offset)
p2 = QgsPointXY(p.x() + offset, p.y() + offset)
aoi = QgsRectangle(p1, p2)

layer.selectByRect(aoi)
print (len(layer.selectedFeatures()))

layer.invertSelection()

iface.actionToggleEditing().trigger()
iface.actionDeleteSelected().trigger()
iface.actionToggleEditing().trigger()

Script to download.

Task

Compare with

layer.startEditing()
layer.deleteSelectedFeatures()
layer.commitChanges()