Maps for Areal Data

https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_{COUNTRY}_{LEVEL}.json

import geopandas as gpd

url   = "https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_SAU_1.json"
saudi = gpd.read_file(url)

print(type(saudi).__name__)
print("rows:", len(saudi))

# Coordinate reference system: GADM always ships lat/lon (EPSG:4326).
print(saudi.crs)

# Identifier + metadata columns alongside `geometry`.
print(saudi.columns.tolist())

print(saudi[['GID_1', 'NAME_1', 'ENGTYPE_1']].head())

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 8))
saudi.plot(ax=ax, edgecolor='black', linewidth=0.5, facecolor='#e8f4f8')

for _, row in saudi.iterrows():
    c = row.geometry.centroid
    ax.annotate(row['NAME_1'], (c.x, c.y),
                ha='center', va='center', fontsize=7,
                bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.75))

ax.set_title('Saudi Arabia: 13 administrative regions (ADM1)', fontweight='bold')
ax.axis('off')
plt.tight_layout(); plt.show()

def get_admin_boundaries(country_code, admin_level=1):
    """ISO 3166-1 alpha-3 country code (e.g. 'SAU'); admin_level in 0..3."""
    url = f"https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_{country_code}_{admin_level}.json"
    return gpd.read_file(url)

saudi  = get_admin_boundaries('SAU', 1)
lebanon = get_admin_boundaries('LBN', 2)  # district-level

13                          // line 1: total number of regions (n)
1 2 4 5                     // region 1 has 2 neighbors: regions 4 and 5
2 3 5 6 12                  // region 2 has 3 neighbors: regions 5, 6, and 12
3 2 6 7
4 3 1 5 10
5 5 1 2 4 6 10
6 5 2 3 5 7 10
7 4 3 6 8 10
8 3 7 9 10
9 2 8 10
10 6 4 5 6 7 8 9
11 2 12 13
12 3 2 11 13
13 2 11 12

from libpysal.weights import Queen

def build_adjacency(gdf):
    """Build adjacency graph and return INLA-compatible graph string.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Administrative boundaries (row order defines region IDs)

    Returns
    -------
    dict with keys:
        'weights': libpysal weights object
        'n': number of regions
        'graph_str': INLA graph file content (string)
        'neighbors': dict mapping region index to neighbor indices
    """
    # Create Queen contiguity (regions sharing any boundary point are neighbors)
    w = Queen.from_dataframe(gdf)

    # Check for islands (regions with no neighbors)
    islands = [i for i, neighs in w.neighbors.items() if len(neighs) == 0]
    if islands:
        print(f"Warning: {len(islands)} region(s) with no neighbors: {islands}")

    # Build INLA graph format (1-indexed)
    n = len(gdf)
    lines = [str(n)]  # First line: number of regions

    for i in range(n):
        # Get neighbors and convert to 1-indexed, sorted
        neighbors_1idx = sorted([j + 1 for j in w.neighbors[i]])
        n_neighbors = len(neighbors_1idx)
        neighbor_str = ' '.join(map(str, neighbors_1idx))
        lines.append(f"{i + 1} {n_neighbors} {neighbor_str}")

    return {
        'weights': w,
        'n': n,
        'graph_str': '\n'.join(lines),
        'neighbors': w.neighbors,
        'mean_neighbors': w.mean_neighbors
    }

# Download boundaries
saudi = get_admin_boundaries('SAU', 1)

# Build adjacency
adj = build_adjacency(saudi)
print(f"Number of regions: {adj['n']}")
print(f"Average neighbors per region: {adj['mean_neighbors']:.1f}")

# Save graph file
with open('saudi.graph', 'w') as f:
    f.write(adj['graph_str'])

print(f"\nGraph file saved to: saudi.graph")
print(f"\nGraph file content:")
print(adj['graph_str'])

import matplotlib.pyplot as plt

w         = adj['weights']
centroids = saudi.geometry.centroid          # one point per region

fig, ax = plt.subplots(figsize=(10, 8))
saudi.plot(ax=ax, facecolor='#e8f4f8', edgecolor='#333', linewidth=0.5)

# Draw one segment per (i, j) neighbor pair.
for i, neighbours in w.neighbors.items():
    xi, yi = centroids.iloc[i].x, centroids.iloc[i].y
    for j in neighbours:
        xj, yj = centroids.iloc[j].x, centroids.iloc[j].y
        ax.plot([xi, xj], [yi, yj], 'b-', alpha=0.5, linewidth=1)

