⚙️ How-To Getting started Data cleaning Updated

How to Explore a Spatial Dataset You Have Never Seen Before

Problem statement

A colleague sends boundaries_final_v3.zip. There is no documentation, no metadata, and no reply to your email.

gdf = gpd.read_file("boundaries_final_v3.shp")
gdf.head()

Five rows of a table with columns called OBJECTID, NAME_1, TYPE_CD and SHAPE_Area. You still do not know what the layer covers, what coordinate system it is in, whether SHAPE_Area means anything, whether the geometry is usable, or how many of the 12,400 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

gdf = gpd.read_file("boundaries_final_v3.shp")

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   (12400, 9)
2. CRS         EPSG:27700
3. bounds      [ 82672. 5342. 655604. 657534.]
4. geometry    {'Polygon': 12103, 'MultiPolygon': 297}
5. health      null=6 empty=214 invalid=31
6. columns     {'OBJECTID': 'int64', 'NAME_1': 'object', 'SHAPE_Area': 'float64', ...}
7. nulls       {'NAME_1': 842, 'TYPE_CD': 0, 'SHAPE_Area': 6}

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
mixed Polygon/MultiPolygon a geometry(Polygon) column will reject 297 rows
214 empty geometries dropna() will not catch them, and they match nothing
842 null names any groupby("NAME_1") silently loses them
a SHAPE_Area column precomputed, in unknown units, possibly stale

The seven questions

Checklist of the seven questions to ask a new spatial layer before analysis.
Twenty minutes here, or a wrong number later. The order matters: CRS before anything measured.

Step-by-step solution

Vertical steps from listing layers through CRS, geometry health, attributes and a first 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.gpkg"))
# ['parcels', 'wards', 'roads_2023', 'roads_2024']

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

with fiona.open("data.gpkg", layer="parcels") as src:
    print(src.crs)
    print(src.schema)
    print(len(src))
    print(src.bounds)
# {'init': 'epsg:27700'}
# {'geometry': 'Polygon', 'properties': OrderedDict([('OBJECTID', 'int:10'), ...])}
# 12400
# (82672.0, 5342.0, 655604.0, 657534.0)

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

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

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

print(gdf.crs)                      # EPSG:27700
print(gdf.crs.is_projected)         # True
print(gdf.crs.axis_info[0].unit_name)  # metre
print(gdf.total_bounds)             # [82672. 5342. 655604. 657534.]

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":    gdf.geometry.isna().sum(),
    "empty":   gdf.geometry.is_empty.sum(),
    "invalid": (~gdf.is_valid).sum(),
    "types":   gdf.geom_type.value_counts().to_dict(),
}

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())
# Self-intersection    28
# Too few points        3

Twenty-eight self-intersections in 12,400 rows is normal digitising noise. Twenty-eight hundred 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. {s.dropna().iloc[0] if s.notna().any() else '—'!r}")
OBJECTID       int64    nulls=    0 unique= 12400  e.g. 1
NAME_1         object   nulls=  842 unique=   312  e.g. 'Leith Walk'
TYPE_CD        object   nulls=    0 unique=     4  e.g. 'RES'
SHAPE_Area     float64  nulls=    6 unique= 12180  e.g. 4211.883

Three things to look for:

# a) a column that should be a key — is it?
gdf["OBJECTID"].is_unique                # True → usable as an id

# b) low-cardinality columns are categories; find out what the codes mean
gdf["TYPE_CD"].value_counts()
# RES    8102 / COM 2911 / IND 1102 / MIX 285

# c) precomputed geometry columns are usually stale
(gdf["SHAPE_Area"] - gdf.geometry.area).abs().describe()
# max  18442.6   ← they disagree; the stored value predates an edit

That last check is worth running every time. SHAPE_Area, Shape_Leng, area_ha and friends were computed at some point in the past, in some CRS, and are frequently wrong after editing. Recompute rather than trust — and if they disagree, that is information about how the file was produced.

5. Look for duplicates, in both senses 

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

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

Forty-four identical shapes usually means a merge happened twice. They double every area total and produce duplicate rows in every join. See how to find and remove duplicate geometries.

6. Check the spatial extent makes sense 

print(gdf.total_bounds)
# [ 82672. 5342. 655604. 657534.]

That bounding box is 573 km wide — the whole of Great Britain. For a layer named "boundaries" that may be right, or it may mean 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 far from the rest")
print(outliers.total_bounds)

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.plot(ax=ax2, column="TYPE_CD", legend=True, categorical=True)
ax2.set_title("by TYPE_CD")
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 214 empty geometries looks exactly like a plot of a layer with none.

Code examples

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, "boundaries_v3"), indent=2, default=str))

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.

Example 2: profiling every layer in a folder 

from pathlib import Path

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"]]
print(summary.to_string(index=False))

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="OBJECTID"):
    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"

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.

Explanation

Triage rows pairing an unchecked assumption with the wrong result it silently produces.
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.

Edge cases or notes

  • 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.

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?

No. Compare it against gdf.geometry.area — they disagree surprisingly often, because the stored value was computed before an edit or in a different CRS.

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.