How to Explore a Spatial Dataset You Have Never Seen Before

On this page

Problem statement

A colleague sends Leeds.zip — in fact HM Land Registry's INSPIRE Index Polygons for Leeds, downloaded exactly as the Land Registry publishes them (the setup is under The data used on this page). Inside is one GML file and a licence PDF. There is no other documentation.

import geopandas as gpd

gdf = gpd.read_file("zip://Leeds.zip!Land_Registry_Cadastral_Parcels.gml")
print(gdf.head(3).drop(columns="geometry").to_string())

                                      gml_id  INSPIREID     LABEL  NATIONALCADASTRALREFERENCE                        VALIDFROM             BEGINLIFESPANVERSION
0  PREDEFINED.fid--4959b20f_186afe90467_-7e8   54154992  54154992                    54154992        2012-05-24 23:00:00+00:00 2012-05-25 11:00:28.895000+00:00
1  PREDEFINED.fid--4959b20f_186afe90467_-7e7   50184973  50184973                    50184973 2011-01-31 14:52:05.922000+00:00 2011-01-31 14:52:05.922000+00:00
2  PREDEFINED.fid--4959b20f_186afe90467_-7e6   50184219  50184219                    50184219 2011-01-31 14:44:18.526000+00:00 2011-01-31 14:44:18.526000+00:00

Three rows of a table with columns called INSPIREID, LABEL, NATIONALCADASTRALREFERENCE and a couple of dates. You still do not know what the layer covers, what coordinate system it is in, whether those three ids mean different things, whether the geometry is usable, or how many of the rows are duplicates of each other.

Jumping straight to analysis on an unknown layer is how a wrong answer gets produced confidently. Twenty minutes of systematic interrogation prevents most of it — and the questions are always the same, which means the process can be a checklist.

Quick answer

Ask seven questions, in this order, before doing anything else:

import geopandas as gpd

print(f"1. rows/cols   {gdf.shape}")
print(f"2. CRS         {gdf.crs}")
print(f"3. bounds      {gdf.total_bounds}")
print(f"4. geometry    {gdf.geom_type.value_counts().to_dict()}")
print(f"5. health      null={gdf.geometry.isna().sum()} "
      f"empty={gdf.geometry.is_empty.sum()} invalid={(~gdf.is_valid).sum()}")
print(f"6. columns     {dict(gdf.dtypes.astype(str))}")
print(f"7. nulls       {gdf.isna().sum().to_dict()}")

1. rows/cols   (293729, 7)
2. CRS         EPSG:27700
3. bounds      [411649.666 422246.489 447962.783 450596.2  ]
4. geometry    {'Polygon': 293729}
5. health      null=0 empty=0 invalid=0
6. columns     {'gml_id': 'str', 'INSPIREID': 'int32', 'LABEL': 'int32', 'NATIONALCADASTRALREFERENCE': 'int32', 'VALIDFROM': 'datetime64[ms, UTC]', 'BEGINLIFESPANVERSION': 'datetime64[ms, UTC]', 'geometry': 'geometry'}
7. nulls       {'gml_id': 0, 'INSPIREID': 0, 'LABEL': 0, 'NATIONALCADASTRALREFERENCE': 0, 'VALIDFROM': 0, 'BEGINLIFESPANVERSION': 0, 'geometry': 0}

Each line answers something you would otherwise have assumed:

Line	Tells you
bounds in the hundreds of thousands	metres, consistent with the declared British National Grid
one geometry type, all Polygon	a strict `geometry(Polygon)` column will accept every row
no null, empty or invalid geometry	this supplier validates before publishing — not something to assume next time
no null attributes	every row has an id and both dates
three id columns	whether they are three things or one is the next question

The seven questions

Twenty minutes here, or a wrong number later. The order matters: CRS before anything measured.

Layers, a truthful CRS, usable geometry, then a plot

Plot last. A map of data you have not checked is a convincing picture of nothing.

