Heritage GIS

Automating Flight Grid Generation with Python

Drawing a survey grid by hand in a ground-station app is slow, unrepeatable, and impossible to audit — the moment a site is re-flown next season, the lines land somewhere slightly different and the multitemporal comparison inherits the drift. Automating the grid turns the flight plan into code: a site polygon in, a set of ordered waypoints out, identical every time. This guide, part of Drone Flight Planning for Photogrammetric Survey, builds a serpentine (boustrophedon) waypoint generator with shapely and geopandas, reprojects it to EPSG:27700, and exports to CSV, KML, and GeoJSON so any ground station can fly it.

The scenario is routine: a heritage unit holds a site boundary as a GeoJSON or shapefile, and needs a flight grid that covers it at a spacing already fixed by the image overlap and GSD calculation. The generator must respect that spacing, minimise turns, and stay inside a safety buffer of the boundary.

Serpentine grid over a site polygon Parallel flight lines clipped to a buffered site polygon, connected in a boustrophedon pattern with alternate lines reversed. Serpentine grid over site polygon Alternate lines reversed; terracotta links are the connecting turns

Context & When to Use

Automate the grid when the survey must be reproducible, when the boundary is irregular enough that hand-drawing wastes coverage, or when a campaign flies many sites and manual planning does not scale. The trade-off against a ground-station wizard is upfront effort for downstream determinism: the wizard is faster for a one-off, but it cannot version-control a plan, cannot guarantee identical lines on a re-fly, and cannot be gated by the automated QA checks described in the section overview. The heading of the grid matters too — align flight lines with the long axis of the site to minimise turns, or across the dominant slope so terrain-following steps stay gentle. This generator lets you set that heading explicitly rather than accepting a default. Feed the resulting waypoints straight into corridor clipping from planning safe flight corridors over heritage sites before export.

Implementation

The generator rotates the polygon so the desired heading becomes horizontal, lays parallel lines at the planned spacing, clips each to the polygon, reverses alternate lines for the serpentine order, rotates back, and reprojects to the site grid.

# requirements.txt
# shapely==2.0.4
# geopandas==1.0.1
# pyproj==3.6.1
# simplekml==1.3.6
import math
import geopandas as gpd
import simplekml
from shapely.geometry import LineString, Point
from shapely.affinity import rotate
from shapely.ops import unary_union


def serpentine_grid(site_gdf: gpd.GeoDataFrame,
                    line_spacing_m: float,
                    heading_deg: float = 0.0,
                    edge_buffer_m: float = 5.0,
                    target_epsg: int = 27700) -> gpd.GeoDataFrame:
    """
    Build an ordered serpentine waypoint grid over a site polygon.

    site_gdf     : single-polygon GeoDataFrame in any CRS (reprojected here)
    line_spacing_m : distance between adjacent flight lines (from GSD planning)
    heading_deg  : flight-line bearing; 0 = grid lines run west-east
    edge_buffer_m: inset from the boundary to keep the aircraft clear of edges
    """
    # Work in the projected site CRS so spacing is in true metres.
    site = site_gdf.to_crs(epsg=target_epsg)          # substitute your site's EPSG
    poly = unary_union(site.geometry).buffer(-edge_buffer_m)
    if poly.is_empty:
        raise ValueError("edge_buffer_m larger than the site; nothing to fly.")

    # Rotate so the requested heading becomes horizontal (lines run in x).
    cx, cy = poly.centroid.x, poly.centroid.y
    work = rotate(poly, -heading_deg, origin=(cx, cy), use_radians=False)
    minx, miny, maxx, maxy = work.bounds

    # Lay horizontal scan lines at the planned spacing.
    waypoints, seq, line_idx = [], 0, 0
    y = miny + line_spacing_m / 2.0
    while y <= maxy:
        scan = LineString([(minx - 1, y), (maxx + 1, y)])
        clipped = scan.intersection(work)
        if not clipped.is_empty:
            segments = list(clipped.geoms) if clipped.geom_type == "MultiLineString" else [clipped]
            for seg in segments:
                pts = list(seg.coords)
                if line_idx % 2 == 1:          # reverse alternate lines
                    pts = pts[::-1]
                for x, yy in pts:
                    # Rotate the waypoint back to true site orientation.
                    p = rotate(Point(x, yy), heading_deg, origin=(cx, cy))
                    waypoints.append({"seq": seq, "line": line_idx,
                                      "geometry": p})
                    seq += 1
            line_idx += 1
        y += line_spacing_m

    if not waypoints:
        raise ValueError("No waypoints generated; check spacing vs site size.")
    return gpd.GeoDataFrame(waypoints, crs=f"EPSG:{target_epsg}")


