Machine learning with spatial data

"""A made-up land-cover scene of 90 by 60 pixels of 20 m (1.8 km by 1.2 km), from a seeded generator: patches of crop,
forest and built-up land, and a river. Each pixel has red, near-infrared and shortwave-infrared reflectance: the
class's own values, a shift shared by its whole patch (soil, crop variety, roof material), a gentle haze across the
scene, and noise. Then 1,200 reference pixels whose true class is known, as from a field survey."""
import math
import random

COLS, ROWS, PIXEL = 90, 60, 20
CLASSES = ["water", "crop", "forest", "built-up"]
MEAN = {  # reflectance in the red, near-infrared and shortwave-infrared bands
    "water": (0.04, 0.03, 0.02),
    "crop": (0.06, 0.34, 0.20),
    "forest": (0.04, 0.27, 0.13),
    "built-up": (0.15, 0.22, 0.24),
}


def river_y(x):
    """The river's centre line, in pixels, at column x."""
    return 34 + 9 * math.sin(x / 14)


def build(seed=21):
    """(classes, bands, patch) for every pixel, row by row; bands are (red, nir, swir)."""
    rng = random.Random(seed)
    centres = [(rng.uniform(0, COLS), rng.uniform(0, ROWS), rng.choices(CLASSES[1:], weights=(5, 3, 2))[0]) for _ in range(45)]
    shift = [tuple(rng.gauss(0, 0.025) for _ in range(3)) for _ in centres]
    classes, bands, patch = [], [], []
    for r in range(ROWS):
        for c in range(COLS):
            if abs(r - river_y(c)) < 2.2:
                k, p = "water", -1
            else:
                p = min(range(len(centres)), key=lambda i: (c - centres[i][0]) ** 2 + (r - centres[i][1]) ** 2)
                k = centres[p][2]
            haze = 0.02 * c / COLS
            offset = shift[p] if p >= 0 else (0, 0, 0)
            pixel = tuple(max(0.0, MEAN[k][b] + offset[b] + haze + rng.gauss(0, 0.02)) for b in range(3))
            classes.append(k)
            bands.append(pixel)
            patch.append(p)
    return classes, bands, patch


def reference(n=1_200, seed=5):
    """Pixel numbers of n reference pixels drawn at random."""
    return random.Random(seed).sample(range(COLS * ROWS), n)

"""Features of the scene's pixels, a k-nearest-neighbour classifier, and two ways to split reference pixels into five
folds: at random, and by spatial blocks, so that a whole block is either trained on or tested."""
import heapq
import math
import random

from scene import COLS, PIXEL, ROWS, river_y


def ndvi(band):
    red, nir, _ = band
    return (nir - red) / (nir + red) if nir + red > 0 else 0.0


def table(bands):
    """For every pixel: its bands, its NDVI, the mean NDVI of its 3 x 3 neighbourhood, its distance to the river in
    metres, and its x and y in metres."""
    v = [ndvi(b) for b in bands]
    rows = []
    for i, b in enumerate(bands):
        r, c = divmod(i, COLS)
        near = [v[(r + dr) * COLS + c + dc] for dr in (-1, 0, 1) for dc in (-1, 0, 1) if 0 <= r + dr < ROWS and 0 <= c + dc < COLS]
        rows.append({"red": b[0], "nir": b[1], "swir": b[2], "ndvi": v[i], "ndvi 3x3": sum(near) / len(near),
                     "river m": abs(r - river_y(c)) * PIXEL, "x": c * PIXEL, "y": r * PIXEL})
    return rows


def scaled(rows, names):
    """The chosen features of each row, each scaled to mean 0 and standard deviation 1 over the rows."""
    cols = [[row[n] for row in rows] for n in names]
    stats = [(sum(col) / len(col), math.sqrt(sum((v - sum(col) / len(col)) ** 2 for v in col) / len(col)) or 1.0) for col in cols]
    return [tuple((row[n] - m) / s for n, (m, s) in zip(names, stats)) for row in rows]


def knn(train, labels, test, k=5):
    """The majority class of the k nearest training rows of each test row (ties: the class first in alphabetical order)."""
    pairs = list(zip(train, labels))
    out = []
    for x in test:
        votes = {}
        for _, label in heapq.nsmallest(k, ((math.dist(x, t), lab) for t, lab in pairs)):
            votes[label] = votes.get(label, 0) + 1
        out.append(max(sorted(votes), key=lambda lab: votes[lab]))
    return out


def random_folds(n, k=5, seed=1):
    order = list(range(n))
    random.Random(seed).shuffle(order)
    return [order[f::k] for f in range(k)]


def block_folds(pixels, k=5, size=15, seed=1):
    """Folds of whole blocks of size x size pixels: every block goes to one fold, the blocks dealt out at random."""
    blocks = {}
    for j, p in enumerate(pixels):
        r, c = divmod(p, COLS)
        blocks.setdefault((c // size, r // size), []).append(j)
    keys = sorted(blocks)
    random.Random(seed).shuffle(keys)
    folds = [[] for _ in range(k)]
    for n, key in enumerate(keys):
        folds[n % k] += blocks[key]
    return folds


def cross_validate(x, y, folds, k=5):
    """The predicted class of every sample, each from a model trained on the other folds."""
    predicted = [None] * len(y)
    for fold in folds:
        held = set(fold)
        train = [j for j in range(len(y)) if j not in held]
        for j, p in zip(fold, knn([x[t] for t in train], [y[t] for t in train], [x[j] for j in fold], k)):
            predicted[j] = p
    return predicted

"""A classifier that knows only where a pixel is: k-nearest neighbours on x and y alone, tested by random five-fold
cross-validation and by five folds of spatial blocks 300 m wide. The pixels' bands are not used at all."""
from features import block_folds, cross_validate, random_folds, scaled, table
from scene import CLASSES, build, reference

classes, bands, _ = build()
rows = table(bands)
pixels = reference()
x = scaled([rows[p] for p in pixels], ["x", "y"])
y = [classes[p] for p in pixels]
print(f"{len(pixels):,} reference pixels: " + ", ".join(f"{c} {y.count(c)}" for c in CLASSES))
print(f"  always the commonest class, crop: accuracy {y.count('crop') / len(y):.3f}")
for name, folds in (("random folds", random_folds(len(y))), ("blocks of 300 m", block_folds(pixels))):
    predicted = cross_validate(x, y, folds)
    print(f"  {name:<16} accuracy {sum(p == t for p, t in zip(predicted, y)) / len(y):.3f}")

1,200 reference pixels: water 103, crop 593, forest 277, built-up 227
  always the commonest class, crop: accuracy 0.494
  random folds     accuracy 0.891
  blocks of 300 m  accuracy 0.561