The Python GIS Data Cleaning Checklist: From Raw Download to Analysis-Ready
On this page
You downloaded a dataset, opened it in Python, and ran your first spatial join — and got an empty result, a TopologyException, or a silent reprojection that placed your features in the ocean. Welcome to spatial data cleaning, the unglamorous work that decides whether the rest of your analysis is trustworthy.
This guide gives you a complete, ordered mental model for turning a raw download into an analysis-ready GeoDataFrame using current GeoPandas (1.x) and Shapely (2.x). Work through it top to bottom the first few times. Once the order is in your head, you will spot which step a given dataset actually needs.
Problem statement
Two real downloads for the Isle of Wight — the Food Standards Agency's register of 1,701 food businesses and the 86,080 polygons in OpenStreetMap — carry most of the defects spatial data arrives with, and none of them raises an error on its own. This page cleans both and counts what it finds at every step:
- Null geometries — 375 of the food businesses have no coordinates at all.
- Coordinates as text — the register gives latitude and longitude as strings.
- Invalid geometries — 19 OpenStreetMap polygons are self-intersecting, and repairing them the default way turns every one into a
GeometryCollection. - Duplicates that are not errors — 10 polygons share an outline with another, and in four of those five pairs the two are genuinely different features.
- Inconsistent attributes — empty strings meaning "missing", a rating column that mixes numbers with
ExemptandAwaitingInspection, and postcodes in more than one format.
Quick answer
Clean spatial data in a fixed order: load and inspect, build or fix the geometry and CRS, remove null and empty geometries, repair invalid ones, normalise attributes, deduplicate on a real key, validate, then export to GeoPackage — counting rows after every step, so that any drop is one you meant.
Here is the whole pipeline for the food register in one block:
import json
import geopandas as gpd
import pandas as pd
records = pd.DataFrame(json.load(open("data/fsa_isle_of_wight.json"))["establishments"])
lon = pd.to_numeric(records["geocode"].map(lambda g: (g or {}).get("longitude")), errors="coerce")
lat = pd.to_numeric(records["geocode"].map(lambda g: (g or {}).get("latitude")), errors="coerce")
points = gpd.points_from_xy(lon, lat)
points[(lon.isna() | lat.isna()).to_numpy()] = None # else POINT (NaN NaN), see below
shops = gpd.GeoDataFrame(records.drop(columns=["geocode", "scores", "links"], errors="ignore"),
geometry=points, crs="EPSG:4326")
counts = {"load": len(shops)}
shops = shops.to_crs(27700) # British National Grid, metres
shops = shops[shops.geometry.notna() & ~shops.geometry.is_empty].copy()
counts["drop null/empty"] = len(shops)
text = shops.select_dtypes(include=["object", "string"]).columns
shops[text] = shops[text].apply(lambda s: s.str.strip().replace({"": None}))
shops["rating"] = pd.to_numeric(shops["RatingValue"], errors="coerce").astype("Int64")
shops = shops.drop_duplicates(subset="FHRSID")
counts["dedupe on FHRSID"] = len(shops)
assert shops.crs.to_epsg() == 27700 and shops.geometry.is_valid.all() and shops["FHRSID"].is_unique
shops.to_file("output/isle_of_wight_clean.gpkg", layer="food_businesses", driver="GPKG")
print(counts)
{'load': 1701, 'drop null/empty': 1326, 'dedupe on FHRSID': 1326}
The data files are downloaded by the block in The data used on this page. The count dictionary is the most useful line in it: a cleaning run that does not report what it removed is indistinguishable from one that lost data.
The cleaning checklist
The checklist, from raw download to export
Load & inspect
Read the files and look before you touch anything:
raw = pd.DataFrame(json.load(open("data/fsa_isle_of_wight.json"))["establishments"])
print(len(raw), "food businesses;", raw["geocode"].map(lambda g: type((g or {}).get("latitude")).__name__).value_counts().to_dict())
print(raw["RatingValue"].value_counts().to_dict())
osm = gpd.read_file("data/isle_of_wight_polygons.gpkg")
print(len(osm), "OSM polygons;", osm.crs, osm.geom_type.value_counts().to_dict())
print("null:", osm.geometry.isna().sum(), "| empty:", osm.geometry.is_empty.sum(), "| invalid:", (~osm.is_valid).sum())
1701 food businesses; {'str': 1326, 'NoneType': 375}
{'5': 1035, '4': 276, '3': 137, 'AwaitingInspection': 135, 'Exempt': 77, '2': 32, '1': 7, '0': 1, 'AwaitingPublication': 1}
86080 OSM polygons; EPSG:4326 {'MultiPolygon': 86080}
null: 0 | empty: 0 | invalid: 19
Two findings already shape the rest: the register's coordinates are strings (and absent for some rows), and its rating column is text that mostly holds digits.
Build the geometry and confirm the CRS
A register like this arrives as a table, not a spatial file. Convert the coordinates to numbers and build points with longitude as x:
lon = pd.to_numeric(raw["geocode"].map(lambda g: (g or {}).get("longitude")), errors="coerce")
lat = pd.to_numeric(raw["geocode"].map(lambda g: (g or {}).get("latitude")), errors="coerce")
shops = gpd.GeoDataFrame(raw.drop(columns=["geocode", "scores"], errors="ignore"),
geometry=gpd.points_from_xy(lon, lat), crs="EPSG:4326")
print("missing coordinates:", lat.isna().sum())
print("geometry of a row without them:", shops.geometry[lat.isna()].iloc[0],
"| isna:", shops.geometry[lat.isna()].isna().all(), "| is_empty:", shops.geometry[lat.isna()].is_empty.all())
print("bounds after to_crs(27700):", shops.to_crs(27700).total_bounds.round(0))
missing coordinates: 375
geometry of a row without them: POINT (NaN NaN) | isna: False | is_empty: False
bounds after to_crs(27700): [211535. 75472. inf inf]
This is the trap in the whole checklist. points_from_xy() turns a missing coordinate pair into POINT (NaN NaN) — a geometry that is neither null nor empty, so both standard checks pass it. Reprojected, its coordinates become infinite, and the layer's bounds with them. Mark those rows as null yourself, then reproject and check the bounds again:
missing = (lon.isna() | lat.isna()).to_numpy()
shops.loc[missing, "geometry"] = None
shops, osm = shops.to_crs(27700), osm.to_crs(27700)
print("food businesses:", shops.total_bounds.round(0), "\nOSM polygons: ", osm.total_bounds.round(0))
food businesses: [211535. 75472. 465726. 205636.]
OSM polygons: [428948. 75228. 466312. 96666.]
The polygons fit the island, but the businesses reach 211 km west and 205 km north of the grid's origin — the island spans roughly 429–466 km east. Bounds are a coarse check; the precise one is distance from the island itself:
island = gpd.read_file("data/isle_of_wight_boundary.gpkg").union_all()
offshore = shops[shops.geometry.notna()].assign(km=lambda d: d.distance(island) / 1000)
offshore = offshore[offshore["km"] > 0].sort_values("km", ascending=False)
print(offshore[["BusinessName", "BusinessType", "PostCode", "km"]].round(1).to_string(index=False))
BusinessName BusinessType PostCode km
Burger King Restaurant/Cafe/Canteen 248.7
Wightlink - Pier Head Terminal Restaurant/Cafe/Canteen PO33 2HF 0.6
The Best Dressed Crab Takeaway/sandwich shop PO35 5NS 0.0
Three points fall outside the island's boundary. Two are right where they should be — a ferry terminal on its pier, and a business on the shoreline. The third, a Burger King with no postcode in the register, was geocoded 249 km away on the mainland: a wrong geocode, which no null, empty or validity check would ever find. Correct it from the postcode or set its geometry to null; here it is set to null, so the later steps drop it with the unlocated rows.
shops.loc[offshore.index[offshore["km"] > 10], "geometry"] = None
print("bounds now:", shops.total_bounds.round(0))
bounds now: [430036. 75472. 465726. 96591.]
Now the bounds put both layers on the Isle of Wight in British National Grid metres.
set_crs() only labels coordinates; use it when a file lost a CRS you know. to_crs() transforms them. For measuring on the Isle of Wight, British National Grid (EPSG:27700) is the right choice. Web Mercator (EPSG:3857), often used as a placeholder, would overstate areas there about 2.5 times — its area scale is 1/cos² of the latitude, and the island is at 50.7° N.
Drop null and empty geometries
Null and empty are distinct: notna() catches None, is_empty catches a geometry with no coordinates. Check both, and count:
import numpy as np
before = len(shops)
located = shops[shops.geometry.notna() & ~shops.geometry.is_empty].copy()
assert np.isfinite(located.geometry.x).all() and np.isfinite(located.geometry.y).all()
print(f"dropped {before - len(located)} of {before} rows without a usable location ({(before - len(located)) / before:.0%})")
print(shops.loc[shops.geometry.isna(), "BusinessType"].value_counts().head(3).to_dict())
dropped 376 of 1701 rows without a usable location (22%)
{'Mobile caterer': 124, 'Other catering premises': 115, 'Restaurant/Cafe/Canteen': 45}
The finite-coordinate assertion is what would have caught the POINT (NaN NaN) rows had they slipped through. Whether to drop the unlocated rows is a decision, not a formality: 22% of the register has no location, and those rows are not a random sample of it, as the business types show. A map or a spatial join cannot use them; a count of businesses by type should.
Check and repair validity
Find invalid geometries, see why they are invalid, then repair them:
from shapely import is_valid_reason
bad = osm[~osm.is_valid]
print(bad.geometry.map(is_valid_reason).str.split("[").str[0].value_counts().to_dict())
default = bad.geometry.make_valid()
typed = bad.geometry.make_valid(method="structure", keep_collapsed=False)
print("default method: ", default.geom_type.value_counts().to_dict())
print("method='structure':", typed.geom_type.value_counts().to_dict(),
"| largest area change:", f"{((typed.area - bad.area) / bad.area).abs().max():.1e}")
{'Self-intersection': 19}
default method: {'GeometryCollection': 19}
method='structure': {'Polygon': 16, 'MultiPolygon': 3} | largest area change: 1.4e-15
The default make_valid() (the "linework" method) returned a GeometryCollection for every one of the 19 — polygon pieces plus the stray lines where the rings crossed. A layer holding GeometryCollections cannot be written to a polygon layer, and many operations handle them poorly. method="structure" with keep_collapsed=False keeps only polygons, and here changed no area measurably. It needs GEOS 3.10 or later, which Shapely 2 wheels include.
osm.loc[bad.index, "geometry"] = typed
osm = osm[osm.geometry.notna() & ~osm.geometry.is_empty]
print("invalid after repair:", (~osm.is_valid).sum(), "| types:", osm.geom_type.value_counts().to_dict())
invalid after repair: 0 | types: {'MultiPolygon': 86064, 'Polygon': 16}
Normalise attribute columns
Standardise the text, turn the missing-value tokens into real missing values, and split mixed columns into typed ones. This comes before deduplication, because two rows that differ only by a trailing space or an empty string are duplicates only after it:
text = located.select_dtypes(include=["object", "string"]).columns
print("empty strings before:", int((located[text] == "").sum().sum()))
located[text] = located[text].apply(lambda s: s.str.strip().replace({"": None}))
print("empty strings after:", int((located[text] == "").sum().sum()))
located["rating"] = pd.to_numeric(located["RatingValue"], errors="coerce").astype("Int64")
located["rating_status"] = located["RatingValue"].where(located["rating"].isna(), "rated")
print(located["rating_status"].value_counts().to_dict())
pc = located["PostCode"].dropna()
print("postcodes without a space:", (~pc.str.contains(" ")).sum(), "e.g.", pc[~pc.str.contains(" ")].head(3).tolist())
empty strings before: 3305
empty strings after: 0
{'rated': 1195, 'Exempt': 74, 'AwaitingInspection': 55, 'AwaitingPublication': 1}
postcodes without a space: 0 e.g. []
select_dtypes(include="object") alone still finds text columns in pandas 3, but with a deprecation warning, because text now has its own string dtype. include=["object", "string"] selects them on pandas 2 and 3 alike; "str" works on pandas 3 but raises TypeError on pandas 2.2. The rating column becomes a number that can be averaged plus a status that says why some businesses have none.
Remove duplicates
Deduplicate on a key that identifies a real-world thing, and treat geometric duplicates as questions rather than deletions:
print("duplicate FHRSID:", located["FHRSID"].duplicated().sum())
shape = osm.geometry.normalize()
same_outline = osm[shape.duplicated(keep=False)].assign(outline=shape[shape.duplicated(keep=False)].to_wkb())
print(len(same_outline), "polygons share an outline with another")
pairs = same_outline.groupby("outline")[["name", "landuse", "natural", "building"]].agg(lambda s: "/".join(s.dropna().astype(str)))
print(pairs.reset_index(drop=True).head(5).to_string())
no_id = osm.drop(columns=["osm_id", "osm_way_id", "geometry"])
print("rows whose attributes repeat another's (ignoring IDs and geometry):", no_id.duplicated().sum())
duplicate FHRSID: 0
10 polygons share an outline with another
name landuse natural building
0 residential/residential
1 scrub/cliff
2 plant_nursery/grass
3 scrub/bare_rock
4 Carisbrooke Castle
rows whose attributes repeat another's (ignoring IDs and geometry): 83299
The register has a real key and no duplicates. The polygons are the instructive case: of the five pairs that share an outline, four carry different tags — scrub and a cliff, a plant nursery and grass — so they are two things mapped on one boundary, and deleting either loses a real feature. Only the residential pair looks like a genuine duplicate, and the place to fix that is OpenStreetMap. An attribute-only drop_duplicates() without the ID columns would have been far worse: 84,183 rows — mostly buildings tagged only building=yes — repeat another row's attributes. Deduplicate on an ID, or on attributes and geometry together, never on attributes alone.
Validate against expectations
Before exporting, assert what must be true:
checks = {
"CRS is British National Grid": located.crs.to_epsg() == 27700 and osm.crs.to_epsg() == 27700,
"no null or empty geometry": bool(located.geometry.notna().all() and (~located.geometry.is_empty).all()),
"all coordinates finite": bool(np.isfinite(located.geometry.x).all() and np.isfinite(located.geometry.y).all()),
"all polygons valid": bool(osm.is_valid.all()),
"only polygon types": set(osm.geom_type) <= {"Polygon", "MultiPolygon"},
"FHRSID unique": located["FHRSID"].is_unique,
"ratings within 0-5": bool(located["rating"].dropna().between(0, 5).all()),
}
for name, ok in checks.items():
print(("PASS " if ok else "FAIL ") + name)
assert all(checks.values())
PASS CRS is British National Grid
PASS no null or empty geometry
PASS all coordinates finite
PASS all polygons valid
PASS only polygon types
PASS FHRSID unique
PASS ratings within 0-5
Export cleanly
Write to GeoPackage, and read it back:
located.to_file("output/isle_of_wight_clean.gpkg", layer="food_businesses", driver="GPKG")
osm.to_file("output/isle_of_wight_clean.gpkg", layer="osm_polygons", driver="GPKG")
for layer in ["food_businesses", "osm_polygons"]:
back = gpd.read_file("output/isle_of_wight_clean.gpkg", layer=layer)
print(layer, len(back), back.crs.to_epsg(), back.geom_type.value_counts().to_dict())
food_businesses 1325 27700 {'Point': 1325}
osm_polygons 86080 27700 {'MultiPolygon': 86080}
GeoPackage keeps both layers in one file, with their CRS, full field names and Unicode. A shapefile would have truncated LocalAuthorityBusinessID to ten characters and could not hold two layers.
Running it on the Isle of Wight food register
The data used on this page
- Food businesses — the Food Standards Agency's food hygiene rating register for the Isle of Wight (local authority 284), from its public API, saved as the JSON the API returns: 1,701 establishments in the copy saved on 10 October 2026. Businesses open and close, so the count changes. The food register is a copy of the API's answer saved on 10 October 2026.
- Island boundary — the Isle of Wight's local authority boundary (December 2024, full resolution) from the ONS, used to check where the points really are.
- OpenStreetMap polygons — the
multipolygonslayer of Geofabrik's Isle of Wight extract, read with GDAL's OSM driver and saved as a GeoPackage with the ID, name and three tag columns. OpenStreetMap is edited daily.
import json
import os
import geopandas as gpd
import pandas as pd
import requests
os.makedirs("data", exist_ok=True)
os.makedirs("output", exist_ok=True)
if not os.path.exists("data/fsa_isle_of_wight.json"):
# live: requests.get("https://api.ratings.food.gov.uk/Establishments", headers={"x-api-version": "2"},
# params={"localAuthorityId": 284, "pageSize": 5000}, timeout=300) # Isle of Wight
# read instead as the API answered on 10 October 2026, so the register stays put
reply = requests.get("https://data.spatialworkflow.io/fhrs-api-284-2026-10-10.json", timeout=300)
reply.raise_for_status()
open("data/fsa_isle_of_wight.json", "w").write(reply.text)
if not os.path.exists("data/isle_of_wight_polygons.gpkg"):
pbf = "isle-of-wight-260101.osm.pbf"
if not os.path.exists(pbf):
with requests.get("https://download.geofabrik.de/europe/united-kingdom/england/isle-of-wight-260101.osm.pbf",
stream=True, timeout=300) as reply:
reply.raise_for_status()
with open(pbf, "wb") as out:
for chunk in reply.iter_content(1 << 20):
out.write(chunk)
polygons = gpd.read_file(pbf, layer="multipolygons",
columns=["osm_id", "osm_way_id", "name", "building", "landuse", "natural"])
polygons.to_file("data/isle_of_wight_polygons.gpkg")
if not os.path.exists("data/isle_of_wight_boundary.gpkg"):
reply = requests.get("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
"Local_Authority_Districts_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
params={"where": "LAD24CD='E06000046'", "outFields": "LAD24NM", "outSR": 27700, "f": "geojson"},
timeout=120)
reply.raise_for_status()
gpd.read_file(reply.text).to_file("data/isle_of_wight_boundary.gpkg")
print(len(json.load(open("data/fsa_isle_of_wight.json"))["establishments"]), "establishments,",
len(gpd.read_file("data/isle_of_wight_polygons.gpkg")), "polygons")
1701 establishments, 86080 polygons
The steps above are the code examples: run in order after this block, they are the complete cleaning run for both layers.
Why the steps come in this order
The order is not arbitrary — each step depends on the one before it being done. Inspecting first tells you which steps matter; there is no point repairing validity on a dataset that has none invalid. CRS comes early because reprojection touches every coordinate, and you want it settled before you compute anything geometric.
Null and empty removal precedes validity repair because make_valid() and is_valid expect real geometries to chew on; feeding them None or empty shapes wastes effort and clutters your masks. And because make_valid() can itself produce empties (when it collapses a degenerate sliver) or GeometryCollections, you re-run the null/empty filter after repair. That second pass is the step beginners forget.
Normalisation comes before deduplication because two rows that differ only by a trailing space, a case difference or an empty string meaning "missing" are duplicates only once those differences are gone; deduplicate first and they survive. Deduplication comes after repair for the same reason on the geometry side: two copies of one broken polygon become comparable only once both are valid. Validation and export close the loop: you assert your invariants — including the ones the standard checks miss, such as finite coordinates and a location inside the study area — then write to a format that preserves them.
Choices the checklist leaves to you
Valid is not the same as correct
make_valid() guarantees a geometry obeys the OGC rules — not that it represents reality. A repaired bowtie polygon is now topologically valid but may have an area you never intended. Validity is a floor, not a guarantee; spot-check repaired features, especially if only a few rows were affected but they carry a lot of weight.
Large files and memory
read_file() loads everything into memory. For large datasets, use the bbox or mask argument to read only your area of interest, the rows argument to sample, or read in chunks via pyogrio. Reproject and clean the subset, not the whole continent.
Re-checking after every transform
Repairs can introduce new defects: make_valid() with its default method returned GeometryCollections for all 19 invalid polygons here, reprojection turned POINT (NaN NaN) into infinite coordinates, and a join inflates duplicates. Treat the inspection block from step one as something you can re-run at any point. Cheap to run, expensive to skip.
CRS for area and distance
EPSG:3857 (Web Mercator) is often used as a placeholder, but its distances and areas are distorted away from the equator — about 2.5 times in area at the Isle of Wight's latitude. For real measurement use an equal-area or local projected CRS for your region (a UTM zone, a national grid, or an Albers equal-area). Pick the CRS for the question you are answering.
Internal links
- How to Fix Invalid Geometries in Python (GeoPandas)
- How to Remove Null and Empty Geometries in GeoPandas
- How to Fix CRS Mismatch in GeoPandas
- How to Reproject Spatial Data in Python (GeoPandas)
- How to Read a Shapefile in Python with GeoPandas
- How to Read and Write GeoPackage Files in Python
- How to run a privacy check before publishing a spatial dataset — the gate before release
FAQ
What is the difference between a null and an empty geometry?
A null geometry is None in the geometry cell — there is no geometry object at all, and notna() catches it. An empty geometry is a real object that contains no coordinates, such as POLYGON EMPTY; it passes a null check but fails area, length, and distance calculations. You must test for both with gdf.geometry.notna() & ~gdf.geometry.is_empty.
Should I use make_valid() or buffer(0) to fix geometries?
Use make_valid() from Shapely 2.x as your default — it is purpose-built and preserves your geometry's intent without distorting it. buffer(0) is an old trick that sometimes repairs polygons but can silently delete parts of a shape or change its area, and it does not help with lines or points. Reach for buffer(0) only as a fallback when make_valid() produces something you cannot use.
Why GeoPackage instead of Shapefile?
GeoPackage is a single SQLite-based file with no 2 GB size limit, no 10-character field-name truncation, reliable CRS and Unicode handling, and support for multiple layers. Shapefile is an older multi-file format with all of those limitations. Read shapefiles when you are handed them, but write your cleaned output to .gpkg with driver="GPKG".
In what order should I clean spatial data?
Inspect, fix and set the CRS, reproject, drop null and empty geometries, repair invalid geometries, re-drop any new empties, remove duplicates, normalise attributes, validate, then export. The order matters because each step assumes the previous one already ran — for example, validity repair expects non-null geometries, and you re-check for empties because make_valid() can create them.
How do I know if my CRS is actually correct?
Check gdf.crs to see the declared CRS, then check gdf.total_bounds and confirm the coordinate ranges match where the data should sit on Earth. Latitude/longitude data should fall within roughly -180 to 180 and -90 to 90; projected data will have large metre values. If the bounds put your features in the wrong place, the declared CRS is wrong — do not "fix" it with set_crs() unless you know the true source projection.
Why do I still have duplicate geometries after drop_duplicates()?
drop_duplicates() on attribute columns only catches rows whose attributes match — two records of the same shape with different IDs survive. Geometrically identical shapes can also differ in vertex order or starting point, so a raw comparison misses them. Use gdf.geometry.normalize().duplicated() to canonicalise geometries before comparing, which collapses logically identical shapes.
Do I need to reproject before calculating area or distance?
Yes, if your data is in a geographic CRS like EPSG:4326. Area and length computed on degrees are meaningless. Reproject to a metric CRS appropriate for your region — a UTM zone or an equal-area projection — before any measurement, and do it early so every downstream calculation uses correct units.