Python for GIS: Automating Workflows with GDAL & Geopandas
Doing repetitive GIS tasks by hand (clipping 50 files, calculating statistics for hundreds of features) is slow and error-prone. **Python** lets you automate these workflows using powerful open-source libraries — turning hours of manual clicking into a script that runs in seconds. ---
📋 Prerequisites
- Python 3.8+ installed (via Conda or pip).
- Basic knowledge of Python scripting and syntax.
Step-by-Step Instructions
Set Up Your Python Environment
1. Install **Python** (3.9+ recommended) 2. Install core GIS libraries using pip: 3. (Optional) Use **Jupyter Notebook** or **VS Code** for an interactive coding experience ---
pip install geopandas rasterio gdal matplotlibLearn the Core Libraries
| Library | Purpose | |---|---| | **GDAL** | The foundational library for reading/writing raster and vector formats | | **Geopandas** | Pandas-like data handling, but for vector/spatial data | | **Rasterio** | Simplified raster reading, writing, and processing | | **Shapely** | Geometric operations (buffer, intersection, etc.) | | **Fiona** | Low-level vector file I/O (used internally by Geopandas) | ---
Read Vector Data with Geopandas
---
import geopandas as gpd
# Load a shapefile
gdf = gpd.read_file("study_area.shp")
# View the first few rows
print(gdf.head())
# Check the CRS
print(gdf.crs)Perform Spatial Operations
---
# Reproject to a projected CRS (for accurate area/distance)
gdf_utm = gdf.to_crs(epsg=32643) # example: UTM Zone 43N
# Calculate area in square meters
gdf_utm["area_m2"] = gdf_utm.geometry.area
# Create a 500m buffer around each feature
gdf_utm["buffer"] = gdf_utm.geometry.buffer(500)
# Save the result
gdf_utm.to_file("output_with_buffer.shp")Work with Raster Data using Rasterio
---
import rasterio
import numpy as np
with rasterio.open("landsat_image.tif") as src:
red = src.read(4).astype(float)
nir = src.read(5).astype(float)
profile = src.profile
# Calculate NDVI
ndvi = (nir - red) / (nir + red)
# Save the output
profile.update(dtype=rasterio.float32, count=1)
with rasterio.open("ndvi_output.tif", "w", **profile) as dst:
dst.write(ndvi.astype(rasterio.float32), 1)Automate a Batch Task
Example: Clip 20 raster files to the same boundary automatically. ---
import glob
import rasterio
from rasterio.mask import mask
import geopandas as gpd
boundary = gpd.read_file("boundary.shp")
geom = [boundary.geometry.iloc[0]]
for file in glob.glob("rasters/*.tif"):
with rasterio.open(file) as src:
out_image, out_transform = mask(src, geom, crop=True)
out_meta = src.meta.copy()
out_meta.update({"height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform})
output_path = file.replace("rasters/", "clipped/")
with rasterio.open(output_path, "w", **out_meta) as dst:
dst.write(out_image)
print("All rasters clipped successfully!")Visualize Results with Matplotlib
---
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8,8))
gdf.plot(ax=ax, column="area_m2", cmap="viridis", legend=True)
plt.title("Study Area by Size")
plt.show()🧠 Quick Recap & Practice Questions
## Quick Recap
- GDAL is the foundation; Geopandas handles vector, Rasterio handles raster
- Python scripts can automate repetitive tasks (clipping, reprojecting, calculating indices) across many files at once
- Always reproject to a projected CRS before calculating area or distance
- Combine libraries (Geopandas + Rasterio + Matplotlib) for a full analysis-to-visualization workflow
---
## Practice Questions
1. Why do you need to reproject data before calculating area in Geopandas?
2. Which library would you use to read and process a satellite raster band in Python?
3. What's one advantage of scripting a batch clip operation instead of doing it manually in QGIS?
---
*This completes the full 20-part GIS & Remote Sensing tutorial series — from GIS fundamentals to Python automation.*