Automating Repetitive GIS Work with Python and GeoPandas

Most GIS work repeats. A new batch of district layers arrives every month, each one needs the same CRS check, the same clip, the same attribute added. Doing that by hand in a desktop GIS takes an hour, and the one time you skip a step is the time it matters.

This post builds a short Python script that does the whole routine in one run.

What is GIS automation?

Automation here means writing the steps down as a script so the computer repeats them identically. You are not replacing the GIS software. QGIS and ArcGIS are still the right tools for looking at data and making maps. A script is better for steps that are the same every time.

The libraries used here:

• GeoPandas reads vector data (shapefile, GeoPackage, GeoJSON) into a table where each row has a geometry. It is built on pandas and Shapely.
• GDAL/OGR is the library underneath that does the actual file reading and writing.

Why it matters

• Repeatable. The same input gives the same output every time.
• Checkable. The script is a record of exactly what was done to the data, which is useful when someone asks how an area figure was produced.
• Faster at scale. Ten layers or a hundred take about the same effort to write.

Not everything is worth automating. If you will do a task once, do it by hand. If you do it monthly, or if mistakes are costly, script it.

How it works

A GeoPandas layer is a GeoDataFrame: a table with one column, geometry, that holds the shapes. Everything you know from tables works (filter, add columns, join), and you gain spatial operations like to_crs() for reprojecting and clip() for cutting to a boundary.

The CRS is stored on the GeoDataFrame as gdf.crs. This matters for the script below, because two checks depend on it:

• Layers must share a CRS before you clip one with another.
• Area in square metres only makes sense in a projected CRS. A geographic CRS such as EPSG:4326 measures in degrees.

Practical GIS example

You receive land-use shapefiles for several talukas in Pune district. Each one might be in a different CRS. You want one GeoPackage, in EPSG:32643 (WGS 84 / UTM zone 43N), clipped to the district boundary, with an area column in hectares.

Workflow

1. Put the shapefiles in a folder called raw/ and the district boundary in district.gpkg.
2. For each file: read it, report its CRS, reproject to EPSG:32643.
3. Clip to the district boundary.
4. Add an area column in hectares.
5. Write every result into one GeoPackage, one layer per input.

# Clean, clip and combine a folder of shapefiles into one GeoPackage.
# Folder layout:  raw/*.shp   and   district.gpkg  (the clip boundary)

from pathlib import Path
import geopandas as gpd

TARGET_CRS = "EPSG:32643"   # WGS 84 / UTM zone 43N, metres
OUTPUT = "clean_layers.gpkg"

# Step 1: Read the boundary once and put it in the target CRS
boundary = gpd.read_file("district.gpkg").to_crs(TARGET_CRS)

# Step 2: Loop over every shapefile in the raw folder
for path in sorted(Path("raw").glob("*.shp")):
    gdf = gpd.read_file(path)

    # Step 3: Stop on a missing CRS. Guessing here would hide an error.
    if gdf.crs is None:
        print(f"SKIPPED {path.name}: no CRS defined")
        continue

    # Step 4: Reproject to the target CRS (recalculates the coordinates)
    original = gdf.crs.to_string()
    gdf = gdf.to_crs(TARGET_CRS)

    # Step 5: Clip to the district boundary
    gdf = gpd.clip(gdf, boundary)

    # Step 6: Area in hectares (1 hectare = 10,000 square metres).
    # This is valid because TARGET_CRS is projected and uses metres.
    gdf["area_ha"] = gdf.geometry.area / 10_000

    # Step 7: Write one layer per input file into the GeoPackage
    gdf.to_file(OUTPUT, layer=path.stem, driver="GPKG")

    print(f"OK {path.name}: {original} -> {TARGET_CRS}, {len(gdf)} features")

Run it with python clean_layers.py. The printed line for each file records what CRS it started in, which is the kind of detail that is easy to forget doing it by hand.

Two decisions in the script are worth noticing. The script skips a layer with no CRS rather than assigning one, because a wrong guess produces a layer in the wrong place that looks fine in a table. And it reprojects before measuring, never after.

Common mistakes

• Assigning a CRS when you should reproject. gdf.set_crs(...) changes the label. gdf.to_crs(...) recalculates coordinates. Use set_crs only on a layer with no CRS whose correct CRS you know.
• Calculating area in EPSG:4326. The answer comes out in square degrees. Reproject first.
• Clipping layers in different CRS. clip() expects both inputs in the same CRS and warns when they differ. Reproject the boundary too, as the script does.
• Overwriting output. Writing to a layer name that already exists in a GeoPackage replaces it, and the file itself is only created on the first write. Decide whether you want a fresh file each run, and delete the old one if so.
• Silent failures. A script that skips a bad file without saying so is worse than one that crashes. Print what happened to each input.
• Invalid geometries. Clipping can fail on self-intersecting polygons. If it does, try gdf.geometry = gdf.geometry.make_valid() before the clip.

Key takeaways

• Automate the steps you repeat and the steps where a skipped check is costly.
• A GeoDataFrame is a table with a geometry column, so pandas habits carry over.
• Keep the order: check CRS, reproject, then clip and measure.
• Print a line per input so the run explains itself.
• Start small. This script is under 40 lines; add steps only when you need them.