ax.scatter(centroids.x, centroids.y, c='red', s=40, zorder=5)
ax.set_title(f'Saudi Arabia adjacency ({len(saudi)} regions, '
              f'avg {w.mean_neighbors:.1f} neighbors)', fontweight='bold')
ax.axis('off')
plt.tight_layout(); plt.show()

def create_region_mapping(gdf, name_col='NAME_1'):
    """Create bidirectional mapping between region names and IDs.

    Parameters
    ----------
    gdf : geopandas.GeoDataFrame
        Administrative boundaries (row order defines IDs)
    name_col : str
        Column containing region names

    Returns
    -------
    dict with keys:
        'id_to_name': {1: 'Region A', 2: 'Region B', ...}
        'name_to_id': {'Region A': 1, 'Region B': 2, ...}
        'names': ['Region A', 'Region B', ...] in order
    """
    names = gdf[name_col].tolist()

    return {
        'id_to_name': {i + 1: name for i, name in enumerate(names)},
        'name_to_id': {name: i + 1 for i, name in enumerate(names)},
        'names': names
    }

# Create mapping
mapping = create_region_mapping(saudi, 'NAME_1')

# Display the mapping
print("Region ID Mapping:")
print("-" * 40)
for id, name in mapping['id_to_name'].items():
    print(f"  ID {id:2d} -> {name}")

# Display neighbor structure with actual names
print("Neighbor Structure:")
print("=" * 60)
for i in range(len(saudi)):
    region_name = mapping['id_to_name'][i + 1]
    neighbor_indices = adj['neighbors'][i]
    neighbor_names = [mapping['id_to_name'][j + 1] for j in neighbor_indices]
    print(f"{region_name}:")
    print(f"    Neighbors: {neighbor_names}")

import pandas as pd
import numpy as np

# Example: Create data with region names
np.random.seed(42)
n_regions = len(saudi)

# Your raw data might have region names
data = pd.DataFrame({
    'region_name': mapping['names'],
    'cases': np.random.poisson(50, n_regions),
    'population': np.random.randint(100000, 5000000, n_regions)
})

# ADD THE SPATIAL ID COLUMN - this links data to the graph
data['region_id'] = data['region_name'].map(mapping['name_to_id'])

# Add offset for Poisson model (log of expected count)
data['log_pop'] = np.log(data['population'])

print("Data with spatial ID:")
print(data)

from pyinla import pyinla

# Model specification for BYM2 (or use 'besag')
model = {
    'response': 'cases',                    # Response variable
    'fixed': [],                             # No fixed effects (intercept only)
    'random': [
        {
            'id': 'region_id',             # Column with spatial IDs (1, 2, 3...)
            'model': 'bym2',               # Spatial model type
            'graph': 'saudi.graph'        # Path to graph file
        }
    ]
}

# Fit the model
# Note: E and offset must be passed as numeric arrays, not column names.
result = pyinla(
    model=model,
    family='poisson',
    data=data,
    E=data['population'].to_numpy(),   # expected count per region (numeric array)
    keep=True
)
# Equivalent log-offset form: offset=data['log_pop'].to_numpy()

# View results
print("Fixed Effects:")
print(result.summary_fixed)

print("\nHyperparameters:")
print(result.summary_hyperpar)

print("\nSpatial Random Effects (first 5):")
print(result.summary_random['region_id'].head())

# Get posterior summaries; for BYM2 take the first n rows (combined b+u).
n = len(saudi)
posterior = result.summary_random['region_id'].iloc[:n]

# Add to GeoDataFrame for plotting
saudi_plot = saudi.copy()
saudi_plot['posterior_mean'] = posterior['mean'].values
saudi_plot['posterior_sd'] = posterior['sd'].values

# For relative risk interpretation (if using Poisson)
saudi_plot['relative_risk'] = np.exp(posterior['mean'].values)

import matplotlib.pyplot as plt

def plot_choropleth(gdf, values, title, cmap='RdYlBu_r'):
    """Plot values on a choropleth map."""
    fig, ax = plt.subplots(figsize=(10, 8))

    gdf_plot = gdf.copy()
    gdf_plot['value'] = values

    gdf_plot.plot(
        column='value',
        ax=ax,
        legend=True,
        legend_kwds={'shrink': 0.7, 'label': title},
        cmap=cmap,
        edgecolor='black',
        linewidth=0.5
    )

    ax.set_title(title, fontsize=14, fontweight='bold')
    ax.axis('off')
    plt.tight_layout()
    return fig, ax

