Remark

Please be aware that these lecture notes are accessible online in an ‘early access’ format. They are actively being developed, and certain sections will be further enriched to provide a comprehensive understanding of the subject matter.

3.4. KML and GML

3.4.1. Learning Objectives

By the end of this section, you should be able to:

• Describe KML’s role in visualization-oriented sharing with Google Earth and Google Maps.

• Describe GML’s role as the OGC’s general-purpose XML encoding for government and Spatial Data Infrastructure (SDI) data exchange.

• Read KML and GML files in Python, list their layers, and recognize when direct XML parsing is needed.

• Explain what is lost when a styled KML is read as data, and recover styling from the XML.

• Explain why GML coordinate order can differ from the longitude-latitude order you expect, and how srsName signals it.

• Recognize geometry-type heterogeneity within a GML layer and explain why it matters downstream.

Section Goal

KML and GML are both XML-based, but they serve different audiences: KML favors visual styling for Google Earth and Google Maps, while GML favors expressive, standards-based feature schemas for spatial data infrastructures. This section shows how to recognize and read each one, and flags a data-quality issue — mixed geometry types — that GML sources commonly introduce.

3.4.2. KML — Keyhole Markup Language

KML (Keyhole Markup Language) is an XML-based format originally developed for Google Earth and Google Maps. It is now an OGC standard [Wilson, 2008] and is widely used for publishing and sharing geographic data on the web.

3.4.2.1. Structure

KML files encode geographic features as XML elements nested within a hierarchical document structure. Features can carry stylistic information (colors, icons, line weights) and temporal data alongside their geometry.

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Placemark>
    <name>Eiffel Tower</name>
    <Point>
      <coordinates>2.2945,48.8584,0</coordinates>
    </Point>
  </Placemark>
</kml>

Table 3.21 KML Building Blocks

Element	Role
<kml>	The root element; its namespace identifies the KML version
<Document>, <Folder>	Containers that organize features into a hierarchy; tools often treat each folder as a layer
<Placemark>	One feature, with a name, a description, and a geometry
<Point>, <LineString>, <Polygon>, <MultiGeometry>	Geometry; a single folder may mix these freely
<Style> and <styleUrl>	Presentation: icon, color, and line width, defined once and referenced by features
<ExtendedData>	Custom attributes, as name-value pairs or typed fields
<TimeStamp>, <TimeSpan>	Time information, used for animation in Google Earth

Three rules about positions explain most KML surprises:

• Order is longitude,latitude,altitude, separated by commas with no spaces inside one position. Spaces separate positions in a line or ring.

• The CRS is fixed: KML always uses WGS 84 (EPSG:4326), as GeoJSON does.

• The third value is altitude in meters (the Eiffel Tower example gives 0). It is stored with the geometry, so a KML point can be a 3D point.

KML files are often delivered as KMZ, a ZIP archive containing the KML (conventionally doc.kml) plus any icons or images it references.

The first cell writes the Eiffel Tower example to disk and reads it with GeoPandas, to check what comes back.

Show code cell source

Hide code cell source

from pathlib import Path

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import pyogrio

DATA_DIR = Path("data")
DATA_DIR.mkdir(exist_ok=True)

eiffel_kml = """<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Placemark>
    <name>Eiffel Tower</name>
    <Point>
      <coordinates>2.2945,48.8584,0</coordinates>
    </Point>
  </Placemark>
</kml>"""

eiffel_path = DATA_DIR / "eiffel.kml"
eiffel_path.write_text(eiffel_kml, encoding="utf-8")

eiffel = gpd.read_file(eiffel_path)
print("CRS            :", eiffel.crs)
print("Geometry       :", eiffel.geom_type.tolist())
print("Has altitude   :", eiffel.geometry.has_z.tolist())
print("First position :", eiffel.geometry.iloc[0].coords[0])
eiffel[["Name", "geometry"]]

CRS            : EPSG:4326
Geometry       : ['Point']
Has altitude   : [True]
First position : (2.2945, 48.8584, 0.0)

	Name	geometry
0	Eiffel Tower	POINT Z (2.2945 48.8584 0)

GeoPandas reports EPSG:4326 and a three-dimensional point: the altitude 0 in the file became a z value, so the geometry prints as POINT Z. The first position is (2.2945, 48.8584), longitude first, as the rules above state.

3.4.2.2. When to Use KML

• Sharing data for visualization in Google Earth or Google Maps.

• Delivering spatial data with embedded styling information.

• Publishing field survey waypoints and tracks.

Limitations: the fixed WGS 84 CRS, a fixed vocabulary of elements (anything else must go into ExtendedData), and styling that is mixed into the file but is not part of the data model that GIS tools read.

3.4.2.3. Reading KML in Python

To see what a KML really contains, we build a small one by hand, styled the way Google Earth files are. It has a shared Style for each size class, three folders, and a mix of geometry types, and it stores a custom attribute in ExtendedData. The cities and the attribute are illustrative only.