1. Find out what is in the file before opening it

A GeoPackage or GDB can hold several layers, and read_file without a layer argument silently gives you the first one.

import fiona

print(fiona.listlayers("data/data.gpkg"))

['parcels', 'wards', 'stations']

Read the schema without loading the data — instant on a 2 GB file:

with fiona.open("data/data.gpkg", layer="parcels") as src:
    print(src.crs)
    print(src.schema)
    print(len(src))
    print(src.bounds)

EPSG:27700
{'properties': {'gml_id': 'str', 'INSPIREID': 'int32', 'LABEL': 'int32', 'NATIONALCADASTRALREFERENCE': 'int32', 'VALIDFROM': 'datetime', 'BEGINLIFESPANVERSION': 'datetime'}, 'geometry': 'Polygon'}
293729
(411649.666, 422246.489, 447962.783, 450596.2)

If the file is large, this is where you decide whether to read it all:

sample = gpd.read_file("data/data.gpkg", layer="parcels", rows=1000)

2. Establish the CRS — and check it is telling the truth

print(gdf.crs)
print(gdf.crs.is_projected)
print(gdf.crs.axis_info[0].unit_name)
print(gdf.total_bounds.round(1).tolist())

EPSG:27700
True
metre
[411649.7, 422246.5, 447962.8, 450596.2]

The declared CRS and the actual coordinates can disagree, and nothing warns you. The bounds are the check:

values between -180 and 180 → degrees, so a geographic CRS
values in the hundreds of thousands → metres, a projected national grid
values in the millions → probably Web Mercator (EPSG:3857)

A layer declaring EPSG:27700 with bounds of [-3.2, 55.9, -3.1, 56.0] is mislabelled — the coordinates are degrees. Fixing that is set_crs, not to_crs, and getting it backwards puts the data in the North Sea.

if gdf.crs is None:
    raise ValueError("no CRS — do not guess, ask the supplier")

3. Check the geometry is usable

health = {
    "null":    int(gdf.geometry.isna().sum()),
    "empty":   int(gdf.geometry.is_empty.sum()),
    "invalid": int((~gdf.is_valid).sum()),
    "types":   gdf.geom_type.value_counts().to_dict(),
}
print(health)

{'null': 0, 'empty': 0, 'invalid': 0, 'types': {'Polygon': 293729}}

Four separate conditions, four separate checks — and dropna() catches exactly one of them. See null, empty, missing and invalid for why each behaves differently.

If anything is invalid, find out how:

from shapely.validation import explain_validity
bad = gdf[~gdf.is_valid]
print(bad.geometry.apply(lambda g: explain_validity(g).split("[")[0]).value_counts())

The Land Registry file has nothing to explain: every one of its polygons is valid. On a layer that has been edited by hand, a few dozen self-intersections in tens of thousands of rows is normal digitising noise; thousands means something upstream is broken.

4. Understand the attributes, not just their names

for col in gdf.columns.drop("geometry"):
    s = gdf[col]
    print(f"{col:14s} {str(s.dtype):8s} nulls={s.isna().sum():5d} "
          f"unique={s.nunique():6d}  e.g. {str(s.dropna().iloc[0])[:44] if s.notna().any() else '—'}")

gml_id         str      nulls=    0 unique=293729  e.g. PREDEFINED.fid--4959b20f_186afe90467_-7e8
INSPIREID      int32    nulls=    0 unique=293729  e.g. 54154992
LABEL          int32    nulls=    0 unique=293729  e.g. 54154992
NATIONALCADASTRALREFERENCE int32    nulls=    0 unique=293729  e.g. 54154992
VALIDFROM      datetime64[ms, UTC] nulls=    0 unique=152845  e.g. 2012-05-24 23:00:00+00:00
BEGINLIFESPANVERSION datetime64[ms, UTC] nulls=    0 unique=193409  e.g. 2012-05-25 11:00:28.895000+00:00

Three things to look for:

