How to Find and Remove Spatial Outliers in a Point Dataset

Problem statement

Sensor.community publishes 24-hour averages from its citizen air-quality sensors, each with the country its owner registered it in. Take the particulate sensors registered in the Netherlands and look at the extent:

print(len(gdf), "sensors registered in NL")
print(gdf.total_bounds.round(2))

1493 sensors registered in NL
[ 0.    0.    7.14 57.21]

A country 300 km across should not span 57 degrees of latitude. Some sensors sit at exactly (0, 0), some in the North Sea, a few just over the border or the coastline, and one reports a daily PM2.5 hundreds of times its neighbours'. 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:

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

gdf["outlier_reason"] = ""

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

# 2. outside the study area
area = netherlands.to_crs(gdf.crs).union_all()
outside = ~gdf.geometry.within(area) & gdf["outlier_reason"].eq("")
gdf.loc[outside, "outlier_reason"] = "outside study area"

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

outlier_reason
                       1478
outside study area       12
missing/null-island       3
Name: count, dtype: int64
1478 of 1493 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. Several of the sensors "outside" here are a few hundred metres over a generalised coastline, as the next sections show.

Five kinds of spatial outlier

From impossible coordinates to local anomalies

Start with the checks that need no statistics

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 | x.isna() | y.isna(),
        "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

flagged = basic_flags(gdf)
print(flagged["flag"].replace("", "(none)").value_counts())
print(flagged.loc[flagged["flag"].ne(""), ["sensor_id", "lat", "lon", "pm25"]])

flag
(none)         1490
null island       3
Name: count, dtype: int64
     sensor_id  lat  lon  pm25
141      35977  0.0  0.0  0.75
142      35981  0.0  0.0  2.99
920      77160  0.0  0.0  2.79

These rules are absolute — a latitude of 95 is not an unusual observation, it is an error — so they can be applied without judgement. Here they catch only the three sensors at Null Island: nothing in this feed has an impossible latitude, and every sensor has some coordinate.

Clip to a reference boundary

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

boundary = netherlands.to_crs(gdf.crs).union_all()
live = flagged[flagged["flag"].eq("")]

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

# how far outside? a tolerance band separates "just over the line" from "somewhere else"
metric_crs = live.estimate_utm_crs()
outside = live[~inside].to_crs(metric_crs)
edge = gpd.GeoSeries([boundary], crs=gdf.crs).to_crs(metric_crs).iloc[0]
outside["km_outside"] = (outside.distance(edge) / 1000).round(2)
print(outside.sort_values("km_outside")[["sensor_id", "lat", "lon", "km_outside"]].to_string(index=False))

12 points outside the study area
 sensor_id       lat      lon  km_outside
     52928 53.173277 5.420108        0.01
     96254 52.111478 4.279996        0.14
     98436 52.098000 4.262000        0.15
     76423 52.104000 4.268000        0.25
     67080 51.041880 5.857735        0.26
     11531 51.000000 5.874000        0.36
     82632 51.444000 3.582000        0.55
     52456 52.840000 7.082000        0.67
     65881 51.035718 5.866678        0.87
     49592 51.482000 4.264000        1.92
     87370 57.060000 6.188000      390.29
     99793 57.212000 3.374000      439.93

Distinguishing "just outside" from "in another country" matters. Ten sensors are within 2 km of the line — on The Hague's seafront, at Harlingen harbour, along the Maas where it is the Belgian border. The boundary is Natural Earth's 1:10 million outline, and at that scale a coastline or a river border is generalised by hundreds of metres. Two sensors are in the North Sea, 390 and 440 km off — a real error, not an imprecision.

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.

from sklearn.neighbors import NearestNeighbors

