Spatial autocorrelation: Peninsular Malaysia self-study companion

# %pip install numpy pandas matplotlib
from pathlib import Path
import io, json, hashlib, urllib.request, urllib.parse
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import font_manager
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib.collections import PatchCollection
from matplotlib.patches import PathPatch, Patch
from matplotlib.path import Path as MplPath
from IPython.display import display

CACHE = Path('peninsula_data_cache')
CACHE.mkdir(exist_ok=True)
REFRESH = False
USE_BOOK_SNAPSHOT = True
# The checking script can supply BOOK_DATA_DIR for an offline run.
book_data_dir = globals().get('BOOK_DATA_DIR')
if book_data_dir is None:
    for parent in [Path.cwd(), *Path.cwd().parents]:
        for candidate in [parent/'src/content/datasets/peninsula_income_2022', parent/'datasets/peninsula_income_2022']:
            if (candidate/'snapshot.json').exists():
                book_data_dir = candidate
                break
        if book_data_dir is not None: break
INCOME_URL = 'https://storage.dosm.gov.my/hies/hh_income_district.csv'
BOUNDARY_SERVICE = 'https://gisdev.planmalaysia.gov.my/server/rest/services/Hosted/Data_Sempadan_Daerah/FeatureServer/0'
params = {'where': "kod_negeri NOT IN ('12','13','15')", 'outFields': 'negeri,kod_negeri,daerah,kod_daerah',
          'returnGeometry': 'true', 'outSR': '3375', 'maxAllowableOffset': '100', 'geometryPrecision': '1', 'f': 'json'}
BOUNDARY_URL = BOUNDARY_SERVICE + '/query?' + urllib.parse.urlencode(params)
def cached_bytes(url, filename):
    p = CACHE / filename
    if REFRESH or not p.exists():
        with urllib.request.urlopen(url, timeout=60) as response:
            payload = response.read()
        p.write_bytes(payload)
    return p.read_bytes()
if USE_BOOK_SNAPSHOT and book_data_dir is not None:
    import gzip
    book_data_dir = Path(book_data_dir)
    income_bytes = (book_data_dir/'hh_income_district.csv').read_bytes()
    boundary_bytes = gzip.decompress((book_data_dir/'boundaries_rso_display.json.gz').read_bytes())
    expected = json.loads((book_data_dir/'snapshot.json').read_text())
    assert hashlib.sha256(income_bytes).hexdigest() == expected['income_sha256']
    assert hashlib.sha256(boundary_bytes).hexdigest() == expected['boundaries_sha256']
    data_route = 'book snapshot, retrieved 2026-10-20'
else:
    income_bytes = cached_bytes(INCOME_URL, 'hh_income_district.csv')
    boundary_bytes = cached_bytes(BOUNDARY_URL, 'peninsula_boundaries_rso_display.json')
    data_route = 'live sources (or local download cache)'
print('Data route:', data_route)
all_income = pd.read_csv(io.BytesIO(income_bytes))
layer = json.loads(boundary_bytes)
if 'error' in layer:
    raise RuntimeError(layer['error'])
assert not layer.get('exceededTransferLimit', False), 'Boundary query was truncated.'
assert layer['spatialReference']['wkid'] == 3375, 'Unexpected boundary CRS.'
features = layer['features']
geo_table = pd.DataFrame([f['attributes'] for f in features])
assert len(features) == 92 and geo_table['kod_daerah'].is_unique
print('Loaded 92 mapped areas; display CRS EPSG:3375 (metres).')
source_hashes = {'income_sha256': hashlib.sha256(income_bytes).hexdigest(),
                 'boundaries_sha256': hashlib.sha256(boundary_bytes).hexdigest()}
print(source_hashes)
display(geo_table[['negeri','daerah','kod_daerah']])

region_income = all_income[~all_income['state'].isin(['Sabah','Sarawak','W.P. Labuan'])].copy()
region_income['year'] = pd.to_datetime(region_income['date']).dt.year
print('Variables: income_median = median; income_mean = mean gross monthly household income (RM).')
display(region_income.groupby(['year','state']).size().unstack(fill_value=0).rename_axis('year'))
display(region_income.sort_values(['year','state','district']))

YEAR = 2022
VARIABLE = 'income_median'         # Or 'income_mean'.
CMAP = 'OrRd'                      # Sequential cream/orange/dark red; try 'YlOrBr' for gold.
N_CLASSES = 5
CLASS_BREAKS = [0, 3000, 5000, 7000, 9000, 15000]  # Or None for quantiles.
SEED = 5230
USE_XKCD = True
LABEL_DISTRICTS = False            # Full labels can be crowded in the small panels.
EXPORT = True
OUTPUT_DIR = Path('peninsula_figures')

