Geospatial Data Formats for Reporting Pipelines

Q: Which vector format is fastest to read into a report pipeline?

For a full-dataset read with pyogrio, FlatGeobuf and GeoPackage are consistently faster than GeoJSON because they are binary and carry a spatial index. FlatGeobuf wins when you only need a bounding-box subset, since its packed Hilbert R-tree lets the reader seek to the relevant byte ranges instead of parsing the whole file. GeoJSON is slowest because the entire text file must be parsed before any feature is available.

Q: Should I use fiona or pyogrio to read spatial data for reports?

Use pyogrio as the default engine in geopandas 0.14+ for batch report generation: it reads whole layers into a GeoDataFrame far faster than fiona because it bypasses per-feature Python object creation. Keep fiona when you need record-by-record streaming with a low, constant memory footprint, since fiona exposes an iterator that never materialises the full table.

Q: Can I read a FlatGeobuf file directly over HTTP without downloading it?

Yes. FlatGeobuf is a cloud-native format: with GDAL's /vsicurl/ virtual file system, pyogrio and fiona can issue HTTP range requests that fetch only the header, spatial index, and the byte ranges covering your bounding box. This lets a report pipeline read a small extent from a multi-gigabyte remote file without transferring the whole dataset.

Every automated spatial report begins with a read: some vector file on disk or behind a URL becomes the feature set that populates maps, attribute tables, and legends. The format you choose for that input silently governs how fast the pipeline runs, how much memory it consumes, and whether you can render a bounded extent without loading a gigabyte of geometry you will never draw. This guide compares the four formats you will actually encounter as report inputs — GeoJSON, GeoPackage (GPKG), FlatGeobuf, and the legacy Shapefile — and shows how to read each one efficiently with geopandas, fiona, and pyogrio inside a report pipeline that feeds the Dynamic Map & Data Embedding Workflows covered across this section.

Prerequisites

Before wiring format reads into a pipeline, confirm your stack meets these baselines:

Python 3.10+ with geopandas>=0.14, pyogrio>=0.7, and fiona>=1.9
GDAL 3.6+ as the backing library — both read engines depend on it, and FlatGeobuf plus /vsicurl/ range reads require a reasonably recent build
Input formats on hand: at least one .geojson, one .gpkg, and one .fgb file so you can benchmark read paths against your own data
A target report CRS decided in advance (for example a national grid such as EPSG:27700, or Web Mercator EPSG:3857 for slippy-map backdrops)
Familiarity with how a rendered feature set becomes template output — see Loop Mapping for Dynamic Attribute Tables for the downstream iteration pattern

Bash

pip install "geopandas>=0.14" "pyogrio>=0.7" "fiona>=1.9"
python -c "import pyogrio, fiona, geopandas; print(pyogrio.__gdal_version__)"

The one-liner prints the GDAL version pyogrio is linked against; if it errors or reports below 3.6, upgrade before relying on FlatGeobuf reads.

Pipeline Architecture

A report ingestion stage is a funnel: heterogeneous vector inputs enter, a read engine selected per format normalises them into a common in-memory representation, and a bounded, reprojected feature set exits toward the template layer. The diagram traces that flow and marks where format choice changes the cost.

Ingestion funnel: each format is read by the engine best suited to it, reprojected to a single report CRS, then reduced to a bounded feature context. FlatGeobuf's spatial index lets the read engine seek an extent instead of parsing the entire file.

Step-by-Step Implementation

Step 1: Inventory Input Formats and Pick a Read Engine

The first decision is which engine reads the file. pyogrio vectorises the read and returns a GeoDataFrame in a single GDAL call, which makes it the correct default for batch report jobs that load whole layers. fiona exposes a per-feature iterator, which is what you want when a dataset is too large to hold in memory and you would rather stream it — the technique detailed in Streaming FlatGeobuf Features into Report Context.

Python

# read_dispatch.py
from __future__ import annotations

import geopandas as gpd

# pyogrio is the fast, vectorised path used for full-layer reads
def read_layer(path: str, layer: str | None = None) -> gpd.GeoDataFrame:
    """Read an entire vector layer into memory using the pyogrio engine."""
    return gpd.read_file(path, layer=layer, engine="pyogrio")

The layer argument matters only for multi-layer containers. GeoJSON, FlatGeobuf, and Shapefile hold exactly one layer per file, so layer stays None; GeoPackage can hold many, so you must name the one you want or pyogrio reads the first and silently ignores the rest.

Step 2: Read the Format with a Bounding-Box Filter

Never read more geometry than the report draws. Both engines accept a bbox argument that pushes the spatial filter down into GDAL, so indexed formats (GeoPackage, FlatGeobuf) return only intersecting features without scanning the whole table.