Show code cell source

Hide code cell source

kml_text = """<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <name>Demonstration places</name>
    <Style id="small"><IconStyle><color>ff00ff00</color><scale>0.8</scale></IconStyle></Style>
    <Style id="large"><IconStyle><color>ff0000ff</color><scale>1.6</scale></IconStyle></Style>
    <Folder>
      <name>Large cities</name>
      <Placemark>
        <name>Tokyo</name><description>&lt;b&gt;Large&lt;/b&gt; city</description><styleUrl>#large</styleUrl>
        <ExtendedData><Data name="size_class"><value>large</value></Data></ExtendedData>
        <Point><coordinates>139.69,35.69,40</coordinates></Point>
      </Placemark>
      <Placemark>
        <name>Lagos</name><styleUrl>#large</styleUrl>
        <ExtendedData><Data name="size_class"><value>large</value></Data></ExtendedData>
        <Point><coordinates>3.39,6.52,0</coordinates></Point>
      </Placemark>
      <Placemark>
        <name>Mexico City</name><styleUrl>#large</styleUrl>
        <ExtendedData><Data name="size_class"><value>large</value></Data></ExtendedData>
        <Point><coordinates>-99.13,19.43,2240</coordinates></Point>
      </Placemark>
    </Folder>
    <Folder>
      <name>Small cities</name>
      <Placemark>
        <name>Reykjavik</name><styleUrl>#small</styleUrl>
        <ExtendedData><Data name="size_class"><value>small</value></Data></ExtendedData>
        <Point><coordinates>-21.94,64.15,0</coordinates></Point>
      </Placemark>
      <Placemark>
        <name>Wellington</name><styleUrl>#small</styleUrl>
        <ExtendedData><Data name="size_class"><value>small</value></Data></ExtendedData>
        <Point><coordinates>174.78,-41.29,0</coordinates></Point>
      </Placemark>
      <Placemark>
        <name>Valletta</name><styleUrl>#small</styleUrl>
        <ExtendedData><Data name="size_class"><value>small</value></Data></ExtendedData>
        <Point><coordinates>14.51,35.90,0</coordinates></Point>
      </Placemark>
    </Folder>
    <Folder>
      <name>Routes and areas</name>
      <Placemark>
        <name>Route</name>
        <LineString><coordinates>14.51,35.90 3.39,6.52</coordinates></LineString>
      </Placemark>
      <Placemark>
        <name>Study area</name>
        <Polygon><outerBoundaryIs><LinearRing>
          <coordinates>-10,0 10,0 10,15 -10,15 -10,0</coordinates>
        </LinearRing></outerBoundaryIs></Polygon>
      </Placemark>
    </Folder>
  </Document>
</kml>"""

kml_path = DATA_DIR / "demo_places.kml"
kml_path.write_text(kml_text, encoding="utf-8")

layers = pyogrio.list_layers(kml_path)
print("Layers (one per folder):", layers[:, 0].tolist())
print()
for name in layers[:, 0]:
    layer = gpd.read_file(kml_path, layer=name)
    print(f"{name:<18} {len(layer)} features, geometry {layer.geom_type.unique().tolist()}, "
          f"has altitude {layer.geometry.has_z.tolist()}")

Layers (one per folder): ['Large cities', 'Small cities', 'Routes and areas']

Large cities       3 features, geometry ['Point'], has altitude [True, True, True]
Small cities       3 features, geometry ['Point'], has altitude [True, True, True]
Routes and areas   2 features, geometry ['LineString', 'Polygon'], has altitude [False, False]

Three things to notice in the output. Folders became layers. One folder mixed geometry types (a line and a polygon together), which a Shapefile could not do. And altitude was kept: the point layers are 3D points (has_z is true), while the line and polygon, written without altitude, are 2D.

Now look at the table itself. GDAL adds a long list of KML bookkeeping columns (timestamp, altitudeMode, tessellate, drawOrder, and others) [GDAL/OGR contributors, 2026], most of them empty. The data you care about are the Name, the description (which still contains HTML), and the custom attribute from ExtendedData, which became an ordinary column.

Show code cell source

Hide code cell source

large = gpd.read_file(kml_path, layer="Large cities")
print("All columns:", large.columns.tolist())
large[["Name", "description", "size_class", "geometry"]]

All columns: ['id', 'Name', 'description', 'timestamp', 'begin', 'end', 'altitudeMode', 'tessellate', 'extrude', 'visibility', 'drawOrder', 'icon', 'size_class', 'geometry']

	Name	description	size_class	geometry
0	Tokyo	<b>Large</b> city	large	POINT Z (139.69 35.69 40)
1	Lagos	NaN	large	POINT Z (3.39 6.52 0)
2	Mexico City	NaN	large	POINT Z (-99.13 19.43 2240)

3.4.2.3.1. What Is Lost: Styling