if VARIABLE not in ['income_mean','income_median']:
    raise ValueError('Select income_mean or income_median.')
selected = region_income[region_income['year'] == YEAR].copy()
name_aliases = {'Ulu Langat':'HULU LANGAT', 'Ulu Selangor':'HULU SELANGOR',
                'W.P. Kuala Lumpur':'KUALA LUMPUR', 'W.P. Putrajaya':'PUTRAJAYA',
                'Perlis':'KANGAR', 'Kecil Lojing':'LOJING', 'Larut dan Matang':'LARUT & MATANG', 'Timur Laut':'TIMOR LAUT'}
selected['join_district'] = selected['district'].map(lambda x:name_aliases.get(x,x.upper()))
selected['join_state'] = selected['state'].str.upper()
geo_table['join_state'] = geo_table['negeri'].str.upper()
geo_table['join_district'] = geo_table['daerah'].str.upper()
if len(selected) != 92:
    raise ValueError(f'Missing {YEAR} income: expected 92 areas; choose a complete year such as 2022.')
df = geo_table.merge(selected,on=['join_state','join_district'],how='left',validate='one_to_one',sort=False)
assert df['kod_daerah'].tolist() == geo_table['kod_daerah'].tolist(), 'Feature order changed.'
if df[VARIABLE].isna().any():
    missing = df.loc[df[VARIABLE].isna(),'daerah'].tolist()
    raise ValueError(f'Missing {YEAR} income for {missing}. Choose a complete common year, such as 2022.')
values = df[VARIABLE].to_numpy(dtype=float)
assert np.isfinite(values).all() and (values>0).all()
if CLASS_BREAKS is None:
    if not isinstance(N_CLASSES,int) or N_CLASSES<2:
        raise ValueError('N_CLASSES must be an integer of at least 2.')
    breaks = np.unique(np.quantile(values,np.linspace(0,1,N_CLASSES+1)))
else:
    breaks = np.asarray(CLASS_BREAKS,dtype=float)
if len(breaks)<3 or not np.isfinite(breaks).all() or not (np.diff(breaks)>0).all():
    raise ValueError('Need at least three increasing, distinct finite breaks.')
if breaks[0]>values.min() or breaks[-1]<values.max():
    raise ValueError('Class breaks do not cover all values.')
cmap = ListedColormap(plt.get_cmap(CMAP)(np.linspace(.05,.95,len(breaks)-1)))
norm = BoundaryNorm(breaks,cmap.N,clip=True)
rng = np.random.default_rng(SEED)
permutations = [rng.permutation(values) for _ in range(4)]
assert all(np.array_equal(np.sort(p),np.sort(values)) for p in permutations)
title = f'{"Median" if VARIABLE=="income_median" else "Mean"} gross monthly household income, {YEAR} (RM)'
assert len(df) == 92 and len(set(df['kod_daerah'])) == 92
print(title)
print('Class breaks:', breaks)
display(df[['daerah','income_median','income_mean']].sort_values(VARIABLE,ascending=False))

handwriting_font = 'DejaVu Sans'
if USE_XKCD:
    try:
        font = globals().get('BOOK_FONT_PATH')
        if font is None and book_data_dir is not None:
            candidate = Path(book_data_dir).parents[3]/'_static/css/fonts/ComicNeue-Regular.ttf'
            if candidate.exists(): font = candidate
        font = Path(font) if font is not None else CACHE/'ComicNeue-Regular.ttf' 
        if not font.exists():
            with urllib.request.urlopen('https://raw.githubusercontent.com/google/fonts/main/ofl/comicneue/ComicNeue-Regular.ttf',timeout=30) as response:
                font.write_bytes(response.read())
        font_manager.fontManager.addfont(font)
        handwriting_font = 'Comic Neue'
    except Exception as exc:
        print('Using DejaVu Sans:',type(exc).__name__)
patches=[];centres=[]
for feature in features:
    vertices=[];codes=[]
    rings=feature['geometry']['rings']
    for ring in rings:
        a=np.asarray(ring,dtype=float)[:,:2]
        assert len(a)>=4 and np.allclose(a[0],a[-1])
        vertices.extend(a)
        codes.extend([MplPath.MOVETO]+[MplPath.LINETO]*(len(a)-2)+[MplPath.CLOSEPOLY])
    patch=PathPatch(MplPath(vertices,codes))
    patches.append(patch)
    a=np.asarray(vertices);centres.append((a.min(axis=0)+a.max(axis=0))/2)
