How to Find and Remove Spatial Outliers in a Point Dataset

Problem statement

A dataset of 12,000 sensor locations should sit inside one county. Plotted, it covers half a continent:

>>> gdf.total_bounds
array([-118.24,  33.79,   2.12,  55.95])

Three points are in California, one is at (0, 0), a dozen are stacked on the same rounded coordinate, and one sits 40 km offshore. None of them raise an error. They quietly stretch the map extent, distort every summary statistic, and inflate the search radius of every nearest-neighbour query you run afterwards.

Spatial outliers come in several flavours, and they need different treatment:

Out of area — coordinates far outside the expected region (wrong CRS, swapped lat/lon, a typo)
Null Island — exactly (0, 0), from missing values filled with zero
Geometric outliers — points far from every other point, but plausibly located
Attribute outliers in space — a value that is normal globally but anomalous for its neighbourhood
Duplicate or gridded coordinates — rounded to a coarse grid, or geocoded to a postcode centroid
In-the-wrong-place — a terrestrial observation that falls in water

Quick answer

Work from cheap, certain checks to expensive, statistical ones:

validate ranges and drop (0, 0) — no statistics needed
clip to a reference boundary and flag what falls outside
flag distance outliers with a k-nearest-neighbour distance threshold
flag local attribute anomalies with a neighbourhood z-score
flag, do not delete — keep a reason column and let the analyst decide

import geopandas as gpd
import numpy as np

gdf = gpd.read_file("data/raw/sensors.gpkg")
gdf["outlier_reason"] = ""

# 1. impossible or missing coordinates
bad = gdf.geometry.is_empty | gdf.geometry.isna()
null_island = gdf.geometry.x.eq(0) & gdf.geometry.y.eq(0)
gdf.loc[bad | null_island, "outlier_reason"] = "missing/null-island"

# 2. outside the study area
area = gpd.read_file("data/ref/county.gpkg").to_crs(gdf.crs).union_all()
outside = ~gdf.geometry.within(area) & (gdf["outlier_reason"] == "")
gdf.loc[outside, "outlier_reason"] = "outside study area"

print(gdf["outlier_reason"].value_counts())
clean = gdf[gdf["outlier_reason"] == ""].copy()
print(f"{len(clean)} of {len(gdf)} points kept")

Flagging rather than deleting is the important habit: a removed row cannot be reviewed, and "outlier" is often a judgement about the data-collection process rather than a fact about the point.

Five kinds of spatial outlier

Three map panels showing an out-of-area point, a distance outlier, and a local attribute anomaly. — Far away, far from its neighbours, or unlike its neighbours — three different problems.

Step-by-step solution

Vertical steps from range checks through boundary clip, kNN distance, local z-score to review. — Cheap and certain first; statistical last. Most datasets never need step four.

Start with the checks that need no statistics

import geopandas as gpd