The colors and sizes in the Style elements are not in the table. KML mixes presentation into the file, but a GIS table has no place for it. The reference to a style (styleUrl) is gone too. If you need the styling, for example to reproduce the intended look of a Google Earth map, you must read the XML directly.

KML colors have one more trap: they are written aabbggrr (alpha, blue, green, red), the reverse of the web’s #rrggbb. The cell below parses the styles and each placemark’s reference with Python’s built-in xml.etree, converts the colors, and draws the points as the author styled them over a Natural Earth [Natural Earth contributors, 2026] country backdrop.

Show code cell source

Hide code cell source

import urllib.request
import xml.etree.ElementTree as ET


def download(url, target, timeout=60):
    """Download a file once; reuse the local copy afterwards."""
    if not target.exists():
        request = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"})
        with urllib.request.urlopen(request, timeout=timeout) as response:
            target.write_bytes(response.read())
    return target


NE_BASE = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/"
countries_path = download(NE_BASE + "ne_110m_admin_0_countries.geojson", DATA_DIR / "ne_110m_admin_0_countries.geojson")
places_path = download(NE_BASE + "ne_110m_populated_places_simple.geojson", DATA_DIR / "ne_110m_populated_places_simple.geojson")
countries = gpd.read_file(countries_path)
places = gpd.read_file(places_path)

KML_NS = {"k": "http://www.opengis.net/kml/2.2"}


def kml_color_to_hex(value):
    """KML colors are aabbggrr; return the usual #rrggbb."""
    value = value.strip().lstrip("#")
    return f"#{value[6:8]}{value[4:6]}{value[2:4]}"


root = ET.parse(kml_path).getroot()
styles = {
    style.get("id"): {
        "color": kml_color_to_hex(style.findtext(".//k:color", namespaces=KML_NS)),
        "scale": float(style.findtext(".//k:scale", namespaces=KML_NS)),
    }
    for style in root.iterfind(".//k:Style", KML_NS)
}
style_of = {
    placemark.findtext("k:name", namespaces=KML_NS): (placemark.findtext("k:styleUrl", namespaces=KML_NS) or "").lstrip("#")
    for placemark in root.iterfind(".//k:Placemark", KML_NS)
}
print("Styles in the file:", styles)

cities = pd.concat([gpd.read_file(kml_path, layer=n) for n in ["Large cities", "Small cities"]], ignore_index=True)
cities["style"] = cities["Name"].map(style_of)
cities["color"] = cities["style"].map(lambda s: styles[s]["color"])
cities["size"] = cities["style"].map(lambda s: styles[s]["scale"] * 60)

Styles in the file: {'small': {'color': '#00ff00', 'scale': 0.8}, 'large': {'color': '#ff0000', 'scale': 1.6}}

Fig. 3.13 The six demonstration cities drawn with the colors and sizes recovered from the KML’s own Style elements.

In Fig. 3.13, the large cities (Tokyo, Lagos, Mexico City) are big red markers and the small ones (Reykjavik, Wellington, Valletta) are small green markers, exactly as the two styles define them. Once reordered from aabbggrr, ff0000ff is opaque red and ff00ff00 is opaque green. None of this information was in the GeoDataFrame; it came entirely from parsing the XML.

3.4.2.3.2. KMZ: The Zipped Form

A KMZ is simply a ZIP archive, and GDAL can read it directly, with no unzipping step. The cell below packs the KML into a KMZ and lists its layers.

Show code cell source

Hide code cell source

import zipfile

kmz_path = DATA_DIR / "demo_places.kmz"
with zipfile.ZipFile(kmz_path, "w", zipfile.ZIP_DEFLATED) as archive:
    archive.write(kml_path, arcname="doc.kml")

print("Layers read straight from the KMZ:", pyogrio.list_layers(kmz_path)[:, 0].tolist())
print(f"KML: {kml_path.stat().st_size:,} bytes     KMZ: {kmz_path.stat().st_size:,} bytes")

Layers read straight from the KMZ: ['Large cities', 'Small cities', 'Routes and areas']
KML: 2,450 bytes     KMZ: 715 bytes

The KMZ holds the same three layers in 715 bytes against 2,450 for the plain KML, because ZIP compression shrinks this small, repetitive XML to under a third of its size.

Note: KML files are human-readable XML, but large KML files with complex styling can be slow to parse. For analysis, convert KML to GeoJSON or GeoPackage first (see conversion workflow).

Because a folder is a layer and a layer may mix geometry types, a robust conversion has two steps: read every layer (not just the first), and keep each geometry type in its own output layer. The function below reads all layers into one table with a layer column. The cell after it writes one GeoPackage layer per folder and geometry type, dropping the empty bookkeeping columns, and exports the whole thing as a single GeoJSON.

Show code cell source

Hide code cell source

