Module 9 · Applications · Free & interactive

Real Data: India, Live

Run spatial statistics on India's 725 districts, pulled live from the free geo-data.io API. Written by Eric Vaz. Read it here, or open the interactive version to write and run the Python yourself, in your browser, on real open data.

Open the interactive version (write & run the code) →

Spatial autocorrelation across India (live API)

Everything so far used teaching datasets. Now point the same methods at real data: the Census of India at the district level, pulled live from geo-data.io, a free, CORS-enabled open-geodata API (81 variables across 725 districts). This whole lesson runs on it, powered by geo-data.io.

Every variable is one URL: https://geo-data.io/v1/india/<key>.json. We'll map and test literacy rate, but you can swap in any of the 81 keys and re-run.

ExamplePull a variable from the geo-data.io API

from pyodide.http import open_url
import json

VAR = "literacy_rate"     # try: density, pm25, sex_ratio, pct_internet, idx_affluence ...
api = json.load(open_url(f"https://geo-data.io/v1/india/{VAR}.json"))
print(api["label"], "|", api["unit"], "|", api["year"], "|", api["count"], "districts")

vals = {r["dt_code"]: r["value"] for r in api["data"] if r["dt_code"] is not None}
print("\nHighest:")
for r in api["data"][:5]:
    print(f"  {r['district']:<22} {r['state']:<16} {r['value']}")

ExampleGet the geometry, join the values by dt_code

import numpy as np

gj = json.load(open_url("https://geo-data.io/india_districts.geojson"))

def feature_polys(f):
    g = f["geometry"]
    return g["coordinates"] if g["type"] == "MultiPolygon" else [g["coordinates"]]

def centroid(f):
    best, best_area = None, -1
    for poly in feature_polys(f):
        r = np.array(poly[0])
        a = abs(0.5*np.sum(r[:-1,0]*r[1:,1] - r[1:,0]*r[:-1,1]))   # shoelace
        if a > best_area:
            best_area, best = a, (r[:,0].mean(), r[:,1].mean())
    return best

feats = [f for f in gj["features"] if f["properties"].get("dt_code") in vals]
y     = np.array([vals[f["properties"]["dt_code"]] for f in feats], float)
cent  = np.array([centroid(f) for f in feats])
print(len(feats), "districts joined to geometry")
print("value range:", round(float(y.min()),1), "to", round(float(y.max()),1))

Try it yourselfNeighbours + Moran's I on India

Build k-nearest weights from the district centroids and test for spatial autocorrelation. Fill two blanks: row-standardise the weights, and feed a reshuffled copy of the values into the permutation test.

from scipy.spatial import cKDTree

def knn_W(coords, k=6):
    n = len(coords); tree = cKDTree(coords)
    _, idx = tree.query(coords, k=k+1)         # k+1: nearest point is itself
    W = np.zeros((n, n))
    for i in range(n):
        for j in idx[i, 1:]: W[i, j] = 1.0
    return W / W.sum(1, keepdims=True)

W = knn_W(cent, 6)

def morans_I(y, W):
    z = y - y.mean(); S0 = W.sum()
    return (len(y)/S0) * np.sum(W * np.outer(z, z)) / np.sum(z**2)

I = morans_I(y, W)
rng = np.random.default_rng(0)
ref = np.array([morans_I(rng.permutation(y), W) for _ in range(199)])
p = (np.sum(ref >= I) + 1) / 200
print(f"Moran's I ({api['label']}) = {I:.3f}   pseudo p = {p:.3f}")

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

ExampleMap India, and the Moran scatterplot

import matplotlib.pyplot as plt, matplotlib as mpl
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib.collections import PatchCollection

norm = mpl.colors.Normalize(y.min(), y.max()); cm = plt.get_cmap("viridis")
fig, ax = plt.subplots(1, 2, figsize=(10, 4.6))
patches, colors = [], []
for f, v in zip(feats, y):
    for poly in feature_polys(f):
        patches.append(PathPatch(Path(np.array(poly[0])))); colors.append(cm(norm(v)))
ax[0].add_collection(PatchCollection(patches, facecolor=colors, edgecolor="white", linewidth=.15))
ax[0].autoscale(); ax[0].set_aspect("equal"); ax[0].axis("off")
ax[0].set_title(api["label"] + " by district")
fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cm), ax=ax[0], shrink=.6)

zy = (y - y.mean())/y.std(); zlag = W @ zy
ax[1].axhline(0, color="#aaa", lw=.8); ax[1].axvline(0, color="#aaa", lw=.8)
ax[1].scatter(zy, zlag, s=10, color="#14233D")
slope = np.polyfit(zy, zlag, 1)[0]; xs = np.linspace(zy.min(), zy.max(), 2)
ax[1].plot(xs, slope*xs, color="#D6552F", lw=2, label=f"Moran's I = {slope:.3f}")
ax[1].set(xlabel="value (standardised)", ylabel="spatial lag", title="Moran scatterplot"); ax[1].legend()
plt.tight_layout()
print("Each district coloured by the variable; the scatterplot slope is Moran's I.")

This lesson runs on the free geo-data.io API

geo-data.io — India API catalogue — 81 variables, 725 districts; swap VAR for any key above and re-run
geo-data.io — free, CORS-enabled open geodata API

Found this useful?It's free and always will be. A coffee keeps new modules coming.

☕ Buy me a coffee

More free modules

Spatial Autocorrelation — Understand how spatial relationships create patterns.Interpolation & Kriging — Predict values at unsampled locations.Hot Spot Analysis — Identify statistically significant hot and cold spots.All modules & the interactive course →