🔧 Fix Spatial statistics Joins & relationships Spatial Statistics · step 12 Notebook By · Updated

Spatial Weights Warn About Islands or Moran's I Returns nan

On this page

Problem statement

Building spatial weights for Cornwall's 337 neighbourhoods (see The data used on this page) produces a warning you were not expecting:

import warnings

from libpysal.weights import Queen

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    w = Queen.from_dataframe(zones, use_index=True)
print("\n".join(str(c.message) for c in caught if "connected" in str(c.message)))
The weights matrix is not fully connected: 
 There are 2 disconnected components.
 There is 1 island with id: E01019077.

Or the statistic itself comes back as nothing at all — here after joining Census 2011 population onto the 2021 neighbourhoods, a join that looks harmless:

from esda.moran import Moran

pop_2011 = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_144_1.data.csv?date=latest"
                       f"&geography={','.join(zones.index)}&rural_urban=0&cell=0&measures=20100"
                       "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["rate_2011"] = 1000 * zones["count"] / zones.index.map(pop_2011)

moran = Moran(zones["rate_2011"].to_numpy(), w, permutations=999)
print(f"I = {moran.I}  p = {moran.p_sim}")
('WARNING: ', 'E01019077', ' is an island (no neighbors)')
I = nan  p = 0.001

Look carefully at that second line. I is nan and the p-value is 0.001 — a normal, significant-looking number. Any pipeline that checks the p-value and reports "significant clustering" will sail straight past this.

These two problems look related and are not. Islands produce a warning and a quietly reduced analysis. nan comes from somewhere else entirely.

Quick answer

import numpy as np

values = zones["rate_2011"].to_numpy()
print(f"islands       {len(w.islands)}  {w.islands[:5]}")
print(f"components    {w.n_components}")
print(f"NaN values    {np.isnan(values).sum()}")
print(f"inf values    {np.isinf(values).sum()}")
print(f"variance      {np.nanvar(values):.6f}")
islands       1  ['E01019077']
components    2
NaN values    22
inf values    0
variance      33.580748
Symptom Cause Fix
"not fully connected" warning disconnected units drop them, or use KNN weights
islands with ids: … units with no neighbours same
I = nan a single NaN or inf in the values drop or impute before computing
I = nan zero variance (all values identical) there is nothing to correlate
p_sim normal while I is nan permutation p is computed regardless check np.isfinite(moran.I) explicitly

Islands do not cause nan. Cornwall with the Isles of Scilly as an island returned I = 0.1689 perfectly happily from the complete 2021 rates — the islands were simply given zero weight and excluded from the calculation without appearing anywhere in the result.

Two independent weights problems: islands, and non-finite values producing nan. two unrelated failures at the same line of codeSymptomCause and fix"not fully connected" warningislands — drop them, or use KNNweightsI = nan, p = 0.001one NaN or inf in the values —drop firstI = nan, all inputs finitezero variance — nothing tocorrelatean inland unit reported as anislandgeometry does not touch — clean itIslands do NOT cause nan — they are silently excluded spatialworkflow.io
Two unrelated failures that get confused because they surface at the same line of code.

The nan first, then zero variance, then islands

1. Find the nan first — it is the one that voids the result

values = zones["rate_2011"].to_numpy()
bad = ~np.isfinite(values)
if bad.any():
    print(f"{bad.sum()} non-finite values at index {zones.index[bad].tolist()[:5]}")
    print(zones.loc[bad, ["LSOA21NM", "count", "population", "rate_2011"]].head())
22 non-finite values at index ['E01034847', 'E01034848', 'E01034849', 'E01034850', 'E01034851']
                LSOA21NM  count  population  rate_2011
LSOA21CD                                              
E01034847  Cornwall 005F     13        1820        NaN
E01034848  Cornwall 005G     23        2136        NaN
E01034849  Cornwall 016E      6        1467        NaN
E01034850  Cornwall 016F     18        1399        NaN
E01034851  Cornwall 016G      3        1209        NaN

Twenty-two neighbourhoods whose 2021 codes did not exist in 2011: ONS split or merged them at the 2021 census, so the 2011 table has no row for them and the division produced NaN. Moran's I is computed from sums over all units, so one non-finite value poisons the entire statistic — nan propagates through every term.