# a) a column that should be a key — is it?
print("INSPIREID unique:", gdf["INSPIREID"].is_unique)

# b) are several columns really one?
same = (gdf["INSPIREID"] == gdf["LABEL"]) & (gdf["INSPIREID"] == gdf["NATIONALCADASTRALREFERENCE"])
print("three id columns identical on every row:", bool(same.all()))

# c) a date column is often the most informative field in the file
print(gdf["VALIDFROM"].dt.year.value_counts().sort_index().tail(5).to_string())

# d) precomputed geometry columns are often stale — the wards layer carries one
wards = gpd.read_file("data/data.gpkg", layer="wards")
print("Shape__Area / computed area:", (wards["Shape__Area"] / wards.geometry.area).round(4).unique())

INSPIREID unique: True
three id columns identical on every row: True
VALIDFROM
2019    3493
2020    2373
2021    2089
2022    2587
2023     365
Shape__Area / computed area: [1.]

Three of the seven columns are the same number, which is worth knowing before anyone builds a join on "the cadastral reference" believing it differs from the id. The dates say when each title's extent was last changed — and the latest year in the file, 2023, with only a few hundred changes, says something the download page does not: this is an extract from early 2023, not from this year. The years are also lumpy: 66,666 titles carry a 2002 date and 52,452 a 2008 one — far more than one city's ordinary changes in a year — which looks like bulk loading or re-capture of the index. For those titles the date probably says when the polygon was drawn, not when the land changed. That is an inference from the shape of the column, and exactly the kind of question to put to the supplier. And the last check is worth running every time. Shape__Area, SHAPE_Area, area_ha and friends were computed at some point in the past, in some CRS, and are often wrong after editing — or computed in Web Mercator by the service that exported them. Here the ONS's stored areas agree with British National Grid to four decimal places, which is the answer you hope for and should not assume.

Steps 3 to 5 on the Land Registry file, side by side.

5. Look for duplicates, in both senses

# attribute duplicates
print("attribute duplicates:", gdf.drop(columns="geometry").duplicated().sum())

# geometry duplicates — the expensive kind to miss
wkb = gdf.geometry.to_wkb()
print("geometry duplicates: ", wkb.duplicated().sum())

attribute duplicates: 0
geometry duplicates:  790

Every row has its own id, and 790 of them have exactly the same polygon as another row. In a land register that is not necessarily an error — the index polygon of a title can coincide exactly with another's — but in any analysis it doubles the area of that land and produces duplicate rows in every join. Identical shapes arriving from an ordinary supplier usually mean a merge happened twice. See how to find and remove duplicate geometries.

6. Check the spatial extent makes sense

minx, miny, maxx, maxy = gdf.total_bounds
print(f"{(maxx - minx) / 1000:.1f} km × {(maxy - miny) / 1000:.1f} km")

36.3 km × 28.3 km

That is about the size of Leeds, which is what the layer is supposed to cover. A bounding box hundreds of kilometres across would mean either the layer is bigger than you were told, or a handful of features have coordinates at the origin or in the wrong hemisphere dragging the extent out:

c = gdf.geometry.centroid
outliers = gdf[(c.x < c.x.quantile(0.001)) | (c.x > c.x.quantile(0.999))]
print(f"{len(outliers)} features in the outer 0.1% on each side, spanning "
      f"{outliers.total_bounds.round().tolist()}")

588 features in the outer 0.1% on each side, spanning [411650.0, 426727.0, 447963.0, 449317.0]

The quantile test always returns something — 0.1% of 293,729 rows is hundreds — so read the extent it prints, not the count. Here the outliers lie within the district's own bounds: just its east and west edges, not a stray point at (0, 0).

A few features at (0, 0) is the classic signature of rows whose coordinates failed to parse. See how to find and remove spatial outliers.

7. Now plot it

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 7))
gdf.plot(ax=ax1, edgecolor="none")
ax1.set_title(f"all {len(gdf):,} features")