def read_all_layers(path):
    """Read every layer of a multi-layer file into one GeoDataFrame, remembering each layer's name."""
    parts = []
    for name, _ in pyogrio.list_layers(path):
        layer = gpd.read_file(path, layer=name)
        layer.insert(0, "layer", name)
        parts.append(layer)
    return pd.concat(parts, ignore_index=True)


kml_all = read_all_layers(kml_path)
display(kml_all.assign(geometry_type=kml_all.geom_type).groupby(["layer", "geometry_type"]).size().rename("features").reset_index())

gpkg_path = DATA_DIR / "demo_places.gpkg"
gpkg_path.unlink(missing_ok=True)
for (layer_name, geom_type), part in kml_all.groupby(["layer", kml_all.geom_type]):
    part.dropna(axis=1, how="all").to_file(gpkg_path, layer=f"{layer_name} ({geom_type})", driver="GPKG")
print("GeoPackage layers:", pyogrio.list_layers(gpkg_path).tolist())

geojson_path = DATA_DIR / "demo_places.geojson"
kml_all.dropna(axis=1, how="all").to_file(geojson_path, driver="GeoJSON")
print(f"GeoJSON: {geojson_path.stat().st_size:,} bytes, {len(gpd.read_file(geojson_path))} features")

	layer	geometry_type	features
0	Large cities	Point	3
1	Routes and areas	LineString	1
2	Routes and areas	Polygon	1
3	Small cities	Point	3

GeoPackage layers: [['Large cities (Point)', 'Point Z'], ['Routes and areas (LineString)', 'LineString'], ['Routes and areas (Polygon)', 'Polygon'], ['Small cities (Point)', 'Point Z']]
GeoJSON: 2,001 bytes, 8 features

The “Routes and areas” folder splits into a LineString layer and a Polygon layer, which gives four GeoPackage layers in all; the point layers keep their altitude (Point Z). The single GeoJSON holds all eight features in 2,001 bytes, so mixed geometry types are fine there, and the folder names survive only in the layer column.

3.4.2.3.3. A Real-World KML: The USGS Earthquake Feed

The earlier version of this section used a USGS earthquake feed [U.S. Geological Survey, n.d.]: the magnitude 2.5+, past-30-days KML, downloadable from the USGS Earthquake Hazards Program feeds at 2.5_month_depth.kml. Each earthquake is a styled Placemark, just like the cities above. The feed updates continuously, so the exact placemarks differ from any snapshot, and the folder (layer) names may change. That is why the cell lists the layers first and reads all of them with the function above, rather than naming one. It is optional: it is skipped if the download fails.

Show code cell source

Hide code cell source

USGS_URL = "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month_depth.kml"
quake_path = DATA_DIR / "2.5_month_depth.kml"

quakes = None
try:
    download(USGS_URL, quake_path, timeout=30)
    print("Layers in the feed:", pyogrio.list_layers(quake_path)[:, 0].tolist())
    quakes = read_all_layers(quake_path)
    print(f"{len(quakes)} placemarks")
except Exception as error:
    print(f"Could not retrieve or read the USGS feed ({type(error).__name__}); skipping this optional step.")

Layers in the feed: ['USGS Magnitude 2.5+ Earthquakes, Past Month', 'Magnitude 7', 'Magnitude 6', 'Magnitude 5', 'Magnitude 4', 'Magnitude 3', 'Magnitude 2']
2006 placemarks

Fig. 3.14 USGS earthquakes of magnitude 2.5 and above over a 30-day window, read from the KML feed (a snapshot).

The 1,899 placemarks in Fig. 3.14 trace the plate boundaries: the Aleutian arc and Alaska, the western margin of the Americas from California to Chile, the island arcs from Japan through the Philippines and Indonesia, and a looser band across the Mediterranean, Iran, and the Himalaya. Isolated clusters at Hawaii and Puerto Rico stand out against the nearly empty plate interiors. The feed is a rolling window, so a rerun gives a different snapshot with the same pattern.

3.4.3. GML — Geography Markup Language

GML (Geography Markup Language) is the OGC’s comprehensive XML encoding standard for geographic features [Portele, 2007]. It is more expressive than KML, supporting complex schemas, coordinate systems, and feature relationships. GML is commonly used in government data portals and national SDI (Spatial Data Infrastructure) systems throughout Europe and elsewhere. It is also the default response format of WFS (Web Feature Service), the OGC standard for requesting features from a server.

3.4.3.1. Characteristics

• Verbose XML syntax with namespace declarations and schema references.

• Supports complex feature schemas with type inheritance.

• Multiple GML profile levels (Simple Features, Full GML).

• Files frequently require the lxml library for robust parsing.

• A single GML file can bundle several layers; list them before loading.

Table 3.22 GML Versions and Where You Meet Them

GML version	Typical context	Notable features
2.1.2	WFS 1.0 responses	Simple; coordinates in gml:coordinates as x,y pairs; MultiPolygon
3.1.1	WFS 1.1 responses	Adds curves, surfaces, and gml:pos and gml:posList
3.2.1	WFS 2.0, INSPIRE, ISO 19136	Current standard; multi-surface geometries; stable identifiers