Common sources:

  • a join on codes that changed between releases — as here, 2011 statistics on 2021 boundaries

  • a rate with a zero denominator — a zone with no population

  • a failed join leaving NaN in the value column

  • suppressed values in official statistics, often -99999 rather than NaN, which do not raise but do wreck the mean

2. Check that zero variance is not the cause

if np.nanvar(values) == 0:
    print("all values identical — Moran's I is undefined")

Moran's I divides by the variance of the values. If every unit has the same value there is no variation to correlate spatially, and the result is 0/0. This happens more often than expected — a filtered subset where every row has the same category, or a count column that is zero everywhere.

3. Then deal with the islands

print(f"{len(w.islands)} islands of {w.n} units")
print(f"{w.n_components} disconnected components")
print(zones.loc[w.islands, ["LSOA21NM"]].head())
1 islands of 337 units
2 disconnected components
                       LSOA21NM
LSOA21CD                       
E01019077  Isles of Scilly 001A

An island is a unit that shares no boundary with any other — a literal island, an exclave, or an artefact of geometry that does not quite touch. Contiguity weights give it no neighbours, so it contributes nothing to the statistic and is silently dropped.

That is not always wrong. But it means your reported n is larger than the number of units the statistic actually used, and nothing in the output says so.

4. Choose a deliberate response to the islands

from libpysal.weights import KNN, Queen

# a) drop them, and say so
connected = zones.drop(index=w.islands)
w_dropped = Queen.from_dataframe(connected, use_index=True)
w_dropped.transform = "r"
print(f"dropped {len(w.islands)}: {w_dropped.n} units, {len(w_dropped.islands)} islands")

# b) use k-nearest neighbours, which guarantees every unit has neighbours
w_knn = KNN.from_dataframe(zones, k=4)
w_knn.transform = "r"
print(f"KNN k=4: {w_knn.n} units, {len(w_knn.islands)} islands")
dropped 1: 336 units, 0 islands
KNN k=4: 337 units, 0 islands

Both are defensible, and they answer slightly different questions:

rate = zones["rate"].to_numpy(dtype="float64")
w.transform = "r"
for label, weights, frame in (("Queen with islands", w, zones), ("Queen, islands dropped", w_dropped, connected),
                              ("KNN k=4", w_knn, zones)):
    mi = Moran(frame["rate"].to_numpy(dtype="float64"), weights, permutations=999)
    print(f"  {label:<24} I = {mi.I:.4f}   (n = {weights.n - len(weights.islands)})")
  Queen with islands       I = 0.1689   (n = 336)
  Queen, islands dropped   I = 0.1687   (n = 336)
  KNN k=4                  I = 0.1859   (n = 337)

Dropping the island changes I in the fourth decimal place, which is the reassuring case: one unit out of 337 was never going to matter. When the treatments disagree, the islands were doing real work and the choice needs justifying.

5. Assert the result is finite

values = zones["rate"].to_numpy(dtype="float64")         # the 2021 rate, complete
moran = Moran(values, w, permutations=999)
if not np.isfinite(moran.I):
    raise ValueError("Moran's I is not finite — check for NaN, inf or zero variance")
print(f"I = {moran.I:+.4f}  p = {moran.p_sim:.4f}  n = {w.n - len(w.islands)}")
I = +0.1689  p = 0.0010  n = 336

Reporting n as the number of units with neighbours, not the number of rows, is what makes the result honest.

A connected grid plus two detached cells that receive zero weight and drop out of the statistic. zero weight, zero contribution, still counted in nwhat you builtIsles of Scillyn = 337 reportedwhat was analysedzero weightn = 336 usedReport w.n − len(w.islands), not the number of rows spatialworkflow.io
Islands do not error. They get zero weight, contribute nothing, and the reported n still counts them.

A strict weights builder on Cornwall and Scilly

The data used on this page

Every example runs on real data: the 337 Lower Layer Super Output Areas (LSOAs) of Cornwall and the Isles of Scilly from the ONS, their Census 2021 population from Nomis, and the street-level crimes Devon and Cornwall Police published for June 2026 on data.police.uk. count is crimes, rate is crimes per 1,000 residents.

import io
import time

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