gdf.assign(changed=gdf["VALIDFROM"].dt.year).plot(ax=ax2, column="changed", legend=True, cmap="viridis")
ax2.set_title("by year the extent was last changed")
plt.tight_layout()

Plot last, not first. A map drawn before the checks is convincing regardless of whether the data is any good — and a plot of a layer with hundreds of empty geometries, or 790 stacked duplicates, looks exactly like a plot of a layer with none.

Profiling Leeds parcels, layer by layer, into assertions

The data used on this page

Leeds.zip is HM Land Registry's INSPIRE Index Polygons for Leeds, downloaded as published: every registered freehold title's extent, 293,729 polygons in one GML file. data/data.gpkg puts those parcels in a GeoPackage next to two more layers — Leeds' 33 wards from the ONS, exported with the service's own Shape__Area column, and West Yorkshire's railway station entrances from NaPTAN — to have a multi-layer file to explore. Land Registry data © Crown copyright and database right 2026, OGL v3.

import io
from pathlib import Path

import geopandas as gpd
import pandas as pd
import pyogrio
import requests

reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)                    # sets a cookie, then redirects to itself
reply.raise_for_status()
Path("Leeds.zip").write_bytes(reply.content)
pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})   # else GDAL waits minutes on the file's schema URL

Path("data").mkdir(exist_ok=True)
gpkg = "data/data.gpkg"
gpd.read_file("zip://Leeds.zip!Land_Registry_Cadastral_Parcels.gml").to_file(gpkg, layer="parcels")
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
gpd.read_file(f"{ONS}/Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query"
              "?where=LAD24NM%3D%27Leeds%27&outFields=WD24CD,WD24NM,Shape__Area&outSR=27700&f=geojson"
              ).set_crs(27700, allow_override=True).to_file(gpkg, layer="wards")
# NaPTAN area 450 as saved on 9 October 2026 (the register changes weekly); live:
# naptan.api.dft.gov.uk/v1/access-nodes?dataFormat=csv&atcoAreaCodes=450
naptan = pd.read_csv(io.StringIO(requests.get("https://data.spatialworkflow.io/naptan-area-450-2026-10-09.csv", timeout=300).text),
                     low_memory=False)
rse = naptan[(naptan["Status"] == "active") & (naptan["StopType"] == "RSE")]
gpd.GeoDataFrame(rse[["ATCOCode", "CommonName"]], crs=27700,
                 geometry=gpd.points_from_xy(rse["Easting"], rse["Northing"])).to_file(gpkg, layer="stations")
print(pyogrio.list_layers(gpkg).tolist())

[['parcels', 'Polygon'], ['wards', 'Polygon'], ['stations', 'Point']]

Example 1: a reusable profiling function

import geopandas as gpd
from shapely.validation import explain_validity

def profile(gdf: gpd.GeoDataFrame, name: str = "layer") -> dict:
    """Everything worth knowing about a layer before using it."""
    geom = gdf.geometry
    invalid = gdf[~gdf.is_valid]

    return {
        "name": name,
        "rows": len(gdf),
        "columns": {c: str(t) for c, t in gdf.dtypes.items() if c != "geometry"},
        "crs": str(gdf.crs),
        "projected": bool(gdf.crs.is_projected) if gdf.crs else None,
        "units": gdf.crs.axis_info[0].unit_name if gdf.crs else None,
        "bounds": [round(v, 1) for v in gdf.total_bounds],
        "geom_types": geom.geom_type.value_counts().to_dict(),
        "null_geom": int(geom.isna().sum()),
        "empty_geom": int(geom.is_empty.sum()),
        "invalid_geom": len(invalid),
        "invalid_reasons": (
            invalid.geometry.apply(lambda g: explain_validity(g).split("[")[0])
            .value_counts().to_dict() if len(invalid) else {}
        ),
        "duplicate_geom": int(geom.to_wkb().duplicated().sum()),
        "attribute_nulls": {c: int(n) for c, n in gdf.isna().sum().items() if n},
        "candidate_keys": [
            c for c in gdf.columns
            if c != "geometry" and gdf[c].is_unique and gdf[c].notna().all()
        ],
    }