kept = live[inside].copy()
metric = kept.to_crs(metric_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

kept["knn_dist_m"] = mean_knn.round()
kept["distance_outlier"] = mean_knn > threshold
print(f"median {np.median(mean_knn):,.0f} m, threshold {threshold:,.0f} m — "
      f"{kept['distance_outlier'].sum()} points flagged")
print(kept[kept["distance_outlier"]].sort_values("knn_dist_m", ascending=False)
      [["sensor_id", "lat", "lon", "knn_dist_m"]].head(6).to_string(index=False))

median 2,053 m, threshold 14,103 m — 14 points flagged
 sensor_id       lat      lon  knn_dist_m
     55692 53.366000 5.212000     21300.0
     87849 52.929890 4.767083     20978.0
     54271 53.362000 5.216000     20888.0
     99057 51.974000 6.708000     20236.0
     30605 51.328000 3.860000     19649.0
     56700 53.447682 5.662346     19552.0

The most isolated sensors are on Terschelling and Ameland, near Den Helder, at Winterswijk on the German border and in Zeeuws-Vlaanderen — correctly located, just a long way from the next sensor. That is what this test finds: isolation, which is evidence for a bad coordinate only when isolation is implausible for the data. Always compute distances in a projected CRS; degrees are not distances.

Flag local attribute anomalies

A reading can be unremarkable nationally and impossible for its street. Here, the median sensor reported a daily PM2.5 of about 4 µg/m³, and a few reported hundreds.

RADIUS = 5000        # metres

points = metric.assign(sensor_id=kept["sensor_id"].values)[["sensor_id", "pm25", "geometry"]]
neighbourhoods = gpd.sjoin(
    points,
    points.assign(geometry=points.buffer(RADIUS)),
    how="inner", predicate="within",
)
neighbourhoods = neighbourhoods[neighbourhoods.index != neighbourhoods["index_right"]]   # not itself

stats = neighbourhoods.groupby(level=0)["pm25_right"].agg(["mean", "median", "std", "count"])
kept["local_median"] = stats["median"]
kept["local_n"] = stats["count"].reindex(kept.index).fillna(0).astype(int)
kept["local_z"] = (kept["pm25"] - stats["mean"]) / stats["std"].replace(0, np.nan)

kept["attribute_outlier"] = kept["local_z"].abs().gt(3) & kept["local_n"].ge(5)
print(f"{kept['attribute_outlier'].sum()} local attribute anomalies")
print(kept[kept["attribute_outlier"]].sort_values("pm25", ascending=False)
      [["sensor_id", "sensor_type", "pm25", "local_median", "local_n"]].head(6).to_string(index=False))

at_ceiling = kept["pm25"].eq(999.9)
print(f"at exactly 999.9: {at_ceiling.sum()} sensors, {(at_ceiling & kept['attribute_outlier']).sum()} of them flagged")
print("reporting exactly 0:", int(kept["pm25"].eq(0).sum()))

47 local attribute anomalies
 sensor_id sensor_type   pm25  local_median  local_n
     78460       SPS30 3459.5         6.295       54
     72827      SDS011  999.9         3.065       26
     65457      SDS011  999.9         5.340       20
     52888      SDS011  999.9         5.505        6
     53273      SDS011  999.9         3.310        8
     38397      SDS011  999.9         5.780       13
at exactly 999.9: 14 sensors, 6 of them flagged
reporting exactly 0: 20

The top of the list is not air quality. Fourteen SDS011 sensors averaged exactly 999.9 over 24 hours, with PM10 at exactly 1999.9: the ceiling the sensor reports, held all day, which is a stuck or saturated sensor. Twenty report exactly 0. And the z-score caught only six of the fourteen stuck sensors. Seven have fewer than five neighbours within 5 km, so the minimum count — there to stop sparse areas producing meaningless z-scores — skipped them; the eighth has a broken neighbour of its own, which inflates the local spread it is measured against. That is masking, and it is the argument for turning what you know about the instrument into an absolute rule in step one — pm25 == 999.9 needs no neighbourhood — and keeping the statistics for what you cannot write down.

Detect gridding and duplicate coordinates

Rounded or geocoded coordinates produce towers of identical points that look like a hotspot, and they limit what the data can support.

counts = kept.groupby(["lat", "lon"]).size()
print("most repeated coordinates:")
print(counts.sort_values(ascending=False).head(3))

# coordinate precision: how many decimals does each sensor actually carry?
decimals = kept["lat_text"].str.split(".").str[1].str.len()
print(pd.crosstab(decimals.clip(upper=6).rename("decimals"), kept["exact_location"].rename("exact_location")))

most repeated coordinates:
lat     lon  
51.660  5.824    4
51.774  4.612    3
51.198  6.010    3
dtype: int64
exact_location    0    1
decimals                
1.0              14    0
2.0             161    0
3.0             575    0
6.0               0  726

Half the sensors carry three decimal places or fewer, and they are exactly the ones whose owners did not tick "exact location": Sensor.community rounds those coordinates for privacy. Three decimals of latitude is about 110 m — so a 50 m buffer analysis, or a "nearest sensor to this street" lookup, is working at a precision half the data does not have. Parse precision from the text as delivered; a float such as 51.66 printed back from a GeoSeries cannot tell you whether the source said 51.66 or 51.660.

Keep the reasons, and review before deleting

flagged["reason"] = flagged["flag"]
flagged.loc[outside.index, "reason"] = np.where(outside["km_outside"] <= 2, "just outside boundary", "outside study area")
flagged.loc[kept.index[kept["distance_outlier"]], "reason"] = "isolated point"
flagged.loc[kept.index[kept["attribute_outlier"]], "reason"] = "local value anomaly"
flagged["is_outlier"] = flagged["reason"].ne("")

print({"total": len(flagged), "flagged": int(flagged["is_outlier"].sum()),
       "by reason": flagged.loc[flagged["is_outlier"], "reason"].value_counts().to_dict()})

Path("out").mkdir(exist_ok=True)
flagged[flagged["is_outlier"]].drop(columns="flag").to_file("out/outliers_for_review.gpkg")
flagged[~flagged["is_outlier"]].drop(columns="flag").to_file("out/sensors_clean.gpkg")

{'total': 1493, 'flagged': 76, 'by reason': {'local value anomaly': 47, 'isolated point': 14, 'just outside boundary': 10, 'null island': 3, 'outside study area': 2}}

Two outputs, both written: the clean layer for analysis and the flagged layer for the person who knows the data. Only the first two reasons are errors; the rest are candidates.

Screening the Sensor.community air-quality feed

The data used on this page

Sensor.community's 24-hour average feed (data.24h.json, open data under the Open Database Licence), saved on 30 September 2026 — the feed is live, so a download today gives different sensors and different counts — and Natural Earth's 1:10m country boundaries (public domain). Natural Earth's Netherlands includes the three Caribbean islands, so the setup keeps only the European part: a study area that spans the Atlantic would make every bounding-box check meaningless. The sensors are the ones with a PM2.5 reading (P2) registered with the country code NL.

import json
from pathlib import Path

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

Path("data").mkdir(exist_ok=True)
feed = Path("data/sensors_24h.json")
if not feed.exists():
    # live: https://data.sensor.community/static/v2/data.24h.json, the last 24 hours of every sensor;
    # read instead as saved on 30 September 2026, so the sensors stay the ones measured here
    feed.write_bytes(requests.get("https://data.spatialworkflow.io/sensor-community-24h-2026-09-30.json",
                                  timeout=120).content)
countries_zip = Path("data/ne_10m_admin_0_countries.zip")
if not countries_zip.exists():
    countries_zip.write_bytes(requests.get(
        "https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip", timeout=120).content)

rows = []
for record in json.loads(feed.read_text()):
    values = {v["value_type"]: v["value"] for v in record["sensordatavalues"]}
    where = record["location"]
    if "P2" in values and where["country"] == "NL":
        rows.append({"sensor_id": record["sensor"]["id"], "sensor_type": record["sensor"]["sensor_type"]["name"],
                     "lat_text": where["latitude"], "lon_text": where["longitude"],
                     "exact_location": where["exact_location"], "indoor": where["indoor"],
                     "pm25": values["P2"]})
sensors = pd.DataFrame(rows)
sensors["lat"] = pd.to_numeric(sensors["lat_text"], errors="coerce")
sensors["lon"] = pd.to_numeric(sensors["lon_text"], errors="coerce")
sensors["pm25"] = pd.to_numeric(sensors["pm25"], errors="coerce")
gdf = gpd.GeoDataFrame(sensors, geometry=gpd.points_from_xy(sensors["lon"], sensors["lat"]), crs=4326)

world = gpd.read_file(countries_zip)
kingdom = world[world["ISO_A2_EH"] == "NL"]
print(len(gdf), "sensors |", kingdom["ADMIN"].iloc[0], "extends to", kingdom.total_bounds.round(1))

# Natural Earth's Netherlands includes Bonaire, Sint Eustatius and Saba; the study area is the European part
parts = kingdom.explode(index_parts=False)
netherlands = parts[parts.geometry.representative_point().x > 0].dissolve()
print("European part:", netherlands.total_bounds.round(1), f"({len(parts)} parts, {len(parts) - len(netherlands.explode())} in the Caribbean)")

1493 sensors | Netherlands extends to [-68.4  12.    7.2  53.6]
European part: [ 3.3 50.7  7.2 53.6] (12 parts, 3 in the Caribbean)

Example 1: a complete, reusable outlier screen

from sklearn.neighbors import NearestNeighbors

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

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

    live = out[out["flag"].eq("")]
    metric_crs = live.estimate_utm_crs()
    metric = live.to_crs(metric_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.distance(b) > 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"].isin(["", "just outside boundary"])]
    if len(live) > k + 1:
        m = live.to_crs(metric_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])
        out.loc[live.index, "knn_dist_m"] = mean_knn.round()
        isolated = live.index[(mean_knn > q3 + fence * (q3 - q1)) & live["flag"].eq("").to_numpy()]
        out.loc[isolated, "flag"] = "isolated point"

    out["is_outlier"] = out["flag"].isin(["null geometry", "null island", "impossible latitude",
                                          "impossible longitude", "outside study area"])
    return out

screened = screen_outliers(gdf, boundary=netherlands.union_all())
print(screened["flag"].replace("", "(none)").value_counts())
print("removed as errors:", int(screened["is_outlier"].sum()))

flag
(none)                   1465
isolated point             13
just outside boundary      10
null island                 3
outside study area          2
Name: count, dtype: int64
removed as errors: 5

The screen only removes what is certainly wrong. "Just outside" and "isolated" stay in the data with their flag, for someone to look at.

Example 2: a quick visual QA sheet

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

fig, axes = plt.subplots(1, 2, figsize=(11, 6))
minx, miny, maxx, maxy = netherlands.total_bounds
for ax, keep in zip(axes, (True, False)):
    netherlands.boundary.plot(ax=ax, linewidth=0.6, color="#94a3b8")
    part = screened[screened["flag"].eq("") == keep]
    part.plot(ax=ax, markersize=4 if keep else 14, color="#0ea5e9" if keep else "#ef4444")
    ax.set_xlim(minx - 0.3, maxx + 0.3); ax.set_ylim(miny - 0.3, maxy + 0.3)   # the far-off ones are off-frame
    ax.set_title(f"{'no flag' if keep else 'flagged'}: {len(part)}")
fig.savefig("out/outlier_qa.png", dpi=130, bbox_inches="tight")
plt.close(fig)
print(Path("out/outlier_qa.png").stat().st_size > 0)

True

A before/after pair takes seconds to produce and answers "did I just delete a whole valid region?" better than any statistic — here, the flagged panel shows the isolated points sitting on the islands and at the edges of the country, where a reader who knows the Netherlands can dismiss them at a glance.

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.

from esda.moran import Moran_Local
from libpysal.weights import KNN

outdoor = gpd.read_file("out/sensors_clean.gpkg")
outdoor = outdoor[outdoor["indoor"].eq(0)].to_crs(outdoor.estimate_utm_crs()).reset_index(drop=True)
w = KNN.from_dataframe(outdoor, k=8)
w.transform = "r"

lisa = Moran_Local(np.log1p(outdoor["pm25"].to_numpy()), w, permutations=999, seed=42)
outdoor["lisa_q"] = lisa.q                 # 1 HH, 2 LH, 3 LL, 4 HL
outdoor["lisa_p"] = lisa.p_sim

spatial_outliers = outdoor[(outdoor["lisa_p"] < 0.05) & outdoor["lisa_q"].isin([2, 4])]
print(f"{len(outdoor)} outdoor sensors: {len(spatial_outliers)} significant high-low / low-high outliers")
print(spatial_outliers["lisa_q"].map({2: "low among high", 4: "high among low"}).value_counts().to_dict())

1379 outdoor sensors: 67 significant high-low / low-high outliers
{'low among high': 36, 'high among low': 31}

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. PM2.5 is heavily skewed, so the statistic runs on log1p of it; on the raw values one extreme sensor dominates the variance. At p < 0.05 with over a thousand tests, some of these are chance: treat the list as candidates, not findings.

Example 4: DBSCAN for cluster-and-noise separation

from sklearn.cluster import DBSCAN

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

for eps in (1000, 3000, 5000):
    labels = DBSCAN(eps=eps, min_samples=5).fit_predict(coords)     # eps in metres
    n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
    print(f"eps {eps:>5} m: {(labels == -1).sum():>4} points labelled noise, {n_clusters} clusters")

eps  1000 m:  987 points labelled noise, 36 clusters
eps  3000 m:  578 points labelled noise, 40 clusters
eps  5000 m:  261 points labelled noise, 46 clusters

DBSCAN treats "noise" as points in no dense region, which suits event or observation data where the real signal is clustered. The answer depends entirely on eps: at 1 km, 987 of the 1,379 outdoor sensors are "noise"; at 5 km, 261 are. That is the parameter to justify and record, not the output.

An outlier is relative to an expectation

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

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 country. 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 hundred metres outside, so a tolerance band separates "just over the line" from "in the North Sea".

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 sensor on Terschelling is legitimately isolated. A local z-score, or Local Moran's I, says a value is unlike its neighbours — which is sometimes a broken sensor 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: "1,493 sensors, 5 removed as errors, the rest flagged for review" is a sentence you can defend, and "1,488 sensors" on its own is not.

Projection, Null Island and real extremes

• 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, or rounded for privacy as half these sensors are, stack up legitimately. Detect them with a duplicate-coordinate count and a precision check, not a distance test.
• Boundary data has its own errors: A generalised coastline can put a genuine coastal point in the sea — ten sensors here fall within 2 km outside Natural Earth's 1:10m outline. Use the finest boundary you have for the clip, or a tolerance band.
• Thresholds need justification: Record the fence multiplier, radius, k and eps in the run summary, or the results are not reproducible — example 4 shows how far one parameter moves the answer.
• Removing outliers changes downstream statistics: Interpolation, density surfaces and kriging are all sensitive to what you dropped. State it in the metadata.
• Check the metadata columns too: Sensor.community marks indoor sensors with indoor = 1; an indoor reading is not an outdoor outlier, it is a different measurement, and example 3 drops them before testing.
• Some outliers are the point: In fraud, leak or fault detection, the anomalies are the signal. Screen them, do not discard them.

Internal links

• How to Handle Missing and Null Values in Spatial Datasets in Python
• How to Clean Messy CSV Coordinates into a Reliable GeoDataFrame
• My Points Plot in the Ocean: Fixing Swapped Latitude and Longitude
• How to Find the Nearest Point in GeoPandas
• How to Select Features by Location in GeoPandas
• The Python GIS Data Cleaning Checklist: From Raw Download to Analysis-Ready
• How to clean yield monitor data in Python — structural artefacts rather than outliers

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 ($Q_3 + 3 \times \mathrm{IQR}$) adapts to the density of your dataset. Record whatever you used. Expect it to flag places that are genuinely remote — on this page, sensors on the Wadden islands.

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.