3.4.3.2. Structure

Unlike KML, GML does not define the features themselves: each data provider defines its own application schema (usually an XSD file), and the GML document contains features of those types. The generic parts come from the GML namespace.

Table 3.23 Main Parts of a GML Feature Collection

Part	Role
gml:featureMember	A wrapper for one feature in the collection
Application element (for example demo:cities)	The feature itself; its name is the feature type, and each child element is a property
Geometry property	Contains a GML geometry such as gml:Point, gml:Polygon, or gml:MultiPolygon (gml:MultiSurface in GML 3)
srsName attribute	Names the coordinate reference system of the geometry
gml:pos, gml:posList, gml:coordinates	The coordinate values, in the form that the GML version uses
.xsd (schema), .gfs	The application schema, and GDAL’s cached schema guess for documents without one

To work with GML without depending on an external service, the next cell builds a WFS-style GML 2 document by hand from two layers we already know: country polygons (countries) and populated places (cities), which together form one document with two feature types, just as a WFS server would return for a request naming two layers. Writing it makes the structure concrete. After writing, the cell prints one complete cities feature.

Show code cell source

Hide code cell source

from xml.sax.saxutils import escape


def coords_xml(positions):
    return " ".join(f"{x:.4f},{y:.4f}" for x, y, *_ in positions)


def ring_xml(ring):
    return f"<gml:LinearRing><gml:coordinates>{coords_xml(ring.coords)}</gml:coordinates></gml:LinearRing>"


def polygon_xml(polygon, srs=""):
    holes = "".join(f"<gml:innerBoundaryIs>{ring_xml(r)}</gml:innerBoundaryIs>" for r in polygon.interiors)
    return f"<gml:Polygon{srs}><gml:outerBoundaryIs>{ring_xml(polygon.exterior)}</gml:outerBoundaryIs>{holes}</gml:Polygon>"


def geometry_xml(geom):
    srs = ' srsName="EPSG:4326"'
    if geom.geom_type == "Point":
        return f"<gml:Point{srs}><gml:coordinates>{geom.x:.4f},{geom.y:.4f}</gml:coordinates></gml:Point>"
    if geom.geom_type == "Polygon":
        return polygon_xml(geom, srs)
    if geom.geom_type == "MultiPolygon":
        members = "".join(f"<gml:polygonMember>{polygon_xml(p)}</gml:polygonMember>" for p in geom.geoms)
        return f"<gml:MultiPolygon{srs}>{members}</gml:MultiPolygon>"
    raise ValueError(f"Unsupported geometry type: {geom.geom_type}")


def feature_type_xml(type_name, gdf):
    """One <gml:featureMember> per row; every non-geometry column becomes a property element."""
    members = []
    for i, (_, row) in enumerate(gdf.iterrows()):
        properties = "".join(f"<demo:{c}>{escape(str(row[c]))}</demo:{c}>" for c in gdf.columns if c != "geometry")
        members.append(
            f'  <gml:featureMember><demo:{type_name} fid="{type_name}.{i}">{properties}'
            f"<demo:geometry>{geometry_xml(row.geometry)}</demo:geometry></demo:{type_name}></gml:featureMember>"
        )
    return "\n".join(members)


gml_text = (
    '<?xml version="1.0" encoding="UTF-8"?>\n'
    '<wfs:FeatureCollection xmlns:demo="http://example.org/demo" '
    'xmlns:gml="http://www.opengis.net/gml" xmlns:wfs="http://www.opengis.net/wfs">\n'
    + feature_type_xml("countries", countries[["NAME", "CONTINENT", "POP_EST", "geometry"]]) + "\n"
    + feature_type_xml("cities", places[["name", "adm0name", "pop_max", "geometry"]]) + "\n"
    "</wfs:FeatureCollection>\n"
)
GML_PATH = DATA_DIR / "wfs_demo.gml"
GML_PATH.write_text(gml_text, encoding="utf-8")
print(f"{GML_PATH.name}: {GML_PATH.stat().st_size / 1024:.0f} KB of XML for {len(countries) + len(places)} features")
print()

start = gml_text.index("<gml:featureMember><demo:cities")
print(gml_text[start:gml_text.index("</gml:featureMember>", start) + len("</gml:featureMember>")])

wfs_demo.gml: 337 KB of XML for 420 features

<gml:featureMember><demo:cities fid="cities.0"><demo:name>Vatican City</demo:name><demo:adm0name>Vatican</demo:adm0name><demo:pop_max>832</demo:pop_max><demo:geometry><gml:Point srsName="EPSG:4326"><gml:coordinates>12.4534,41.9033</gml:coordinates></gml:Point></demo:geometry></demo:cities></gml:featureMember>