import json

print(json.dumps(profile(gdf, "leeds_inspire"), indent=2, default=str))

{
  "name": "leeds_inspire",
  "rows": 293729,
  "columns": {
    "gml_id": "str",
    "INSPIREID": "int32",
    "LABEL": "int32",
    "NATIONALCADASTRALREFERENCE": "int32",
    "VALIDFROM": "datetime64[ms, UTC]",
    "BEGINLIFESPANVERSION": "datetime64[ms, UTC]"
  },
  "crs": "EPSG:27700",
  "projected": true,
  "units": "metre",
  "bounds": [
    411649.7,
    422246.5,
    447962.8,
    450596.2
  ],
  "geom_types": {
    "Polygon": 293729
  },
  "null_geom": 0,
  "empty_geom": 0,
  "invalid_geom": 0,
  "invalid_reasons": {},
  "duplicate_geom": 790,
  "attribute_nulls": {},
  "candidate_keys": [
    "gml_id",
    "INSPIREID",
    "LABEL",
    "NATIONALCADASTRALREFERENCE"
  ]
}

Save the output next to the data. Three months later it is the only record of what the file looked like on arrival — and the fastest way to prove that a supplier's "we didn't change anything" is wrong.

A date column describes the file’s history as well as the land’s. Two single years hold 40% of these dates.

Example 2: profiling every layer in a folder

from pathlib import Path

import pandas as pd


def profile_folder(folder: Path) -> list[dict]:
    out = []
    for path in sorted(folder.rglob("*.gpkg")):
        for layer in fiona.listlayers(path):
            try:
                gdf = gpd.read_file(path, layer=layer)
                out.append({"file": str(path), **profile(gdf, layer)})
            except Exception as exc:
                out.append({"file": str(path), "name": layer,
                            "error": f"{type(exc).__name__}: {exc}"})
    return out

rows = profile_folder(Path("data"))
summary = pd.DataFrame(rows)[["file", "name", "rows", "crs", "invalid_geom", "duplicate_geom"]]
print(summary.to_string(index=False))

          file     name   rows        crs  invalid_geom  duplicate_geom
data/data.gpkg  parcels 293729 EPSG:27700             0             790
data/data.gpkg    wards     33 EPSG:27700             0               0
data/data.gpkg stations    155 EPSG:27700             0               0

Running this before a batch job answers the question that batch jobs fail on: are all these files actually in the same CRS with the same schema? See how to build an inventory of a GIS data folder.

Example 3: turning the profile into assertions

Once you know what the layer should look like, encode it:

def assert_expected(gdf, *, crs=27700, geom_types=("Polygon", "MultiPolygon"),
                    min_rows=1000, key="INSPIREID", max_duplicate_geom=1000):
    assert gdf.crs and gdf.crs.to_epsg() == crs, f"CRS is {gdf.crs}"
    assert len(gdf) >= min_rows, f"only {len(gdf)} rows"
    unexpected = set(gdf.geom_type.dropna()) - set(geom_types)
    assert not unexpected, f"unexpected geometry types: {unexpected}"
    assert gdf[key].is_unique, f"{key} is not unique"
    assert not gdf.geometry.is_empty.any(), f"{gdf.geometry.is_empty.sum()} empty"
    dup = int(gdf.geometry.to_wkb().duplicated().sum())
    assert dup <= max_duplicate_geom, f"{dup} duplicate geometries — more than this supplier's usual"


assert_expected(gdf)
print("delivery matches expectations")

delivery matches expectations

That is the moment exploration turns into a pipeline: the things you checked by hand once become the things checked automatically on every delivery. See how to validate a GeoDataFrame against a schema.

Spatial data fails quietly