np.random.seed(0)                     # every permutation test on this page draws from here

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Cornwall %' OR LSOA21NM LIKE 'Isles of Scilly %'", "outFields": "LSOA21CD,LSOA21NM",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).set_index("LSOA21CD").sort_index()

population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = population

west, south, east, north = zones.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in reply.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
incidents = gpd.sjoin(incidents, zones[["geometry"]], predicate="within").rename(columns={"index_right": "lsoa", "LSOA21CD": "lsoa"})
zones["count"] = incidents["lsoa"].value_counts().reindex(zones.index, fill_value=0)
zones["rate"] = 1000 * zones["count"] / zones["population"]
print(f"{len(zones)} LSOAs, {zones['population'].sum():,} residents, {len(incidents):,} crimes; "
      f"rate median {zones['rate'].median():.1f}, max {zones['rate'].max():.1f} per 1,000")
337 LSOAs, 572,373 residents, 3,577 crimes; rate median 4.7, max 56.1 per 1,000

Example 1 — a weights builder that refuses to hide anything

import geopandas as gpd
import numpy as np
import warnings
from esda.moran import Moran
from libpysal.weights import KNN, Queen, Rook


def build_weights(zones, scheme="queen", *, k=4, transform="r", islands="report"):
    """Build spatial weights and report exactly what happened to disconnected units."""
    builders = {
        "queen": lambda g: Queen.from_dataframe(g, use_index=True),
        "rook": lambda g: Rook.from_dataframe(g, use_index=True),
        "knn": lambda g: KNN.from_dataframe(g, k=k),
    }

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")          # we report this ourselves, in detail
        w = builders[scheme](zones)

    if w.islands and islands == "drop":
        kept = zones.drop(index=w.islands)
        print(f"dropped {len(w.islands)} island(s): {list(w.islands)[:5]}")
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            w = builders[scheme](kept)
        zones = kept
    elif w.islands:
        print(f"{len(w.islands)} island(s) will get zero weight and be excluded "
              f"from the statistic: {list(w.islands)[:5]}")

    w.transform = transform
    print(f"{scheme}: n={w.n}, {w.n_components} component(s), "
          f"mean neighbours {w.mean_neighbors:.2f}, "
          f"min {min(w.cardinalities.values())}, max {max(w.cardinalities.values())}")
    return w, zones


def safe_moran(zones, column, w, *, permutations=999):
    values = zones[column].to_numpy(dtype="float64")

    bad = ~np.isfinite(values)
    if bad.any():
        raise ValueError(
            f"{bad.sum()} non-finite value(s) in {column!r} at index "
            f"{zones.index[bad].tolist()[:5]} — Moran's I would return nan"
        )
    if np.var(values) == 0:
        raise ValueError(f"{column!r} has zero variance — Moran's I is undefined")

    moran = Moran(values, w, permutations=permutations)
    if not np.isfinite(moran.I):
        raise ValueError(f"Moran's I is {moran.I} despite finite inputs")

    effective = w.n - len(w.islands)
    print(f"I = {moran.I:+.4f}  E[I] = {moran.EI:+.4f}  z = {moran.z_sim:+.2f}  "
          f"p = {moran.p_sim:.4f}  (n = {effective})")
    return moran


w, zones_used = build_weights(zones, "queen")
try:
    moran = safe_moran(zones_used, "rate_2011", w)
except ValueError as error:
    print(f"ValueError: {error}")
1 island(s) will get zero weight and be excluded from the statistic: ['E01019077']
('WARNING: ', 'E01019077', ' is an island (no neighbors)')
queen: n=337, 2 component(s), mean neighbours 4.94, min 0, max 13
ValueError: 22 non-finite value(s) in 'rate_2011' at index ['E01034847', 'E01034848', 'E01034849', 'E01034850', 'E01034851'] — Moran's I would return nan

The exception names the column, the count and the offending index. Compare that with I = nan, p = 0.001 slipping through a pipeline unnoticed.

Example 2 — fixing the rate that produced the inf

def safe_rate(zones, count_col, denom_col, *, min_denom=1, per=1000):
    """Compute a rate, making suppressed and zero denominators explicit."""
    counts = zones[count_col].to_numpy(dtype="float64")
    denom = zones[denom_col].to_numpy(dtype="float64")

    # official statistics often use negative sentinels for suppressed values
    suppressed = (counts < 0) | (denom < 0)
    too_small = denom < min_denom

    rate = np.full(len(zones), np.nan)
    usable = ~suppressed & ~too_small
    rate[usable] = counts[usable] / denom[usable] * per

    print(f"{usable.sum()} usable, {suppressed.sum()} suppressed, "
          f"{too_small.sum()} with denominator < {min_denom}")
    return rate, usable