One cities feature shows the anatomy: an fid attribute on the feature element, one element for each property (name, adm0name, pop_max), and a geometry property wrapping a gml:Point that carries its srsName. The coordinates 12.4534,41.9033 are longitude first, the GML 2 convention. The 420 features take 337 KB of XML, which is the price of that explicitness.

3.4.3.3. Reading GML in Python

A GML document can contain more than one layer, one per feature type, so inspect it first. GDAL’s GML driver has to scan the document to discover the feature types and their properties. To avoid repeating the scan on every read, it saves what it learns in a small sidecar file with the extension .gfs [GDAL/OGR contributors, 2026]. The cell below lists the files before and after the first listing, so you can see the sidecar appear. If you later edit a GML file by hand and the structure changes, delete the .gfs so GDAL scans the document again.

Show code cell source

Hide code cell source

import fiona

print("Before listing layers:", sorted(p.name for p in DATA_DIR.glob("wfs_demo.*")))

# A GML document can contain more than one layer — inspect first
layers = fiona.listlayers(GML_PATH)         # GDAL's GML driver scans the document
print("Layers               :", layers)
print("After listing layers :", sorted(p.name for p in DATA_DIR.glob("wfs_demo.*")))
print()

# GDAL/Fiona handles GML via the GML driver
gdf = gpd.read_file(GML_PATH, layer=layers[0])
print(f"Layer '{layers[0]}': {len(gdf)} features, CRS {gdf.crs}")
gdf.drop(columns="geometry").head()

Before listing layers: ['wfs_demo.gml']
Layers               : ['countries', 'cities']
After listing layers : ['wfs_demo.gfs', 'wfs_demo.gml']

Layer 'countries': 177 features, CRS EPSG:4326

	gml_id	NAME	CONTINENT	POP_EST
0	countries.0	Fiji	Oceania	889953.0
1	countries.1	Tanzania	Africa	58005463.0
2	countries.2	W. Sahara	Africa	603253.0
3	countries.3	Canada	North America	37589262.0
4	countries.4	United States of America	North America	328239523.0

Listing the layers triggers the scan: before it only wfs_demo.gml exists, and afterwards GDAL’s .gfs sidecar sits beside it. The two feature types appear as layers, and the gml_id column exposes each feature’s fid. The population values were written as 889953.0, so GDAL inferred a floating-point type rather than an integer.

For files with complex namespaces or non-standard schemas, direct XML parsing with lxml may be necessary to extract features reliably. For very large documents, parse as a stream instead of loading the whole tree: Python’s built-in ElementTree.iterparse (and lxml’s equivalent) hands you each element as it is completed, and you can discard it afterwards. The cell below counts the features of each type in the document that way.

Show code cell source

Hide code cell source

from collections import Counter

GML_NS = "{http://www.opengis.net/gml}"
counts = Counter()
for event, element in ET.iterparse(GML_PATH, events=("end",)):
    if element.tag == GML_NS + "featureMember":
        counts[element[0].tag.split("}")[1]] += 1         # the first child names the feature type
        element.clear()                                    # release memory as we go

print("Features per type, counted by streaming the XML:", dict(counts))

Features per type, counted by streaming the XML: {'countries': 177, 'cities': 243}

The streamed counts (177 countries and 243 cities) match what GDAL read, and the whole tree was never held in memory.

3.4.3.4. Coordinate Order and srsName

GML records the CRS in an srsName attribute, and the form of that name can change the order of the coordinates. Under the rules for GML 3 [Portele, 2007], a CRS written as an EPSG URN follows the axis order defined by EPSG, and EPSG defines WGS 84 (EPSG:4326) as latitude first. This is the opposite of GeoJSON and KML.

Table 3.24 GML srsName Forms and Coordinate Order

srsName	Coordinate order
EPSG:4326	longitude, latitude (the older, GML 2 convention)
urn:ogc:def:crs:EPSG::4326	latitude, longitude
http://www.opengis.net/def/crs/EPSG/0/4326	latitude, longitude
urn:ogc:def:crs:OGC:1.3:CRS84	longitude, latitude

GDAL knows these rules and reads such files correctly [GDAL/OGR contributors, 2026]. The risk appears when you parse the XML yourself and assume the first number is the longitude. To demonstrate, the cell below has GDAL write six places as GML 3.2 (its default form), then reads the coordinates twice: naively, with the first number as x, and with GDAL.

Show code cell source

Hide code cell source

from shapely.geometry import Point

AXIS_GML = DATA_DIR / "places_gml32.gml"
sample = places[["name", "geometry"]].head(6)
sample.to_file(AXIS_GML, driver="GML")           # GDAL writes GML 3.2 with an EPSG URN srsName

ns = {"gml": "http://www.opengis.net/gml/3.2"}
tree_root = ET.parse(AXIS_GML).getroot()
first_point = tree_root.find(".//gml:Point", ns)
print("srsName :", first_point.get("srsName"))
print("gml:pos :", first_point.findtext("gml:pos", namespaces=ns), "  <- first number is the LATITUDE")