Every row is an assumption that produces a plausible wrong answer rather than an error.

The reason this checklist works is that spatial data fails quietly. A CSV with a missing column raises a KeyError the first time you use it. A shapefile with the wrong CRS produces valid geometry, plots without complaint, joins to other layers returning zero rows, and reports areas that are wrong by a factor of ten billion — all without a single exception.

Each of the seven questions targets one of those silent failures:

CRS and bounds catch mislabelled coordinates, which make every measurement wrong.
Geometry health catches rows that will disappear from joins without being counted.
Attribute nulls catch rows that groupby will drop from totals.
Duplicates catch double-counting, which is invisible in a map and obvious in a sum.
Extent outliers catch parse failures that put features at (0, 0).
Precomputed columns catch stale values that were correct when someone else computed them.

There is a second reason to do it systematically rather than ad hoc: the profile is evidence. Data arrives from suppliers, gets edited by colleagues, and changes between deliveries. A saved profile from the day the file arrived turns "this used to work" from an argument into a diff.

The ordering is not arbitrary either. CRS comes before anything measured, because area and length are meaningless until it is settled. Geometry health comes before joins, because empty geometries match nothing. Plotting comes last, because a map is the least reliable diagnostic on this page — it renders a broken layer and a clean one identically.

Wrong CRSs, first-N reads and truncated shapefile names

gdf.crs can be present and wrong. Always cross-check against total_bounds.
read_file(rows=N) reads the first N in file order, which is not random. For a genuine sample use gdf.sample(n) after a full read, or read with a bbox.
centroid on a geographic CRS warns and gives a slightly wrong point. Fine for outlier detection, not for analysis.
Shapefile column names are truncated to 10 characters — population_density arrives as populatio. See truncated column names.
Shapefile has no true date type, so date columns arrive as strings. Check dtypes before comparing.
A .prj file can be missing from a shapefile, giving crs=None even though the data is projected.
fiona.listlayers on a GDB needs the OpenFileGDB driver, which is present in modern GDAL builds.
Encoding matters for text attributes. Garbled names usually mean the file is cp1252 and was read as UTF-8 — see shapefile encoding errors.

Internal links

How to set up a Python GIS environment that actually works — before any of this can run
Null, empty, missing and invalid: four kinds of broken geometry — question three in depth
How to build an inventory of a GIS data folder in Python — profiling at folder scale
The Python GIS data cleaning checklist — what to do about what you find
How to validate a GeoDataFrame against a schema before analysis — turning findings into checks
Coordinate reference systems (CRS) explained for Python GIS — question two, properly
How to find and remove duplicate geometries in GeoPandas — question five
Spatial data quality: the six dimensions that matter — the framework behind the checklist

FAQ

What is the very first thing to check?

The CRS, cross-checked against total_bounds. Everything measured depends on it, and a mislabelled CRS never raises an error.

How do I know if a file has multiple layers?

fiona.listlayers(path). GeoPackages and File Geodatabases routinely hold several; read_file without a layer argument takes the first.

The file is 4 GB. How do I look at it without loading it?

fiona.open() gives you CRS, schema, bounds and feature count without reading geometry. Then read_file(..., rows=1000) or a bbox for a sample.

Should I trust a `SHAPE_Area` column?

Not without checking. Compare it against gdf.geometry.area: a stored value can predate an edit, or have been computed in a different CRS by the software that exported it. The ONS ward areas checked on this page agreed exactly — which is the answer you want, and not one any file promises.

What if there is no CRS?

Do not guess. Check total_bounds for a clue about the units, then ask the supplier. set_crs with the wrong value silently relocates the data.

How do I spot a mislabelled CRS?

Bounds between -180 and 180 mean degrees. If the CRS claims a projected system in metres and the bounds look like degrees, the label is wrong.

Is plotting a good first check?

It is a good last check. A map of a layer with empty geometries, duplicates or a wrong CRS often looks completely normal.