def basic_flags(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    out = gdf.copy()
    out["flag"] = ""

    geographic = out.crs and out.crs.is_geographic
    x, y = out.geometry.x, out.geometry.y

    rules = {
        "null geometry": out.geometry.isna() | out.geometry.is_empty,
        "null island": x.eq(0) & y.eq(0),
    }
    if geographic:
        rules["impossible latitude"] = y.abs() > 90
        rules["impossible longitude"] = x.abs() > 180

    for reason, mask in rules.items():
        out.loc[mask & out["flag"].eq(""), "flag"] = reason
    return out

gdf = basic_flags(gpd.read_file("data/raw/sensors.gpkg"))
print(gdf["flag"].value_counts())

These rules are absolute — a latitude of 95 is not an unusual observation, it is an error — so they can be applied without judgement.

Clip to a reference boundary

The single most effective check is "is it where it should be?".

import geopandas as gpd

area = gpd.read_file("data/ref/county.gpkg").to_crs(gdf.crs)
boundary = area.union_all()

inside = gdf.geometry.within(boundary)
print(f"{(~inside).sum()} points outside the study area")

# a tolerance band catches points just over the line — often a digitising artefact
metric = gdf.to_crs(gdf.estimate_utm_crs())
buffered = gpd.GeoSeries([boundary], crs=gdf.crs).to_crs(metric.crs).buffer(250).iloc[0]
near_edge = ~inside & metric.geometry.within(buffered)
print(f"{near_edge.sum()} of those are within 250 m of the boundary")

Distinguishing "just outside" from "in another country" matters: the first is usually a coordinate-precision issue worth keeping, the second is a genuine error.

Flag distance outliers with k-nearest-neighbour distances

A point that is far from every other point is suspicious even when it is inside the study area.

import numpy as np
import geopandas as gpd
from sklearn.neighbors import NearestNeighbors

metric = gdf.to_crs(gdf.estimate_utm_crs())
coords = np.column_stack([metric.geometry.x, metric.geometry.y])

k = 5
nn = NearestNeighbors(n_neighbors=k + 1).fit(coords)      # +1: the point itself
dist, _ = nn.kneighbors(coords)
mean_knn = dist[:, 1:].mean(axis=1)                        # drop the self-distance

q1, q3 = np.percentile(mean_knn, [25, 75])
threshold = q3 + 3 * (q3 - q1)                             # Tukey fence, generous multiplier

gdf["knn_dist_m"] = mean_knn
gdf["distance_outlier"] = mean_knn > threshold
print(f"threshold {threshold:,.0f} m — {gdf['distance_outlier'].sum()} points flagged")

Without scikit-learn, shapely.STRtree or geopandas.sjoin_nearest give the same answer with a little more code. Always compute distances in a projected CRS; degrees are not distances.

Flag local attribute anomalies

A reading of 240 µg/m³ may be normal citywide and impossible for its street.

import numpy as np
import geopandas as gpd

metric = gdf.to_crs(gdf.estimate_utm_crs())
RADIUS = 2000        # metres

neighbourhoods = gpd.sjoin(
    metric[["value", "geometry"]],
    metric[["value", "geometry"]].assign(geometry=lambda d: d.buffer(RADIUS)),
    how="inner", predicate="within",
)

stats = neighbourhoods.groupby(level=0)["value_right"].agg(["mean", "std", "count"])
gdf["local_mean"] = stats["mean"]
gdf["local_z"] = (gdf["value"] - stats["mean"]) / stats["std"].replace(0, np.nan)
gdf["local_n"] = stats["count"]

gdf["attribute_outlier"] = gdf["local_z"].abs().gt(3) & gdf["local_n"].ge(5)
print(f"{gdf['attribute_outlier'].sum()} local attribute anomalies")

Requiring a minimum neighbour count stops sparse edges of the study area from producing meaningless z-scores.

Detect gridding and duplicate coordinates

Geocoding to postcode centroids produces towers of identical points that look like a hotspot.

coords = gdf.geometry.apply(lambda p: (round(p.x, 6), round(p.y, 6)))
counts = coords.value_counts()

print("most repeated coordinates:")
print(counts.head(5))
gdf["stacked_here"] = coords.map(counts)
gdf["stacked"] = gdf["stacked_here"] > 10

# coordinate precision: how many decimals are actually used?
decimals = gdf.geometry.x.astype(str).str.split(".").str[1].str.len()
print("x decimal places:", decimals.value_counts().head())

Two decimal places in a geographic CRS is roughly a kilometre of precision — worth knowing before you draw conclusions from a 200 m analysis.

Keep the reasons, and review before deleting

gdf["is_outlier"] = (
    gdf["flag"].ne("") | gdf["distance_outlier"] | gdf["attribute_outlier"]
)

summary = {
    "total": len(gdf),
    "flagged": int(gdf["is_outlier"].sum()),
    "by reason": (gdf.loc[gdf["is_outlier"], "flag"]
                  .replace("", "statistical").value_counts().to_dict()),
}
print(summary)

gdf[gdf["is_outlier"]].to_file("data/out/outliers_for_review.gpkg", driver="GPKG")
gdf[~gdf["is_outlier"]].to_file("data/out/sensors_clean.gpkg", driver="GPKG")

Two outputs, both written: the clean layer for analysis and the flagged layer for the person who knows the data.

Code examples

Example 1: a complete, reusable outlier screen

import numpy as np
import geopandas as gpd
from sklearn.neighbors import NearestNeighbors

def screen_outliers(gdf, boundary=None, k=5, fence=3.0, edge_tolerance_m=250):
    out = gdf.copy().reset_index(drop=True)
    out["flag"] = ""

    empty = out.geometry.isna() | out.geometry.is_empty
    out.loc[empty, "flag"] = "null geometry"

    valid = out.loc[~empty].copy()
    x, y = valid.geometry.x, valid.geometry.y
    out.loc[valid.index[x.eq(0) & y.eq(0)], "flag"] = "null island"
    if out.crs and out.crs.is_geographic:
        out.loc[valid.index[y.abs() > 90], "flag"] = "impossible latitude"
        out.loc[valid.index[x.abs() > 180], "flag"] = "impossible longitude"

    live = out[out["flag"].eq("") & ~empty]
    metric = live.to_crs(live.estimate_utm_crs())

    if boundary is not None:
        b = gpd.GeoSeries([boundary], crs=gdf.crs).to_crs(metric.crs).iloc[0]
        outside = ~metric.geometry.within(b)
        far = outside & ~metric.geometry.within(b.buffer(edge_tolerance_m))
        out.loc[live.index[far], "flag"] = "outside study area"
        out.loc[live.index[outside & ~far], "flag"] = "just outside boundary"

    live = out[out["flag"].eq("") & ~empty]
    if len(live) > k + 1:
        m = live.to_crs(live.estimate_utm_crs())
        coords = np.column_stack([m.geometry.x, m.geometry.y])
        dist, _ = NearestNeighbors(n_neighbors=k + 1).fit(coords).kneighbors(coords)
        mean_knn = dist[:, 1:].mean(axis=1)
        q1, q3 = np.percentile(mean_knn, [25, 75])
        threshold = q3 + fence * (q3 - q1)
        out.loc[live.index, "knn_dist_m"] = mean_knn
        out.loc[live.index[mean_knn > threshold], "flag"] = "isolated point"

    out["is_outlier"] = out["flag"].ne("")
    return out

screened = screen_outliers(
    gpd.read_file("data/raw/sensors.gpkg"),
    boundary=gpd.read_file("data/ref/county.gpkg").union_all(),
)
print(screened["flag"].value_counts())

Example 2: a quick visual QA sheet

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(13, 6))