def export_waypoints(wp: gpd.GeoDataFrame, stem: str):
    """Write CSV (lon/lat), GeoJSON (WGS84), and KML for ground stations."""
    wgs = wp.to_crs(epsg=4326)
    wgs["lon"], wgs["lat"] = wgs.geometry.x, wgs.geometry.y

    # CSV for Litchi-style importers: ordered lon/lat with sequence.
    wgs[["seq", "line", "lat", "lon"]].to_csv(f"{stem}.csv", index=False)

    # GeoJSON in WGS84 for QGIS / web review.
    wgs.to_file(f"{stem}.geojson", driver="GeoJSON")

    # KML linestring so the route is reviewable in Google Earth.
    kml = simplekml.Kml()
    coords = list(zip(wgs["lon"], wgs["lat"]))
    kml.newlinestring(name=stem, coords=coords)
    for _, r in wgs.iterrows():
        kml.newpoint(name=str(int(r["seq"])), coords=[(r["lon"], r["lat"])])
    kml.save(f"{stem}.kml")


if __name__ == "__main__":
    site = gpd.read_file("site_boundary.geojson")
    grid = serpentine_grid(site, line_spacing_m=13.1, heading_deg=25.0)
    export_waypoints(grid, "hillfort_survey_grid")
    print(f"{len(grid)} waypoints across {grid['line'].nunique()} lines")

The domain logic is in three choices: the negative edge_buffer_m insets the flight area so the aircraft never turns beyond the site edge; the heading rotation lets you align lines with the long axis of an enclosure or across a slope; and reversing alternate lines produces the continuous serpentine path that spares battery on transits.

Verification

Confirm the grid covers the site and that spacing is honoured. Load the output back and check the line count against the site width divided by spacing, and that every waypoint sits inside the buffered polygon:

grid = serpentine_grid(gpd.read_file("site_boundary.geojson"),
                       line_spacing_m=13.1, heading_deg=25.0)

# Every waypoint must lie within the site (in the site CRS).
site = gpd.read_file("site_boundary.geojson").to_crs(epsg=27700)
poly = site.union_all() if hasattr(site, "union_all") else site.unary_union
assert grid.geometry.within(poly.buffer(0.5)).all(), "waypoint outside site"

# Sequence must be gap-free and monotonic.
seqs = grid["seq"].tolist()
assert seqs == list(range(len(seqs))), "waypoint sequence broken"
print("grid valid:", len(grid), "waypoints")

Expected output is grid valid: followed by the waypoint count, with no assertion failures. Open the .kml in Google Earth as a visual check — the route should sweep back and forth with no line straying outside the boundary.

Field note. Ground stations import waypoint order literally, so a grid that looks correct in QGIS can still fly as a chaotic zig-zag if the sequence column is dropped on export. Always carry an explicit `seq` field and sort on it immediately before writing — never rely on the row order of a GeoDataFrame surviving a reprojection or a `to_file` round-trip.

Common Errors & Fixes

  • ValueError: edge_buffer_m larger than the site; nothing to fly — the negative buffer erased a small or thin site. Reduce edge_buffer_m, or confirm the boundary is in a projected CRS (a buffer of 5 “degrees” on an unprojected polygon collapses everything).
  • AttributeError: 'GeoDataFrame' object has no attribute 'union_all' — older geopandas. Use unary_union (as in the generator) or upgrade to geopandas==1.0.1; the verification snippet guards for both.
  • Waypoints appear rotated 90° from the intended heading — heading_deg follows the mathematical convention where 0 runs the lines west-east. Add or subtract 90 to swap between along-axis and across-axis sweeps, and confirm against the KML before flying.

Part of Photogrammetry & 3D Site Mapping Pipelines.