Delivering GeoJSON and GeoPackage to BI Tools

This page solves one exact task: exporting scored trade-area polygons from a GeoDataFrame into two spatial deliverables — a WGS84, RFC 7946-compliant GeoJSON for BI tools and a multi-layer, indexed GeoPackage for desktop GIS — so Tableau, Power BI, and QGIS ingest them without a coordinate-reference surprise.

BI and GIS tools are where a scored trade area becomes an interactive map a stakeholder can filter and drill into. These two formats are the machine-readable end of the automated reporting and GIS export stage: GeoJSON is the lingua franca of web and BI ingestion, and GeoPackage is the portable container desktop GIS prefers. The task is narrow but unforgiving — GeoJSON mandates a specific coordinate reference system and polygon winding order, and getting either wrong produces a file that loads but renders inside-out or in the wrong place.

Prerequisites #

Before running this task you need:

• Python packages: geopandas (with the pyogrio or fiona engine) and shapely for the geometry, plus pyproj for CRS handling. Install with pip install geopandas shapely pyproj pyogrio.
• A scored trade-area frame. A GeoDataFrame with one polygon per candidate site, a suitability_score, attribute columns, and a CRS set. This comes from building weighted site suitability scores and the ranking and shortlisting stage.
• The parent context. This page implements the write_geojson and write_gpkg writers from Automated Reporting & GIS Export; read that first for the assembler and the coordinate-precision budget.

Configuration and execution parameters #

The two writers share a CRS policy but diverge on structure. The table below lists the parameters that determine whether a BI tool ingests the file cleanly.

The one absolute is GeoJSON’s CRS: RFC 7946 defines coordinates as WGS84 longitude/latitude, full stop. A GeoJSON in any other CRS is technically malformed, and BI tools that trust the spec will place its features in the wrong hemisphere. GeoPackage is more forgiving — it records whatever CRS you give it in its metadata — but keeping a single, documented CRS per layer avoids confusing desktop analysts.

Annotated implementation #

The writers below reproject explicitly to EPSG:4326 for GeoJSON, enforce the right-hand rule, round coordinates to the configured precision, and write a multi-layer GeoPackage with an R-tree index. Neither mutates the source frame.

from __future__ import annotations

from pathlib import Path

import geopandas as gpd
from shapely.geometry.polygon import orient
from shapely.geometry import MultiPolygon
from pyproj import CRS

GEOJSON_CRS = CRS.from_epsg(4326)      # RFC 7946 mandates WGS84