area.boundary.plot(ax=axes[0], linewidth=0.6, color="#94a3b8")
screened[~screened["is_outlier"]].plot(ax=axes[0], markersize=4, color="#0ea5e9")
axes[0].set_title(f"kept: {(~screened['is_outlier']).sum()}")

area.boundary.plot(ax=axes[1], linewidth=0.6, color="#94a3b8")
screened[screened["is_outlier"]].plot(ax=axes[1], markersize=14, color="#ef4444")
axes[1].set_title(f"flagged: {screened['is_outlier'].sum()}")

fig.savefig("data/out/outlier_qa.png", dpi=130, bbox_inches="tight")
plt.close(fig)

A before/after pair takes seconds to produce and answers "did I just delete a whole valid region?" better than any statistic.

Example 3: local spatial autocorrelation with PySAL

For attribute outliers, a Local Moran's I identifies statistically significant spatial outliers — a high value surrounded by low ones and vice versa.

import geopandas as gpd
from libpysal.weights import KNN
from esda.moran import Moran_Local

gdf = gpd.read_file("data/out/sensors_clean.gpkg").to_crs(3857)
w = KNN.from_dataframe(gdf, k=8)
w.transform = "r"

lisa = Moran_Local(gdf["value"].values, w)
gdf["lisa_q"] = lisa.q                 # 1 HH, 2 LH, 3 LL, 4 HL
gdf["lisa_p"] = lisa.p_sim

spatial_outliers = gdf[(gdf["lisa_p"] < 0.05) & gdf["lisa_q"].isin([2, 4])]
print(f"{len(spatial_outliers)} significant high-low / low-high outliers")
spatial_outliers.to_file("data/out/lisa_outliers.gpkg", driver="GPKG")

Quadrants 2 and 4 are the spatial outliers proper — a value unlike its neighbourhood, which is precisely the anomaly a global z-score cannot see.

Example 4: DBSCAN for cluster-and-noise separation

import numpy as np
import geopandas as gpd
from sklearn.cluster import DBSCAN

metric = gdf.to_crs(gdf.estimate_utm_crs())
coords = np.column_stack([metric.geometry.x, metric.geometry.y])

labels = DBSCAN(eps=750, min_samples=5).fit_predict(coords)     # eps in metres
gdf["cluster"] = labels
noise = labels == -1
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
print(f"{noise.sum()} points labelled noise, {n_clusters} clusters")

DBSCAN treats "noise" as points in no dense region, which suits event or observation data where the real signal is clustered.

Explanation

"Outlier" is not a property of a point; it is a statement about a point relative to an expectation. That is why the useful checks fall into a hierarchy rather than a single test.

Grid comparing range checks, boundary clip, kNN distance, local z-score and DBSCAN on cost, certainty and what they catch. — Five methods, ordered by how much judgement they require of you.