Python

# bbox_read.py
from __future__ import annotations

import geopandas as gpd

def read_extent(
    path: str,
    bbox: tuple[float, float, float, float],
    layer: str | None = None,
) -> gpd.GeoDataFrame:
    """Read only features intersecting (minx, miny, maxx, maxy)."""
    gdf = gpd.read_file(path, layer=layer, bbox=bbox, engine="pyogrio")
    if gdf.empty:
        # A zero-feature extent is valid, not an error — hand it downstream
        return gdf
    return gdf

For GeoJSON the bbox filter still works, but it is a convenience, not a speed-up: the parser must still read the entire text file because GeoJSON carries no index. That asymmetry is the core of the GeoPackage vs FlatGeobuf for Reporting Pipelines comparison.

Step 3: Normalise the CRS to a Report Projection

GeoJSON is WGS84 by the RFC 7946 specification; GeoPackage, FlatGeobuf, and Shapefile can each carry any CRS. A report that overlays two inputs must reproject both to one projection or the geometries will not align.

Python

# crs_normalise.py
from __future__ import annotations

import geopandas as gpd

def to_report_crs(gdf: gpd.GeoDataFrame, report_epsg: int) -> gpd.GeoDataFrame:
    """Reproject to the report CRS, failing loudly on an undefined source CRS."""
    if gdf.crs is None:
        raise ValueError(
            "Layer has no CRS; refusing to assume EPSG:4326. "
            "Set the source CRS explicitly before reprojecting."
        )
    if gdf.crs.to_epsg() == report_epsg:
        return gdf
    return gdf.to_crs(epsg=report_epsg)

The explicit None check is deliberate. Shapefiles missing their .prj sidecar and GeoPackages written by loose tooling both surface as crs is None, and silently assuming EPSG:4326 is the single most common cause of maps that render in the wrong hemisphere.

Step 4: Serialise a Bounded Context for the Template

The template never touches geometry objects. Convert the reprojected frame into plain records plus the metadata the layout needs, exactly as the downstream table loop expects.

Python

# build_context.py
from __future__ import annotations

import geopandas as gpd

def build_context(gdf: gpd.GeoDataFrame, columns: list[str]) -> dict:
    """Reduce a GeoDataFrame to a JSON-safe report context."""
    minx, miny, maxx, maxy = (
        gdf.total_bounds if not gdf.empty else (0.0, 0.0, 0.0, 0.0)
    )
    return {
        "feature_count": len(gdf),
        "has_features": not gdf.empty,
        "crs_epsg": gdf.crs.to_epsg() if gdf.crs else None,
        "bbox": {"minx": minx, "miny": miny, "maxx": maxx, "maxy": maxy},
        "features": gdf[columns].to_dict("records"),
    }

The has_features boolean is what a template uses to suppress an empty map frame — the exact contract described in Conditional Rendering for Missing Spatial Data. Passing an empty list rather than None keeps the template’s for loop safe.

Step 5: Validate Schema and Geometry Before Rendering

A format read that succeeds can still hand you unusable data: an attribute column renamed by an export tool, or a self-intersecting polygon that Matplotlib will refuse to fill. Validate before the render, not after.

Python

# validate.py
from __future__ import annotations

import geopandas as gpd

REQUIRED = {"site_id", "status", "area_ha"}

def validate(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    missing = REQUIRED - set(gdf.columns)
    if missing:
        raise KeyError(f"Input is missing required columns: {sorted(missing)}")
    # Drop null geometry rows; repair invalid rings with a zero-width buffer
    gdf = gdf[gdf.geometry.notna()].copy()
    invalid = ~gdf.geometry.is_valid
    if invalid.any():
        gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].buffer(0)
    return gdf

Production-Ready Script

The following script ties the five steps into one CLI entry point with logging, error handling, and argparse. It auto-selects the read strategy from the file extension and emits the JSON-safe context a template layer consumes.

Python

#!/usr/bin/env python3
"""
ingest_features.py

Read a vector input (GeoJSON, GeoPackage, or FlatGeobuf) into a bounded,
reprojected report context.

Usage:
    python ingest_features.py --input sites.gpkg --layer parcels \\
        --report-epsg 27700 --bbox 530000 180000 545000 195000 \\
        --columns site_id status area_ha --out context.json
"""
from __future__ import annotations

import argparse
import json
import logging
import sys
from typing import Any

import geopandas as gpd

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s %(levelname)s %(message)s",
    datefmt="%H:%M:%S",
)
logger = logging.getLogger("ingest")

REQUIRED_DEFAULT = {"site_id"}


