التحليل المكاني GIS للمؤسسات: خرائط وإحصاء مكاني بلغتي R وPython

import numpy as np

# Synthetic 10 x 10 grid of areas with a spatial trend
rng = np.random.default_rng(1)
i, j = np.indices((10, 10))
x = (i + j + rng.normal(0, 3, (10, 10))).ravel()

# Rook contiguity weights, row-standardised
idx = np.arange(100).reshape(10, 10)
W = np.zeros((100, 100))
pairs = [(idx[:-1], idx[1:]), (idx[:, :-1], idx[:, 1:])]
for a, b in pairs:
    W[a.ravel(), b.ravel()] = W[b.ravel(), a.ravel()] = 1
W /= W.sum(axis=1, keepdims=True)

# Global Moran's I with a permutation test
z = x - x.mean()
moran = z @ W @ z / (z @ z)
sims = np.array([p @ W @ p / (p @ p) for p in
                 (rng.permutation(z) for _ in range(999))])
p_val = (np.sum(sims >= moran) + 1) / 1000
print(f"Moran's I = {moran:.3f}, pseudo p = {p_val:.3f}")

library(sf)
library(spdep)

nc <- st_read(system.file("shape/nc.shp", package = "sf"),
              quiet = TRUE)
nc$rate <- 1000 * nc$SID79 / nc$BIR79

nb <- poly2nb(nc, queen = TRUE)
lw <- nb2listw(nb, style = "W")

# Global and local spatial autocorrelation
moran.test(nc$rate, lw)
nc$lisa <- localmoran(nc$rate, lw)[, "Ii"]

# Getis-Ord Gi* hot spots (neighbours include self)
lw_self <- nb2listw(include.self(nb), style = "W")
nc$gi_star <- as.numeric(localG(nc$rate, lw_self))

plot(nc["gi_star"], main = "Gi* z-scores")