zones["rate"], usable = safe_rate(zones, "count", "population", min_denom=50)
analysis = zones[usable].copy()
print(f"analysing {len(analysis)} of {len(zones)} zones")
337 usable, 0 suppressed, 0 with denominator < 50
analysing 337 of 337 zones

min_denom is doing real work beyond avoiding the division by zero. A zone with three residents and one event has a rate of 333 per 1,000, which will be the maximum on your map and is entirely an artefact of the denominator. Excluding tiny denominators is the crude fix; Empirical Bayes smoothing via esda.moran.Moran_Rate is the proper one.

Example 3 — comparing weights schemes when islands exist

import pandas as pd


def islands_sensitivity(zones, column):
    rows = []
    for label, scheme, kwargs in [
        ("queen, islands kept", "queen", {"islands": "report"}),
        ("queen, islands dropped", "queen", {"islands": "drop"}),
        ("knn k=4", "knn", {"k": 4}),
        ("knn k=8", "knn", {"k": 8}),
    ]:
        w, used = build_weights(zones, scheme, **kwargs)
        moran = safe_moran(used, column, w)
        rows.append({
            "weights": label,
            "n_used": w.n - len(w.islands),
            "I": round(moran.I, 4),
            "p": moran.p_sim,
        })
    frame = pd.DataFrame(rows)
    print(f"\nI spread: {frame['I'].max() - frame['I'].min():.4f}")
    return frame


print(islands_sensitivity(zones, "rate").to_string(index=False))
1 island(s) will get zero weight and be excluded from the statistic: ['E01019077']
('WARNING: ', 'E01019077', ' is an island (no neighbors)')
queen: n=337, 2 component(s), mean neighbours 4.94, min 0, max 13
I = +0.1689  E[I] = -0.0030  z = +4.62  p = 0.0010  (n = 336)
dropped 1 island(s): ['E01019077']
queen: n=336, 1 component(s), mean neighbours 4.95, min 1, max 13
I = +0.1687  E[I] = -0.0030  z = +4.77  p = 0.0010  (n = 336)
knn: n=337, 2 component(s), mean neighbours 4.00, min 4, max 4
I = +0.1859  E[I] = -0.0030  z = +5.71  p = 0.0010  (n = 337)
knn: n=337, 1 component(s), mean neighbours 8.00, min 8, max 8
I = +0.0947  E[I] = -0.0030  z = +4.21  p = 0.0020  (n = 337)

I spread: 0.0912
               weights  n_used      I     p
   queen, islands kept     336 0.1689 0.001
queen, islands dropped     336 0.1687 0.001
               knn k=4     337 0.1859 0.001
               knn k=8     337 0.0947 0.002

The islands are not driving the result — keeping or dropping the Isles of Scilly moves I by 0.0002. The neighbourhood size is: eight nearest neighbours give an I of 0.095, about half the queen value, so the spread of 0.091 is a statement about scale. A spread that size is the finding, and it needs reporting in the text rather than burying in a default.

Why nan arrives with a plausible p-value

Why nan comes with a plausible p-value

Moran computes the statistic and the permutation p-value through separate paths. The observed I is a ratio of sums over the values; one nan makes every sum nan.

The permutation p-value is the proportion of shuffled arrangements whose I was at least as extreme as the observed one. Comparisons against nan are always False in NumPy, so no permutation ever "beats" the observed value, and the count of exceedances stays at zero — which produces the smallest possible p-value, 1/(permutations+1) = 0.001.

So a broken statistic yields the most significant-looking p-value available. This is why np.isfinite(moran.I) belongs in the code and not in a comment.

Why islands are silently excluded rather than an error

Zero neighbours means zero weight in every term. Mathematically the unit simply drops out; there is no division by zero and nothing to raise about.

libpysal warns because the situation is usually unintended, but the warning is easy to lose in a notebook and impossible to see in a scheduled job that discards stderr. The number to report is w.n - len(w.islands) — the units that actually contributed.