def _normalize(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Coerce every geometry to a valid, right-hand-oriented MultiPolygon."""
    def fix(geom):
        if geom is None or geom.is_empty:
            return geom
        if geom.geom_type == "Polygon":
            geom = MultiPolygon([orient(geom, sign=1.0)])   # CCW exterior
        elif geom.geom_type == "MultiPolygon":
            geom = MultiPolygon([orient(p, sign=1.0) for p in geom.geoms])
        return geom
    out = gdf.copy()
    out["geometry"] = out.geometry.apply(fix)
    return out


def write_geojson(sites: gpd.GeoDataFrame, cfg) -> Path:
    """Write an RFC 7946 GeoJSON in WGS84 with bounded coordinate precision."""
    assert sites.crs is not None, "scored sites have no CRS; refusing to export"
    # Explicit reprojection to WGS84 — never assume the source is already 4326.
    wgs = sites.to_crs(GEOJSON_CRS)
    assert CRS.from_user_input(wgs.crs) == GEOJSON_CRS
    wgs = _normalize(wgs)
    assert wgs.geometry.is_valid.all(), "invalid geometry before GeoJSON write"

    out_path = cfg.out_dir / f"trade_areas_{cfg.run_id}.geojson"
    out_path.parent.mkdir(parents=True, exist_ok=True)
    # COORDINATE_PRECISION caps decimals -> smaller, deterministic output.
    wgs.to_file(out_path, driver="GeoJSON",
                COORDINATE_PRECISION=cfg.coord_precision, RFC7946="YES")
    return out_path


def write_gpkg(sites: gpd.GeoDataFrame, cfg) -> Path:
    """Write a multi-layer, R-tree-indexed GeoPackage for desktop GIS."""
    assert sites.crs is not None, "scored sites have no CRS; refusing to export"
    sites = _normalize(sites)
    assert sites.geometry.is_valid.all(), "invalid geometry before GeoPackage write"

    out_path = cfg.out_dir / f"site_selection_{cfg.run_id}.gpkg"
    out_path.parent.mkdir(parents=True, exist_ok=True)

    # Point layer (centroids) and polygon layer (trade areas), same container.
    centroids = sites.copy()
    centroids["geometry"] = sites.geometry.centroid
    centroids.to_file(out_path, layer="sites", driver="GPKG",
                      SPATIAL_INDEX="YES")
    sites.to_file(out_path, layer="trade_areas", driver="GPKG",
                  SPATIAL_INDEX="YES")
    return out_path

write_geojson reprojects to WGS84 unconditionally and asserts the result — it never trusts that the source frame is already in EPSG:4326, because a scored frame often arrives in the equal-area CRS used for the ranking math. The RFC7946="YES" driver option and the explicit orient(..., sign=1.0) together guarantee counter-clockwise exterior rings. COORDINATE_PRECISION caps decimals at write time, which is both a size optimization and a determinism guarantee. write_gpkg writes two layers into one SQLite container with SPATIAL_INDEX="YES", so GDAL builds the R-tree that makes QGIS pan and query responsive.

Both formats sit alongside the columnar analytical copy of the data — for repeated querying at scale, scored trade areas are also written as partitioned GeoParquet to object storage, covered in configuring AWS S3 for geospatial data lakes.

Failure modes and debugging #

The winding-order bug is the sneaky one. Pre-RFC 7946 GeoJSON did not constrain ring orientation, so many older files have clockwise exteriors that most viewers tolerate — until a strict RFC 7946 consumer (some Power BI and Mapbox paths) interprets a clockwise exterior as a hole and renders the polygon inside-out. Enforcing counter-clockwise exteriors at write time removes the ambiguity.

Verification #

Round-trip every file — read it back and check what a consumer will actually see:

1. CRS round-trip: the re-read GeoJSON reports EPSG:4326; the GeoPackage reports its intended CRS.
2. Geometry validity: every re-read geometry passes is_valid.
3. Feature and layer counts: the GeoJSON feature count matches the source; the GeoPackage exposes both expected layers.

import geopandas as gpd
import pyogrio
from pyproj import CRS

def verify_exports(geojson_path, gpkg_path, expected_n: int) -> None:
    gj = gpd.read_file(geojson_path)
    assert CRS.from_user_input(gj.crs) == CRS.from_epsg(4326), "GeoJSON not WGS84"
    assert len(gj) == expected_n, f"GeoJSON {len(gj)} features, expected {expected_n}"
    assert gj.geometry.is_valid.all(), "invalid geometry in GeoJSON round-trip"

    layers = {name for name, _ in pyogrio.list_layers(gpkg_path)}
    assert {"sites", "trade_areas"} <= layers, f"missing GPKG layers: {layers}"
    areas = gpd.read_file(gpkg_path, layer="trade_areas")
    assert len(areas) == expected_n, f"GPKG {len(areas)} rows, expected {expected_n}"
    assert areas.geometry.is_valid.all(), "invalid geometry in GeoPackage round-trip"

Re-reading and re-asserting the CRS is the check that catches the most damaging silent failure — a projected frame written where WGS84 was promised. A file that loads is not a file that is correct; only the round-trip confirms a BI tool will place the features where they belong.

Related #

• Automated Reporting & GIS Export — the assembler and precision budget these writers implement.
• Configuring AWS S3 for Geospatial Data Lakes — partitioned GeoParquet for the analytical copy of these polygons.
• Ranking and Shortlisting Candidate Sites — supplies the scored polygons these files carry.

← Back to Automated Reporting & GIS Export