At the top are checks that need no statistics at all. A latitude of 95° is impossible; a pair of exact zeros is a missing value that was filled with zero; an empty geometry is not a location. These are facts about the coordinate system and the data model, so they can be applied automatically and their results trusted.

Next comes the reference boundary, which encodes an expectation you already hold: these sensors are in this county. This single check catches wrong-CRS data, swapped axes, and typos in one pass, and it is cheap. The one subtlety is the boundary itself — coordinate precision and boundary generalisation mean genuine points can sit a few metres outside, so a small tolerance band separates "just over the line" from "in another country".

Below that, the tests become statistical, and their output is evidence rather than proof. A k-nearest-neighbour distance says a point is unusually isolated for this dataset; whether that makes it wrong depends on whether isolation is meaningful in your domain. A remote weather station is legitimately isolated. A local z-score, or Local Moran's I, says a value is unlike its neighbours — which is sometimes an error and sometimes the most interesting record in the file.

This is why flagging beats deleting. A flag column with a reason keeps the decision reversible, makes the cleaning auditable, and lets a domain expert overrule a threshold you picked. It also keeps the counts honest: "12,000 points, 143 flagged, 9 removed after review" is a sentence you can defend, and "11,857 points" on its own is not.

Edge cases or notes

Always project before measuring: Distances and buffers in EPSG:4326 are in degrees. estimate_utm_crs() gives a sensible metric CRS per dataset.
(0, 0) is a real place: It is in the Gulf of Guinea. Only the exact pair is suspicious; longitude 0 alone is Greenwich and perfectly normal.
A tight cluster is not an outlier: Points geocoded to a postcode centroid stack up legitimately. Detect them with a duplicate-coordinate count, not a distance test.
Boundary data has its own errors: A generalised coastline can put a genuine coastal point in the sea. Use the finest boundary you have for the clip.
Thresholds need justification: Record the fence multiplier, radius and k in the run summary, or the results are not reproducible.
Removing outliers changes downstream statistics: Interpolation, density surfaces and kriging are all sensitive to what you dropped. State it in the metadata.
Some outliers are the point: In fraud, leak or fault detection, the anomalies are the signal. Screen them, do not discard them.

Internal links

FAQ

Should I delete spatial outliers or flag them?

Flag them, with a reason, and write them to their own file. Deletion is irreversible and hides a decision that often needs domain knowledge to make correctly.

What is a good distance threshold for isolated points?

Derive it from the data rather than picking a number: a Tukey fence on the mean k-nearest-neighbour distance (Q3 + 3 × IQR) adapts to the density of your dataset. Record whatever you used.

How do I tell a data error from a genuinely remote observation?

Cross-check with a second signal — a study-area boundary, a land mask, or an attribute that should co-vary. A remote sensor with a plausible reading and a valid id is probably real; a remote point at exactly (0, 0) is not.

Why do my distances look wrong?

You are almost certainly measuring in degrees. Reproject to a metric CRS with gdf.to_crs(gdf.estimate_utm_crs()) before computing distances, buffers or areas.

What is the difference between a global and a local outlier?

A global outlier has an extreme value compared with the whole dataset. A local (spatial) outlier has a value unlike its neighbours, even if it is unremarkable globally — that is what Local Moran's I detects.

Do I need scikit-learn or PySAL for this?

No. Range checks and a boundary clip need only GeoPandas and catch most problems. Reach for sklearn for k-nearest-neighbour distances or DBSCAN, and esda for LISA statistics, when the simple checks are not enough.

How do I stop outliers reappearing next month?

Put the screen in the pipeline, write the flagged rows to a review file each run, and assert on the flagged count so a sudden jump fails the job instead of passing quietly.

How to Find and Remove Spatial Outliers in a Point Dataset #

Problem statement #

Quick answer #

Five kinds of spatial outlier #

Step-by-step solution #

Start with the checks that need no statistics #

Clip to a reference boundary #

Flag distance outliers with k-nearest-neighbour distances #

Flag local attribute anomalies #

Detect gridding and duplicate coordinates #

Keep the reasons, and review before deleting #

Code examples #

Example 1: a complete, reusable outlier screen #

Example 2: a quick visual QA sheet #

Example 3: local spatial autocorrelation with PySAL #

Example 4: DBSCAN for cluster-and-noise separation #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

Should I delete spatial outliers or flag them? #

What is a good distance threshold for isolated points? #

How do I tell a data error from a genuinely remote observation? #

Why do my distances look wrong? #

What is the difference between a global and a local outlier? #

Do I need scikit-learn or PySAL for this? #

How do I stop outliers reappearing next month? #