CAD/GIS/BIM Interop

Mapping GML Geometry to PostGIS with OGR

To load GML or CityGML into PostGIS from Python, open the source with GDAL’s GML driver, create the target layer with an explicitly assigned coordinate reference system, copy features inside a single transaction, and build the spatial index afterwards. The driver does the parsing; what it will not do is guarantee that the schema it inferred and the SRID it resolved are the ones you wanted. This page belongs to CityGML and GML Interchange.

How the GML Driver Builds a Schema

GML has no fixed schema — a document declares its own feature types — so OGR infers one by scanning the file and recording the feature classes it finds, their attributes and their geometry types. The result is cached as a .gfs file next to the source.

Where the inferred schema comes from A call sequence. The driver scans the document to discover feature classes and their attributes, caches the result as a schema file next to the source, and reuses that cache on subsequent opens. A schema inferred from a partial scan therefore persists — including any attribute the scan never reached — until the cache is deleted or replaced. pipeline GML driver the .gfs cache Open(gml) cache present? no — scan the document write inferred schema layer with fields

Two consequences follow. First, the inference is only as complete as the scan: an attribute that appears for the first time in a late feature may be missing from the schema, and the values with it. Second, the cached .gfs is reused on subsequent opens, so a schema inferred once from a partial scan persists until the file is deleted. Supplying a .gfs under your own control turns both from hazards into configuration.

The coordinate reference system is inferred separately, from srsName. Where that is a URN rather than a simple authority code, resolution varies between GDAL builds, and an unresolved system produces a layer with SRID 0 — which loads successfully and then matches nothing in a spatial join.

Production-Ready Script

# GDAL>=3.6 (osgeo), Python 3.9+
from __future__ import annotations

from dataclasses import dataclass
from osgeo import gdal, ogr, osr

gdal.UseExceptions()
ogr.UseExceptions()


@dataclass(frozen=True)
class LoadReport:
    source_features: int
    written: int
    skipped_null_geometry: int
    srid: int

    def yield_rate(self) -> float:
        return self.written / self.source_features if self.source_features else 0.0


def load_gml_to_postgis(
    gml_path: str, pg_dsn: str, table: str, epsg: int, *, layer_index: int = 0,
) -> LoadReport:
    """Copy one GML layer into PostGIS with an explicitly assigned SRID."""
    gdal.SetConfigOption("GML_EXPOSE_GML_ID", "YES")   # keep the source identifier
    gdal.SetConfigOption("GML_READ_MODE", "SEQUENTIAL_LAYERS")

    src = ogr.Open(gml_path)
    if src is None:
        raise RuntimeError(f"cannot open {gml_path}")
    src_layer = src.GetLayer(layer_index)
    source_features = src_layer.GetFeatureCount()

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)                           # explicit, not inferred
    srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)

    dst = ogr.Open(pg_dsn, update=1)
    if dst is None:
        raise RuntimeError("cannot open the PostGIS connection")

    dst_layer = dst.CreateLayer(
        table, srs=srs, geom_type=src_layer.GetGeomType(),
        options=["OVERWRITE=YES", "GEOMETRY_NAME=geom", "SPATIAL_INDEX=NONE"],
    )
    defn = src_layer.GetLayerDefn()
    for i in range(defn.GetFieldCount()):
        dst_layer.CreateField(defn.GetFieldDefn(i))

    written = skipped = 0
    dst_layer.StartTransaction()
    try:
        for feature in src_layer:
            geom = feature.GetGeometryRef()
            if geom is None or geom.IsEmpty():
                skipped += 1
                continue
            geom.AssignSpatialReference(srs)
            out = ogr.Feature(dst_layer.GetLayerDefn())
            out.SetFrom(feature)
            out.SetGeometry(geom)
            dst_layer.CreateFeature(out)
            written += 1
        dst_layer.CommitTransaction()
    except Exception:
        dst_layer.RollbackTransaction()
        raise

    dst.ExecuteSQL(f'CREATE INDEX ON "{table}" USING GIST (geom)')
    return LoadReport(source_features, written, skipped, epsg)