vertices=np.concatenate([p.get_path().vertices for p in patches])
lo,hi=vertices.min(axis=0),vertices.max(axis=0)
assert layer['spatialReference']['wkid'] == 3375, 'Expected Peninsula RSO coordinates.'

def map_panel(ax,data,label):
    collection=PatchCollection(patches,facecolors=[cmap(norm(v)) for v in data],
                               edgecolors='#444444',linewidths=.6,rasterized=True)
    collection.set_sketch_params(None);collection.set_path_effects([])
    ax.add_collection(collection)
    margin=(hi-lo)*.04
    ax.set_xlim(lo[0]-margin[0],hi[0]+margin[0]);ax.set_ylim(lo[1]-margin[1],hi[1]+margin[1])
    ax.set_aspect('equal');ax.axis('off');ax.set_title(label,fontsize=15)
    ax.text(.95,.96,'N',transform=ax.transAxes,ha='center')
    ax.annotate('',xy=(.95,.93),xytext=(.95,.83),xycoords='axes fraction',arrowprops={'arrowstyle':'->'})
    for i,(x,y) in enumerate(centres):
        district=df.iloc[i]['daerah']
        if LABEL_DISTRICTS:
            ax.text(x,y,district.title(),ha='center',fontsize=8)

def draw_figure(panels,name):
    with plt.rc_context():
        if USE_XKCD:plt.xkcd(scale=.5,length=100,randomness=1)
        plt.rcParams.update({'font.family':handwriting_font if USE_XKCD else 'DejaVu Sans','svg.fonttype':'path'})
        np.random.seed(SEED)
        count=len(panels)
        width = {1:6.0, 2:10.5, 3:14.5}[count]
        map_ratio = float((hi[0]-lo[0])/(hi[1]-lo[1]))
        map_height = .92*width/(count+(count-1)*.045)/map_ratio
        height = map_height+1.35  # Fixed room for panel headings and the shared legend.
        fig,axes=plt.subplots(1,count,figsize=(width,height),squeeze=False)
        for ax,(data,label) in zip(axes[0],panels):map_panel(ax,data,label)
        labels=[f'{a:,.0f} to < {b:,.0f}' for a,b in zip(breaks[:-1],breaks[1:])]
        labels[0]=f'Below {breaks[1]:,.0f}'
        labels[-1]=f'{breaks[-2]:,.0f}+'
        handles=[Patch(facecolor=cmap(i),edgecolor='#444444',label=l) for i,l in enumerate(labels)]
        fig.legend(handles=handles,loc='lower center',bbox_to_anchor=(.5,.035),ncol=3,frameon=False,fontsize=9,title=title,title_fontsize=10)
        fig.text(.5,.012,'Income: DOSM · Display boundaries: PLANMalaysia · EPSG:3375',ha='center',fontsize=8)
        fig.subplots_adjust(left=.035,right=.985,bottom=1.1/height,top=1-.25/height,wspace=.045)
        if EXPORT:
            OUTPUT_DIR.mkdir(exist_ok=True)
            for ext in ['png','svg']:fig.savefig(OUTPUT_DIR/f'{YEAR}_{VARIABLE}_{name}.{ext}',dpi=200,bbox_inches='tight')
            print('Exported:',OUTPUT_DIR/f'{YEAR}_{VARIABLE}_{name}.png','and .svg')
        plt.show()
    return fig

fig1=draw_figure([(values,'Observed values')],'01_observed')

fig2=draw_figure([(values,'Observed values'),(permutations[0],'One random permutation')],'02_observed_vs_permuted')

fig3=draw_figure([(p,f'Permutation {i}') for i,p in enumerate(permutations[1:],1)],'03_three_permutations')

settings_record={'year':YEAR,'variable':VARIABLE,'class_breaks':breaks.tolist(),'palette':CMAP,
                 'seed':SEED,'xkcd':USE_XKCD,'crs':'EPSG:3375','areas':92,'display_generalization_metres':100,
                 'data_route':data_route,'income_url':INCOME_URL,'boundary_query_url':BOUNDARY_URL,**source_hashes}
print(json.dumps(settings_record,indent=2))
if EXPORT:(OUTPUT_DIR/f'{YEAR}_{VARIABLE}_settings.json').write_text(json.dumps(settings_record,indent=2))

chain_edges = [(0,1),(1,2),(2,3)]
def neighbour_products(x):
    z = np.asarray(x,dtype=float)-np.mean(x)
    products = [z[i]*z[j] for i,j in chain_edges]
    return products, sum(products)
for example in [[2,2,8,8],[2,8,2,8]]:
    print(example, neighbour_products(example))
assert neighbour_products([2,2,8,8]) == ([9,-9,9],9)
assert neighbour_products([2,8,2,8]) == ([-9,-9,-9],-27)