The deeper issue is what an island means for your model. A genuine island has no land neighbours, so a contiguity model says nothing influences it. If your process is economic or social rather than physical, that is probably wrong, and KNN or a distance band is a better model of influence than shared boundaries.

Why disconnected components matter beyond islands

1 islands of 337 units
2 disconnected components

Two components, one of which is the Isles of Scilly — so mainland Cornwall is one connected block. But a study area split into two large disconnected halves (an archipelago, a region separated by an international boundary you excluded) is a different situation: the statistic then measures autocorrelation within each half with no comparison between them.

Check w.n_components. More than one plus the island count means you have genuinely separate sub-regions, and a single global statistic across them is answering an odd question.

A single NaN making Moran’s I nan while the permutation p-value comes out at its minimum. a broken statistic with the best-looking p-value availableone failed joinrate = NaNevery sum becomesnannan comparisonsare all Falseno permutation "beats"itp = 1/(999+1) =0.001the minimumCheck np.isfinite(moran.I) — the p-value will not tell you spatialworkflow.io
Comparisons against nan are always false, so no permutation beats the observed value and p comes out at its minimum.

Why "not fully connected" is sometimes a geometry bug

Not every island is a real island. Two polygons that should share a boundary may fail to touch because of:

  • precision — coordinates that differ in the eighth decimal place
  • slivers — a hairline gap between neighbours from different sources
  • invalid geometry — a self-intersection that breaks the touch test

If a unit you know is inland comes back as an island, check the geometry before changing the weights scheme:

suspect = zones.loc[w.islands]
print(suspect.geometry.is_valid)
neighbours = zones[zones.geometry.buffer(1).intersects(suspect.geometry.iloc[0])]
print(f"{len(neighbours) - 1} units within 1 m but not touching")
LSOA21CD
E01019077    True
dtype: bool
0 units within 1 m but not touching

Here the island's geometry is valid and nothing lies within 1 m of it: E01019077 really is detached, so the warning is correct. A unit with neighbours 1 m away but no shared boundary is a topology problem, not a spatial-weights problem, and the fix is to snap or clean the geometry.

Index labels, asymmetric KNN and sentinels

  • w.islands is a list of index labels, not positions. Use .drop(index=...), not .drop(...) with integers, unless your index is positional.
  • use_index=True is required by current libpysal when building from a GeoDataFrame, and keeps the weights aligned with your index.
  • KNN weights are asymmetric. A being B's nearest neighbour does not make B A's. That is fine for Moran's I, which symmetrises internally, but surprising if you inspect the matrix.
  • Row-standardising weights with islands leaves those rows summing to zero, not one — another reason they contribute nothing.
  • Moran_Rate applies an Empirical Bayes correction and is the right tool for rates with varying denominators.
  • Zero variance also breaks Geary's C and Getis-Ord. The check is worth applying once, before any spatial statistic.
  • A -99999 sentinel is finite, so np.isfinite will not catch it. Check the minimum of every numeric column from an official source.
  • w.n_components > 1 + len(w.islands) means genuinely separate sub-regions, which is worth handling explicitly.

FAQ

Does an island cause Moran's I to return nan?

No. Islands get zero weight and are silently excluded. nan comes from a non-finite value in the data or from zero variance.

Why is the p-value significant when I is nan?

Comparisons against nan are always false, so no permutation exceeds the observed value and the p-value comes out at its minimum, 1/(permutations+1). Always check np.isfinite(moran.I).

Should I drop islands or use KNN weights?

Either, deliberately. Drop them when a unit with no land neighbours genuinely has no spatial relationship; use KNN when influence travels by distance rather than shared boundaries. Report which you chose.

What causes a inf in a rate column?

A zero denominator. Guard with a minimum population, and remember that very small denominators produce extreme rates that are artefacts rather than findings.

A unit I know is inland is reported as an island. Why?

Its geometry does not quite touch its neighbours — precision, a sliver, or an invalidity. Fix the geometry rather than the weights.

What does "3 disconnected components" mean?

The adjacency graph is in three separate pieces. If that is more than the number of islands, you have genuinely separate sub-regions, and one global statistic across them may not be meaningful.

How do I report n honestly?

As w.n - len(w.islands) — the units that actually contributed to the statistic, not the number of rows in the frame.