naive = [
    Point(*map(float, point.findtext("gml:pos", namespaces=ns).split()))      # assumes x = first number
    for point in tree_root.iterfind(".//gml:Point", ns)
]
naive = gpd.GeoSeries(naive, crs="EPSG:4326")
correct = gpd.read_file(AXIS_GML)

print()
print(f"Vatican City, naive parse: ({naive.iloc[0].x:.2f}, {naive.iloc[0].y:.2f})")
print(f"Vatican City, via GDAL   : ({correct.geometry.iloc[0].x:.2f}, {correct.geometry.iloc[0].y:.2f})")
print(f"Naive points with an impossible latitude (beyond +/-90): {int((naive.y.abs() > 90).sum())} of {len(naive)}")

srsName : urn:ogc:def:crs:EPSG::4326
gml:pos : 41.903282 12.453387   <- first number is the LATITUDE

Vatican City, naive parse: (41.90, 12.45)
Vatican City, via GDAL   : (12.45, 41.90)
Naive points with an impossible latitude (beyond +/-90): 1 of 6

Fig. 3.15 The same six GML points read by GDAL (green) and by a naive parse that takes the first number as longitude (red).

In Fig. 3.15, the green points are where the places really are: four small European states, Eswatini, and Micronesia. The red points are the same numbers read in the wrong order. The European ones land in East Africa and the Arabian Sea, Eswatini’s lands in the Atlantic, and Micronesia’s falls off the map entirely because its longitude of about 158° becomes an impossible latitude. A wrong order produces plausible-looking coordinates for most places, and only an out-of-range value gives it away.

3.4.3.5. Watch for Mixed Geometry Types

GML sources frequently mix Polygon and MultiPolygon geometries within what is conceptually a single layer — for example, a district boundary made of one contiguous polygon next to a district split across two disjoint parcels. This is legal GML but becomes a problem the moment you try to load the layer into a database table that expects one geometry type per column.

Our countries layer shows this directly. Countries made of separate landmasses or islands are encoded as MultiPolygon, and single contiguous ones as Polygon. GDAL cannot declare one geometry type for such a layer, so it reports it as Unknown. The cell counts the types and maps them.

Show code cell source

Hide code cell source

countries_gml = gpd.read_file(GML_PATH, layer="countries")

print("Geometry type GDAL reports for the layer:", dict(pyogrio.list_layers(GML_PATH))["countries"])
print(countries_gml.geom_type.value_counts())

Geometry type GDAL reports for the layer: Unknown
Polygon         148
MultiPolygon     29
Name: count, dtype: int64

Fig. 3.16 Country boundaries in the GML layer, colored by geometry type.

Fig. 3.16 shows that geometry type follows geography. The 29 MultiPolygon countries (green) are those with islands or separate parts, such as the United States, Canada, Russia, Australia, Indonesia, and Angola with its Cabinda exclave, while the 148 plain Polygon countries (orange) are mostly single landmasses, including Greenland and Madagascar, which are one polygon each at this resolution. Neither type is an error, which is why GDAL can only report Unknown for the layer.

Detecting this early — right after loading, before any spatial join or database load — avoids failures much later in a pipeline. The Spatial Database Systems section shows how to normalize this heterogeneity before loading the result into PostGIS. The core idea is to convert toward the more general type: wrapping each Polygon in a one-part MultiPolygon loses nothing, while the reverse would discard parts. The cell shows the idea in three lines.

Show code cell source

Hide code cell source

from shapely.geometry import MultiPolygon


def to_multipolygon(geom):
    return MultiPolygon([geom]) if geom.geom_type == "Polygon" else geom


uniform = countries_gml.set_geometry(countries_gml.geometry.apply(to_multipolygon))
print("Before:", countries_gml.geom_type.value_counts().to_dict())
print("After :", uniform.geom_type.value_counts().to_dict())
print("Total area unchanged:", abs(uniform.to_crs(6933).area.sum() - countries_gml.to_crs(6933).area.sum()) < 1)

Before: {'Polygon': 148, 'MultiPolygon': 29}
After : {'MultiPolygon': 177}
Total area unchanged: True

Wrapping each Polygon in a one-part MultiPolygon makes all 177 countries uniform, and the total area is unchanged (compared in an equal-area CRS, EPSG:6933, to within a square meter). Nothing was lost, which is not true of the reverse conversion.

3.4.3.6. Requesting GML from a WFS

Most real GML arrives as a WFS response, and a WFS request is just a URL with parameters. Everything you learned about narrow requests applies: ask for specific feature types, a bounding box, and a feature limit, and cache the raw response. The cell below builds a GetFeature request for the public MapServer demo WFS that the earlier version of this section used, asking for its two layers (world continents and world cities) in one response, and then optionally fetches it. The continents layer is a good place to check geometry types, since continents made of separate landmasses could be encoded as MultiPolygon. The fetch is skipped if the service is unreachable.

