Exporting Contours to GeoPackage

TL;DR: Load raw contours into GeoPandas, drop closed loops shorter than a threshold, merge line pieces with the same elevation that touch (shapely.line_merge per elevation), add elev, index and type attributes, give each line a Z coordinate equal to its elevation, and write a GeoPackage layer with the project CRS and a descriptive layer name. Export the same frame to DXF for CAD users.

# Context and Motivation

This guide is part of Contour Generation from LiDAR DTMs. Raw contour output is a technical intermediate: lines broken at every block boundary, thousands of trivial loops, floating-point elevations like 212.49999999, no distinction between index and intermediate contours, no Z coordinates. A client opening that in GIS or CAD sees clutter and wonders about quality. A delivery-ready contour layer needs a small, repeatable cleaning step, sensible attributes and self-describing metadata. GeoPackage is the natural container: a single file, open standard, CRS stored properly, readable by every GIS and by GDAL-based converters to anything else.

# Prerequisites and Assumptions

• Raw contours with an elevation attribute, from generating contours with gdal_contour or the topic’s Python route.
• GeoPandas 0.14+ with Shapely 2.x, and GDAL’s GPKG and DXF drivers.
• The contour interval, index spacing and the client’s naming conventions.

# Step-by-Step Implementation

# Step 1 — Normalize elevations

Round elev to the interval’s precision, which removes floating-point noise and makes grouping by elevation reliable.

# Step 2 — Remove short closed loops

Drop rings shorter than a threshold (for example 40 m at 0.5 m interval); they are micro-relief, not landform.

# Step 3 — Merge pieces

Group by elevation and merge touching segments with shapely.line_merge so each contour becomes as few features as possible.

# Step 4 — Attribute and add Z

Add index (every fifth level), type (“index” or “intermediate”), and set each vertex’s Z to the elevation for 3D-aware consumers.

# Step 5 — Write with metadata

Write the GeoPackage layer with a clear name, set the layer description to the interval, datum and source DTM, and optionally export DXF.

# Complete Working Example

"""Clean raw contours and write a delivery GeoPackage plus a DXF copy."""
from __future__ import annotations

from pathlib import Path

import geopandas as gpd
import numpy as np
import shapely
from shapely.geometry import LineString

INTERVAL = 0.5
INDEX_EVERY = 5
MIN_LOOP_M = 40.0


def to_3d(line: LineString, z: float) -> LineString:
    xy = np.asarray(line.coords)[:, :2]
    return LineString(np.column_stack([xy, np.full(len(xy), z)]))


def clean(raw: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    g = raw.explode(index_parts=False).reset_index(drop=True)
    g["elev"] = (np.round(g["elev"] / INTERVAL) * INTERVAL).round(3)
    rings = shapely.is_ring(g.geometry.values)
    g = g[~(rings & (g.length < MIN_LOOP_M))]
    merged = []
    for elev, grp in g.groupby("elev"):
        m = shapely.line_merge(shapely.union_all(grp.geometry.values))
        parts = list(m.geoms) if hasattr(m, "geoms") else [m]
        merged += [{"elev": elev, "geometry": p} for p in parts]
    out = gpd.GeoDataFrame(merged, crs=raw.crs)
    step = np.round(out["elev"] / INTERVAL).astype(int)
    out["index"] = (step % INDEX_EVERY == 0).astype(int)
    out["type"] = np.where(out["index"] == 1, "index", "intermediate")
    out["geometry"] = [to_3d(geom, z) for geom, z in zip(out.geometry, out["elev"])]
    return out[["elev", "index", "type", "geometry"]]


def deliver(raw_gpkg: Path, out_gpkg: Path, dxf: Path | None = None) -> gpd.GeoDataFrame:
    raw = gpd.read_file(raw_gpkg, layer="contours")
    final = clean(raw)
    layer = f"contours_{int(INTERVAL * 100):03d}cm"
    final.to_file(out_gpkg, layer=layer, driver="GPKG",
                  layer_options={"DESCRIPTION": f"{INTERVAL} m contours from 1 m LiDAR DTM; "
                                                f"elev in metres, NAVD88 (GEOID18); index every "
                                                f"{INDEX_EVERY * INTERVAL} m",
                                 "SPATIAL_INDEX": "YES"})
    if dxf is not None:
        final.rename(columns={"type": "Layer"})[["Layer", "elev", "geometry"]].to_file(dxf, driver="DXF")
    print(f"{len(raw):,} raw -> {len(final):,} delivered features in layer {layer}")
    return final


if __name__ == "__main__":
    deliver(Path("out/contours_raw.gpkg"), Path("delivery/contours.gpkg"), Path("delivery/contours.dxf"))

The DXF export uses the Layer field to put index and intermediate contours on separate CAD layers — the convention CAD users expect. The DXF driver writes only geometry and a few standard fields, which is why elevation is carried as Z.

# Key Parameter Table

# Verification

• Feature counts before and after each step, as printed, so reviewers see what cleaning removed.
• Elevations on the series and no crossings between different elevations.
• Open in two consumers. Load the GeoPackage in a GIS and the DXF in a CAD viewer; confirm CRS, Z values and layers appear as intended.

final = gpd.read_file("delivery/contours.gpkg", layer="contours_050cm")
assert final.has_z.all()
assert np.allclose((final.elev / 0.5) - np.round(final.elev / 0.5), 0)
print(final.groupby("type").size())

# Gotchas and Edge Cases

Merging too aggressively. line_merge joins only lines that share end points; it will not join pieces separated by a gap. That is correct — gaps usually mean NoData — but check that tile-edge breaks from separate contouring runs actually share vertices; if not, contour the mosaic instead.

Loop threshold on steep terrain. Real hilltops and depressions can be small closed contours. On steep ground, a 40 m threshold may remove genuine summits; lower it in hilly areas or keep loops that contain a local maximum in the DTM.

Z and 2D consumers. Some older tools reject 3D geometries. Keep a 2D version (shapely.force_2d) available if the client’s software is old.

Attribute names. Shapefile limits field names to ten characters; GeoPackage does not. If the client also wants Shapefiles, keep names short from the start.

# Frequently Asked Questions

Why deliver contours in GeoPackage?

It is a single, open, well-supported file that stores the CRS properly, supports long field names and 3D geometry, and carries a spatial index. Almost every GIS reads it, and GDAL converts it to other formats when needed.

Should contour lines have Z values?

It helps. Setting each vertex’s Z to the contour elevation lets CAD and 3D tools use the lines directly, and it costs little. Keep the elevation attribute as well for GIS styling and labelling.

How do I merge contour pieces split at tile edges?

Group lines by elevation and merge touching segments with Shapely’s line_merge after a union. This joins pieces that share end points; for pieces that do not meet, contour a mosaic instead of individual tiles.

How do I export contours for CAD?

Write the same GeoDataFrame with GDAL’s DXF driver, using a Layer field to separate index and intermediate contours and 3D line geometry to carry elevation.

# Related

• Contour Generation from LiDAR DTMs — the full workflow
• Generating Contours from a DTM with gdal_contour — producing the raw lines
• Smoothing a LiDAR DTM Before Contouring — fewer loops to remove
• Extracting Building Footprints from LiDAR — another vector deliverable from LiDAR
• Building a Tile Index with pdal tindex — GeoPackage for coverage layers