# Example: Plot simulated posterior means
np.random.seed(123)
posterior_mean = np.random.normal(0, 0.4, len(saudi))

fig, ax = plot_choropleth(saudi, posterior_mean, "Spatial Random Effect")
plt.savefig('spatial_effect.png', dpi=150, bbox_inches='tight')
plt.show()

def plot_map_with_labels(gdf, values, name_col='NAME_1', title="Map"):
    """Plot choropleth with region name labels."""
    fig, ax = plt.subplots(figsize=(12, 10))

    gdf_plot = gdf.copy()
    gdf_plot['value'] = values

    gdf_plot.plot(column='value', ax=ax, legend=True, cmap='RdYlBu_r',
                  edgecolor='black', linewidth=0.5)

    # Add labels at region centroids
    for idx, row in gdf_plot.iterrows():
        centroid = row.geometry.centroid
        ax.annotate(
            row[name_col],
            xy=(centroid.x, centroid.y),
            ha='center', va='center',
            fontsize=7, fontweight='bold',
            bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7)
        )

    ax.set_title(title, fontsize=14, fontweight='bold')
    ax.axis('off')
    return fig, ax

# Plot relative risk with labels
np.random.seed(456)
relative_risk = np.exp(np.random.normal(0, 0.3, len(saudi)))

fig, ax = plot_map_with_labels(saudi, relative_risk, 'NAME_1', "Relative Risk by Region")
plt.show()

def plot_mean_and_sd(gdf, mean_values, sd_values, title_prefix=""):
    """Plot posterior mean and SD side by side."""
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))

    gdf_plot = gdf.copy()
    gdf_plot['mean'] = mean_values
    gdf_plot['sd'] = sd_values

    # Left: Posterior mean
    gdf_plot.plot(column='mean', ax=axes[0], legend=True, cmap='RdYlBu_r',
                  edgecolor='black', linewidth=0.3)
    axes[0].set_title(f'{title_prefix}Posterior Mean', fontweight='bold')
    axes[0].axis('off')

    # Right: Posterior SD
    gdf_plot.plot(column='sd', ax=axes[1], legend=True, cmap='YlOrRd',
                  edgecolor='black', linewidth=0.3)
    axes[1].set_title(f'{title_prefix}Posterior SD (Uncertainty)', fontweight='bold')
    axes[1].axis('off')

    plt.tight_layout()
    return fig, axes

# Example
np.random.seed(789)
mean_values = np.random.normal(0, 0.5, len(saudi))
sd_values = np.abs(np.random.normal(0.15, 0.03, len(saudi)))

fig, axes = plot_mean_and_sd(saudi, mean_values, sd_values, "Spatial Effect: ")
plt.show()

def plot_adjacency(gdf, w, title="Adjacency Structure"):
    """Plot map with neighbor connections."""
    fig, ax = plt.subplots(figsize=(10, 8))

    # Plot regions
    gdf.plot(ax=ax, facecolor='#e8f4f8', edgecolor='#333', linewidth=0.5)

    # Get centroids
    centroids = gdf.geometry.centroid

    # Draw edges between neighbors
    for i, neighbors in w.neighbors.items():
        xi, yi = centroids.iloc[i].x, centroids.iloc[i].y
        for j in neighbors:
            xj, yj = centroids.iloc[j].x, centroids.iloc[j].y
            ax.plot([xi, xj], [yi, yj], 'b-', alpha=0.5, linewidth=1)

    # Draw centroids
    ax.scatter(centroids.x, centroids.y, c='red', s=30, zorder=5)

    ax.set_title(f"{title}\n({len(gdf)} regions, avg {w.mean_neighbors:.1f} neighbors)",
                 fontweight='bold')
    ax.axis('off')
    return fig, ax

# Visualize the neighbor structure
fig, ax = plot_adjacency(saudi, adj['weights'], "Saudi Arabia Regions")
plt.show()

import matplotlib.pyplot as plt

lebanon_gov  = get_admin_boundaries('LBN', 1)   # 8 governorates
lebanon_dist = get_admin_boundaries('LBN', 2)   # 30 districts

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 7))
lebanon_gov.plot(ax=ax1, edgecolor='black', linewidth=0.7, facecolor='#e8f4f8')
ax1.set_title(f'Level 1: Governorates (n={len(lebanon_gov)})', fontweight='bold')
ax1.axis('off')