def read_extent(
    path: str,
    bbox: tuple[float, float, float, float] | None,
    layer: str | None,
) -> gpd.GeoDataFrame:
    logger.info("Reading %s (layer=%s) with pyogrio", path, layer)
    gdf = gpd.read_file(path, layer=layer, bbox=bbox, engine="pyogrio")
    logger.info("Read %d features", len(gdf))
    return gdf


def to_report_crs(gdf: gpd.GeoDataFrame, report_epsg: int) -> gpd.GeoDataFrame:
    if gdf.crs is None:
        raise ValueError("Undefined source CRS; refusing to assume EPSG:4326.")
    if gdf.crs.to_epsg() == report_epsg:
        return gdf
    logger.info("Reprojecting %s -> EPSG:%d", gdf.crs.to_epsg(), report_epsg)
    return gdf.to_crs(epsg=report_epsg)


def validate(gdf: gpd.GeoDataFrame, required: set[str]) -> gpd.GeoDataFrame:
    missing = required - set(gdf.columns)
    if missing:
        raise KeyError(f"Missing required columns: {sorted(missing)}")
    gdf = gdf[gdf.geometry.notna()].copy()
    invalid = ~gdf.geometry.is_valid
    if invalid.any():
        logger.warning("Repairing %d invalid geometries", int(invalid.sum()))
        gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].buffer(0)
    return gdf


def build_context(gdf: gpd.GeoDataFrame, columns: list[str]) -> dict[str, Any]:
    minx, miny, maxx, maxy = (
        gdf.total_bounds if not gdf.empty else (0.0, 0.0, 0.0, 0.0)
    )
    return {
        "feature_count": len(gdf),
        "has_features": not gdf.empty,
        "crs_epsg": gdf.crs.to_epsg() if gdf.crs else None,
        "bbox": {"minx": minx, "miny": miny, "maxx": maxx, "maxy": maxy},
        "features": gdf[columns].to_dict("records"),
    }


def main(argv: list[str] | None = None) -> int:
    parser = argparse.ArgumentParser(description="Ingest vector input to report context")
    parser.add_argument("--input", required=True, help="Path to .geojson/.gpkg/.fgb")
    parser.add_argument("--layer", default=None, help="Layer name (GeoPackage only)")
    parser.add_argument("--report-epsg", type=int, default=4326)
    parser.add_argument("--bbox", type=float, nargs=4, default=None,
                        metavar=("MINX", "MINY", "MAXX", "MAXY"))
    parser.add_argument("--columns", nargs="+", required=True)
    parser.add_argument("--out", default="context.json")
    args = parser.parse_args(argv)

    try:
        bbox = tuple(args.bbox) if args.bbox else None
        gdf = read_extent(args.input, bbox, args.layer)
        gdf = to_report_crs(gdf, args.report_epsg)
        gdf = validate(gdf, REQUIRED_DEFAULT | set(args.columns))
        context = build_context(gdf, args.columns)
    except (ValueError, KeyError) as exc:
        logger.error("Ingestion failed: %s", exc)
        return 2
    except Exception as exc:  # noqa: BLE001 — surface GDAL/driver errors clearly
        logger.exception("Unexpected read error: %s", exc)
        return 1

    with open(args.out, "w", encoding="utf-8") as fh:
        json.dump(context, fh, indent=2, default=str)
    logger.info("Wrote %s (%d features)", args.out, context["feature_count"])
    return 0


if __name__ == "__main__":
    sys.exit(main())

Edge Cases & Advanced Configuration

Null and Absent Data

An extent query that returns zero features is normal, not exceptional — a flood-risk layer simply has no parcels in a dry district. Propagate has_features: False downstream and let the template render an empty-state panel rather than raising. Null attribute values inside otherwise-valid features are separate: coerce them to a sentinel ("—") during context building so the table loop never emits the literal string None.

Scale and Read Performance

For layers above roughly 100,000 features, the difference between formats becomes the pipeline’s dominant cost. GeoJSON forces a full-text parse every run; GeoPackage and FlatGeobuf let the bbox filter skip untouched features via their spatial index. When the whole dataset genuinely must be read, prefer pyogrio over fiona — the vectorised path avoids constructing a Python object per feature and typically reads several times faster. When the dataset is larger than memory, invert the model and stream feature-by-feature instead of materialising a frame.

Multi-Format Ingestion

Real pipelines rarely get one format. A dispatch table keyed on the file suffix keeps the reader honest and rejects anything unmapped:

Python

from pathlib import Path

READERS = {".geojson": "single", ".json": "single",
           ".fgb": "single", ".shp": "single", ".gpkg": "multi"}

def requires_layer(path: str) -> bool:
    kind = READERS.get(Path(path).suffix.lower())
    if kind is None:
        raise ValueError(f"Unsupported input format: {Path(path).suffix}")
    return kind == "multi"