Show code cell source

Hide code cell source

from urllib.parse import urlencode

params = {"SERVICE": "WFS", "VERSION": "1.0.0", "REQUEST": "GetFeature", "TYPENAME": "continents,cities"}
WFS_URL = "https://demo.mapserver.org/cgi-bin/wfs?" + urlencode(params, safe=",")
print(WFS_URL)

wfs_path = DATA_DIR / "wfs.gml"
try:
    download(WFS_URL, wfs_path, timeout=30)
    wfs_layers = pyogrio.list_layers(wfs_path)[:, 0].tolist()
    print("Layers in the WFS response:", wfs_layers)
    if "continents" in wfs_layers:
        continents = gpd.read_file(wfs_path, layer="continents")
        print("Continent geometry types:", continents.geom_type.value_counts().to_dict())
except Exception as error:
    print(f"Could not retrieve or read the WFS response ({type(error).__name__}); skipping this optional step.")

https://demo.mapserver.org/cgi-bin/wfs?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=continents,cities
Layers in the WFS response: ['continents', 'cities']
Continent geometry types: {'Polygon': 27458}

In this response all 27,458 continents features are plain Polygon: the service publishes the layer as thousands of single-part pieces, so there is nothing to normalize here. The check costs one line, and the answer depends on the source, so look rather than assume.

3.4.4. Comparing KML and GML

Table 3.25 Comparison of KML and GML

Aspect	KML	GML
Main purpose	Visual sharing in Google Earth and Google Maps	Standards-based exchange for government and SDI systems
Schema	Fixed vocabulary; custom attributes in ExtendedData	Application schema defined by the provider (XSD)
CRS	Always WGS 84, longitude-latitude-altitude	Any CRS, named by srsName; coordinate order depends on its form
Styling	Built in (Style, styleUrl)	None
Layers	Folders	Feature types
Geometry types in one layer	May be mixed	May be mixed (for example, Polygon with MultiPolygon)
Typical sources	Google Earth, GPS and survey apps, public feeds	WFS servers, INSPIRE, national mapping agencies
Packaging	Single file, or KMZ (zipped)	One or more files, often with .xsd and GDAL’s .gfs
Best use	Presentation	Interchange of rich, typed data

Practice

Export one layer of your choice from a GeoPackage to GML with GDAL (driver="GML"), then open the .gml and .xsd files in a text editor. Find the srsName and the first coordinate pair. Read the file back with GeoPandas and check that the coordinates match the original. Then write a short KML by hand with two styles and three placemarks, read it with GeoPandas, and recover each placemark’s color from the XML.

3.4.5. Key Terms

Table 3.26 Key Terms for KML and GML Concepts

Term	Meaning
KML	Keyhole Markup Language: an XML format for visual sharing of geographic data.
KMZ	A ZIP archive containing a KML file and its resources.
Placemark	A single KML feature: a name, a description, and a geometry.
ExtendedData	The KML element for custom attributes.
aabbggrr	The byte order of KML color values (alpha, blue, green, red).
GML	Geography Markup Language: the OGC’s general-purpose XML encoding for geographic features.
Application schema	The feature types and properties a data provider defines for its GML (usually an XSD).
srsName	The attribute that names a GML geometry’s CRS, and signals its coordinate order.
WFS	Web Feature Service: an OGC service whose default response is GML.
.gfs	GDAL’s cached description of a GML document’s feature types.

3.4.6. Review Questions

1. In what order are the three values of a KML coordinate, and what is the CRS?

2. A KML has three folders. How many layers does GDAL see, and can a layer mix geometry types?

3. After reading a styled KML with GeoPandas, the icon colors are missing. Why, and how would you recover them?

4. A KML color is ff0000ff. What color is it?

5. Why does a GML point at <gml:pos>41.9 12.5</gml:pos> with an EPSG URN srsName lie in Italy, not in East Africa?

6. What is a .gfs file, and when should you delete it?

7. A GML layer contains 118 Polygon and 7 MultiPolygon features. What problem does this cause for a database load, and what is the lossless fix?

8. When would you choose KML over GML to deliver a dataset, and when the reverse?

3.4.7. Summary

KML optimizes for visual sharing in consumer mapping tools, while GML optimizes for expressive, standards-based interchange in government and SDI systems — both trade verbosity for those benefits, so convert them to GeoPackage or GeoJSON for analysis once you have the data you need. KML folders become layers, and its styling is presentation that GIS tables do not carry, so recover it from the XML when you need it. GML in particular can bundle several layers, record coordinates in latitude-first order depending on its srsName, and mix Polygon/MultiPolygon geometries for one conceptual feature type, a heterogeneity that must be detected and normalized before database loading.

The next section turns to OpenStreetMap, a global, tag-based data source rather than a single downloadable file format.