Python for GIS: The QGIS Python Console, PyQGIS, GeoPandas and GDAL

# Run this in the QGIS Python Console (Plugins > Python Console). "iface" is QGIS itself, already available there.

# Step 1: load a file into the QGIS project as a layer called "places".
# The r before the quote just means "read the backslashes and slashes literally".
layer = iface.addVectorLayer(r"C:/Users/bhatt/Documents/code/timeandspatial-tutorial-data/vector/places_full.gpkg", "places", "ogr")

# Step 2: ask the layer about itself: name, number of features, coordinate system
print(layer.name(), layer.featureCount(), layer.crs().authid())

# Step 3: list the field (column) names
print([f.name() for f in layer.fields()])

# Step 4: go through the features (rows) one at a time and print three values from each
for f in layer.getFeatures():
    print(f["id"], f["name"], f["category"])

# Step 1: select the features whose category is "park" (the same expression language as the QGIS select box)
layer.selectByExpression('"category" = \'park\'')
print(layer.selectedFeatureCount())          # 3

# Step 2: edit the layer. First switch editing on, and find which column number "category" is.
layer.startEditing()
idx = layer.fields().indexFromName("category")

# Step 3: change the category of every selected feature
for f in layer.selectedFeatures():
    layer.changeAttributeValue(f.id(), idx, "green space")

# Step 4: undo everything. Use commitChanges() instead if you want to keep the edit.
layer.rollBack()                              # use commitChanges() to keep the edit

# Processing = the same tools you find in the QGIS Processing Toolbox, but called from code.
import processing

# Step 1: reproject the layer to metres (UTM 43N), so the buffer distance is in metres.
# run(tool name, settings)  ->  returns a dictionary. ["OUTPUT"] picks the result layer.
# "memory:" means "keep the result in memory, do not save a file".
utm = processing.run("native:reprojectlayer", {
    "INPUT": layer,
    "TARGET_CRS": "EPSG:32643",
    "OUTPUT": "memory:",
})["OUTPUT"]

# Step 2: draw a 2 km buffer around every place
buffers = processing.run("native:buffer", {
    "INPUT": utm,                  # the layer from step 1
    "DISTANCE": 2000,              # metres
    "SEGMENTS": 16,                # how smooth the circle is
    "DISSOLVE": False,             # keep one buffer per place
    "OUTPUT": "memory:",
})["OUTPUT"]

# Step 3: name the result and add it to the project so it shows on the map
buffers.setName("buffer_2km")
QgsProject.instance().addMapLayer(buffers)

# This version runs OUTSIDE QGIS, as a normal Python script. It needs a little setup that the console does for you.
import sys
from qgis.core import QgsApplication, QgsVectorLayer

# Step 1: tell Python where QGIS is installed, then start QGIS in the background (no window)
QgsApplication.setPrefixPath(r"C:/Program Files/QGIS 3.44.13/apps/qgis-ltr", True)
qgs = QgsApplication([], False)       # False = no screen
qgs.initQgis()

# Step 2: turn on the Processing tools
sys.path.append(r"C:/Program Files/QGIS 3.44.13/apps/qgis-ltr/python/plugins")
import processing
from processing.core.Processing import Processing
Processing.initialize()

# Step 3: open the layer straight from the file (there is no QGIS window, so no iface)
layer = QgsVectorLayer("C:/Users/bhatt/Documents/code/timeandspatial-tutorial-data/vector/places_full.gpkg", "places", "ogr")
print(layer.isValid(), layer.featureCount())      # True means the file opened

# Step 4: run the same reproject tool as before
result = processing.run("native:reprojectlayer",
                        {"INPUT": layer, "TARGET_CRS": "EPSG:32643", "OUTPUT": "memory:"})
print(result["OUTPUT"].crs().authid())

# Step 5: shut QGIS down cleanly when finished
qgs.exitQgis()

# The same job in GeoPandas: no QGIS needed, just a few lines.
import geopandas as gpd

root = "C:/Users/bhatt/Documents/code/timeandspatial-tutorial-data"

# Step 1: read the file
places = gpd.read_file(root + "/vector/places_full.gpkg")
print(places.crs, len(places))                    # EPSG:4326 10

# Step 2: convert to metres and keep only the parks
places_utm = places.to_crs(32643)                 # metres
parks = places_utm[places_utm["category"] == "park"]
print(parks[["name", "practice_visitors_per_day"]])

# Step 3: distance from the 4th place (position 3, counting from 0) to every place, in metres
print(places_utm.geometry.distance(places_utm.geometry.iloc[3]).round(0).tolist())

# Step 4: 2 km buffers. Copy first so the original points are untouched.
buffers = places_utm.copy()
buffers["geometry"] = places_utm.buffer(2000)

# Step 5: save the parks to a new file
parks.to_file(root + "/outputs/parks_demo.gpkg", driver="GPKG")

# The same job with GDAL, on a raster (an NDVI image)
from osgeo import gdal
import numpy as np

gdal.UseExceptions()

root = "C:/Users/bhatt/Documents/code/timeandspatial-tutorial-data"

# Step 1: open the file and take band 1
ds = gdal.Open(root + "/outputs/city/ndvi_2025.tif")
band = ds.GetRasterBand(1)

# Step 2: read the values, and drop the cells marked as "no data"
arr = band.ReadAsArray().astype("float32")
nodata = band.GetNoDataValue()
valid = arr[arr != nodata]

# Step 3: summary numbers
print(arr.shape)                                   # (1600, 2400)  rows x columns
print(round(float(valid.mean()), 3))               # 0.278         average NDVI
print(round(float((valid >= 0.5).mean()) * 100, 1))   # 17.7       percent of cells with NDVI 0.5 or more

# Step 4: where the grid sits on the map
gt = ds.GetGeoTransform()
print(gt[1], gt[0], gt[3])                         # cell size, left, top