Headless and Remote Reads

In containerised CI, inputs often live in object storage rather than on disk. GDAL’s /vsicurl/ and /vsis3/ handlers let pyogrio read a FlatGeobuf extent over HTTP range requests without downloading the file — a capability GeoJSON cannot match, since it has no index to seek with. Set CPL_VSIL_CURL_ALLOWED_EXTENSIONS=.fgb and pass the /vsicurl/https://... path straight to read_file.

Troubleshooting

Symptom	Likely Cause	Resolution
`read_file` returns an empty frame for a known-good `bbox`	Coordinates given in the wrong CRS for the layer	Express the `bbox` in the layer’s native CRS, or reproject the box before filtering
Map renders in the ocean off West Africa	Layer read with an assumed EPSG:4326 when it was actually a projected grid	Read the true CRS from metadata; raise on `crs is None` instead of defaulting
GeoPackage read returns the wrong table	`layer` argument omitted on a multi-layer container	Pass `layer="..."`; list layers with `pyogrio.list_layers(path)` first
FlatGeobuf remote read downloads the entire file	`/vsicurl/` not enabled or the server lacks range-request support	Verify `Accept-Ranges: bytes` on the host; prefix the path with `/vsicurl/`
`TopologyException` when Matplotlib fills polygons	Self-intersecting rings in the source geometry	Repair with `geometry.buffer(0)` during validation before rendering
Read is far slower than expected on a large file	GeoJSON input forces a full-text parse regardless of `bbox`	Convert the input to GeoPackage or FlatGeobuf so the index can be used

Frequently Asked Questions

Which vector format is fastest to read into a report pipeline?

For a full-dataset read with pyogrio, FlatGeobuf and GeoPackage are consistently faster than GeoJSON because they are binary and carry a spatial index. FlatGeobuf wins when you only need a bounding-box subset, since its packed Hilbert R-tree lets the reader seek to the relevant byte ranges instead of parsing the whole file. GeoJSON is slowest because the entire text file must be parsed before any feature is available.

Should I use fiona or pyogrio to read spatial data for reports?

Use pyogrio as the default engine in geopandas 0.14+ for batch report generation: it reads whole layers into a GeoDataFrame far faster than fiona because it bypasses per-feature Python object creation. Keep fiona when you need record-by-record streaming with a low, constant memory footprint, since fiona exposes an iterator that never materialises the full table.

How do I keep coordinate reference systems consistent across mixed inputs?

GeoJSON is defined as WGS84 (EPSG:4326) by RFC 7946, while GeoPackage and FlatGeobuf can store any CRS. Read the source CRS from the layer metadata, then call to_crs() to reproject every input to a single report projection before merging. Never assume 4326 for GeoPackage or Shapefile inputs, and always fail loudly when a layer reports an undefined CRS.

Can I read a FlatGeobuf file directly over HTTP without downloading it?

Yes. FlatGeobuf is a cloud-native format: with GDAL’s /vsicurl/ virtual file system, pyogrio and fiona can issue HTTP range requests that fetch only the header, spatial index, and the byte ranges covering your bounding box. This lets a report pipeline read a small extent from a multi-gigabyte remote file without transferring the whole dataset.

Detailed Guides in This Section

GeoPackage vs FlatGeobuf for Reporting Pipelines — a decision guide weighing GeoPackage’s random-access, multi-layer, editable model against FlatGeobuf’s streamable, cloud-native, append-only design.
Streaming FlatGeobuf Features into Report Context — read FlatGeobuf incrementally with bounding-box and HTTP range filters to build a bounded template context without loading the whole dataset.
Choosing GeoJSON vs GeoPackage for Report Inputs — pick between human-readable text interchange and indexed binary storage based on dataset size, editing, and CRS needs.

Automated Static Map Generation from GeoJSON — the downstream stage that turns the feature set you read here into a publication-ready static map
Table Pagination Strategies for Large Attribute Tables — paginate the attribute records emitted by the context builder across PDF pages
Loop Mapping for Dynamic Attribute Tables — iterate the serialised feature records into template rows
Conditional Rendering for Missing Spatial Data — branch the template on the has_features flag when an extent returns no features

Parent: Dynamic Map & Data Embedding Workflows

Conclusion

Format choice is not a storage detail — it is the first performance decision in a report pipeline, and it decides whether reading an extent costs a byte-range seek or a full-file parse. Standardise on pyogrio for full reads, reach for fiona streaming when data outgrows memory, reproject every input to one report CRS, and prefer indexed binary formats (GeoPackage, FlatGeobuf) over GeoJSON whenever a dataset is large enough for the difference to matter.

Explore this section