lebanon_dist.plot(ax=ax2, edgecolor='black', linewidth=0.5, facecolor='#fff7ed')
ax2.set_title(f'Level 2: Districts (n={len(lebanon_dist)})', fontweight='bold')
ax2.axis('off')

plt.tight_layout()
plt.show()

# Level 1: Governorates (Muhafazat)
lebanon_gov = gpd.read_file(
    "https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_LBN_1.json"
)
print(f"Level 1: {len(lebanon_gov)} governorates")

# Level 2: Districts (Caza)
lebanon_dist = gpd.read_file(
    "https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_LBN_2.json"
)
print(f"Level 2: {len(lebanon_dist)} districts")

fig, ax = plt.subplots(figsize=(8, 9))
lebanon_gov.plot(ax=ax, edgecolor='black', linewidth=0.7, facecolor='#e8f4f8')

for _, row in lebanon_gov.iterrows():
    c = row.geometry.centroid
    ax.annotate(row['NAME_1'], (c.x, c.y),
                ha='center', va='center', fontsize=8, fontweight='bold',
                bbox=dict(boxstyle='round,pad=0.25', facecolor='white', alpha=0.85))

ax.set_title(f'Lebanon: {len(lebanon_gov)} governorates (ADM1)', fontweight='bold')
ax.axis('off')
plt.tight_layout(); plt.show()

from libpysal.weights import Queen

w         = Queen.from_dataframe(lebanon_dist)
centroids = lebanon_dist.geometry.centroid

fig, ax = plt.subplots(figsize=(11, 9))
lebanon_dist.plot(ax=ax, facecolor='#fff7ed', edgecolor='#333', linewidth=0.4)

for i, nbrs in w.neighbors.items():
    xi, yi = centroids.iloc[i].x, centroids.iloc[i].y
    for j in nbrs:
        xj, yj = centroids.iloc[j].x, centroids.iloc[j].y
        ax.plot([xi, xj], [yi, yj], 'b-', alpha=0.45, linewidth=0.8)

ax.scatter(centroids.x, centroids.y, c='red', s=20, zorder=5)
ax.set_title(f'Lebanon ADM2: {len(lebanon_dist)} districts '
              f'(avg {w.mean_neighbors:.1f} neighbors)', fontweight='bold')
ax.axis('off')
plt.tight_layout(); plt.show()

import numpy as np

rng = np.random.default_rng(0)
sim = rng.normal(0, 0.4, len(lebanon_dist))   # replace with result.summary_random[...]

fig, ax = plt.subplots(figsize=(8, 9))
lebanon_dist.assign(effect=sim).plot(
    column='effect', ax=ax, legend=True, cmap='RdYlBu_r',
    edgecolor='black', linewidth=0.4,
    legend_kwds={'shrink': 0.6, 'label': 'Spatial effect'},
)
ax.set_title('Lebanon ADM2: simulated spatial random effect', fontweight='bold')
ax.axis('off')
plt.tight_layout(); plt.show()

pip install geopandas libpysal matplotlib pandas numpy

import geopandas as gpd
import pandas as pd
import numpy as np
from libpysal.weights import Queen
from pyinla import pyinla

# 1. Download boundaries
gdf = gpd.read_file("https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_SAU_1.json")

# 2. Build and save graph
w = Queen.from_dataframe(gdf)
n = len(gdf)
lines = [str(n)]
for i in range(n):
    neighbors = sorted([j + 1 for j in w.neighbors[i]])
    lines.append(f"{i + 1} {len(neighbors)} {' '.join(map(str, neighbors))}")
with open('graph.graph', 'w') as f:
    f.write('\n'.join(lines))

# 3. Prepare data with spatial ID
data = pd.DataFrame({
    'y': np.random.poisson(50, n),
    'region_id': range(1, n + 1)
})

# 4. Fit model
model = {
    'response': 'y',
    'fixed': [],
    'random': [{'id': 'region_id', 'model': 'bym2', 'graph': 'graph.graph'}]
}
result = pyinla(model=model, family='poisson', data=data)

# 5. Plot results (BYM2 returns 2n rows; first n hold the combined effect b+u)
gdf['effect'] = result.summary_random['region_id']['mean'].values[:n]
gdf.plot(column='effect', legend=True, cmap='RdYlBu_r')