if __name__ == "__main__":
    report = load_gml_to_postgis(
        "city.gml", "PG:dbname=city user=loader", "buildings_lod1", 25832)
    print(report, f"yield {report.yield_rate():.1%}")
The load order that keeps the index cheap Four stages. The target layer is created with an explicitly assigned coordinate reference system and no spatial index; features are copied inside one transaction; the transaction commits; the index is built once over the loaded rows. Creating the index first makes every insert pay index maintenance, which on a bulk load costs far more than the single build. Create layer explicit SRID, no index 1 SPATIAL_INDEX=NONE Copy features inside a transaction 2 null geometry counted, not written Commit or roll back cleanly 3 no partial table survives Build the GiST index once, over the rows 4 not maintained per insert

Key implementation notes:

  • SPATIAL_INDEX=NONE on creation and an explicit CREATE INDEX afterwards. Maintaining a GiST index across a bulk insert costs far more than building it once.
  • GML_EXPOSE_GML_ID keeps the source identifier as an attribute. Without it the load produces rows that cannot be matched back to the source document.
  • The SRID is assigned from an EPSG code you supply and pushed onto every geometry, rather than being whatever the driver resolved.
  • SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) fixes GDAL 3’s axis handling to easting-northing, the equivalent of always_xy in pyproj.
  • Null and empty geometries are counted separately from failures, because they mean different things and both should be visible in the report.

Compatibility Matrix

Component Supported range Notes
GDAL / OGR >=3.6 axis mapping strategy required from GDAL 3
PostGIS 3.x GiST index, GEOMETRY and typed columns
GML driver GML, CityGML CityGML read through the GML driver
Schema inference .gfs cache delete or supply deliberately
Geometry types 2D and 3D declare a Z geometry type where elevation matters

Fallback Strategies

1. SRID comes out as 0. The driver did not resolve srsName. Assign it explicitly, as above, and assert the result with Find_SRID after the load rather than assuming the assignment took.

How an SRID of zero happens and what it costs Four stages of a silent failure. The document names its coordinate system in a URN form; the driver cannot resolve that form to an authority code; the layer is created with an SRID of zero; and every spatial join against a properly referenced layer then returns no matches. The load succeeds at every stage and the data is unusable. URN srsName in the document SRID 0 on the layer No matches nothing raised unresolved join Assign the SRID explicitly and assert it after the load rather than reading it back hopefully.

2. A missing attribute column. The .gfs was inferred from a partial scan. Delete the cached file, force a full scan, or — better — commit a .gfs you control alongside the source so the schema is not file-order dependent.

3. Mixed geometry types in one layer. CityGML feature classes routinely mix surfaces and solids. A strictly typed PostGIS column then rejects half the load. Declare a generic geometry column deliberately, or split the load by geometry type.

4. The load succeeds and returns nothing useful. Nested features — a building with boundary surfaces as children — flatten in a way the driver chooses, not the way your model needs. This is the case for parsing it yourself; see the sibling guide on parsing CityGML with lxml and Shapely.

5. Memory growth on a large file. SEQUENTIAL_LAYERS read mode keeps the driver from holding the whole document, but a single enormous feature still has to fit. Tile the source before loading if one feature does not.

FAQ

What is a .gfs file and do I need one?

It is the schema OGR infers from a GML file — the feature classes, their attributes and geometry types — cached next to the source. It matters because inference is based on a scan, and a scan that stops early can miss an attribute that appears only in later features. Supplying a .gfs you control makes the resulting table deterministic instead of dependent on file order.

Why does the loaded SRID come out as 0?

Because OGR could not resolve the srsName in the file to an authority code. URN-style names in particular are resolved inconsistently across GDAL versions. Assign the CRS explicitly on the output layer rather than letting it be inferred, and assert the SRID after the load.

Should I load with OGR or parse and insert myself?

Use OGR for bulk loading a well-formed file — it is faster and handles the driver detail. Parse it yourself when the mapping is not one-to-one: when nested features have to be flattened, when attributes need reshaping, or when features must be filtered on something the driver cannot express. The two coexist; loading a raw staging table with OGR and reshaping in SQL is often the shortest path.