Chapter 7: Stochastic Cellular Automata for Urban Growth

7.1 From PDEs to Discrete Cells

While PDEs provide elegant analytical results, real urban landscapes are discrete: individual parcels of land are either developed or not, and the state of each parcel depends on its neighbours. Cellular automata (CA) provide a natural framework for modelling such discrete spatial dynamics.

In a stochastic CA for urban growth, each cell$(i,j)$ on a lattice has a state$s_{ij} \in \{0, 1\}$ (undeveloped or developed). At each time step, the transition probability depends on the states of neighbouring cells, land suitability, and a stochastic component.

7.2 Neighbourhood Definitions

Von Neumann Neighbourhood (4-connected)

The four orthogonal neighbours at Manhattan distance 1:

$$\mathcal{N}_{\text{vN}}(i,j) = \{(i\pm 1, j), (i, j\pm 1)\}, \qquad |\mathcal{N}_{\text{vN}}| = 4$$

Moore Neighbourhood (8-connected)

All eight neighbours including diagonals:

$$\mathcal{N}_{\text{M}}(i,j) = \{(i+\delta_i, j+\delta_j) : \delta_i, \delta_j \in \{-1,0,1\}, (\delta_i,\delta_j) \neq (0,0)\}, \qquad |\mathcal{N}_{\text{M}}| = 8$$

Extended Moore with Distance Weighting

For radius $R$, the neighbourhood includes all cells within Chebyshev distance $R$, with weights decaying with distance:

$$w_k = \frac{1}{d_k^\alpha}, \qquad d_k = \max(|\delta_{i,k}|, |\delta_{j,k}|)$$

where $\alpha$ controls the distance decay (typically$\alpha = 1$ or$\alpha = 2$).

7.3 Transition Probability

The neighbourhood influence on cell $(i,j)$ is the weighted average of neighbouring states:

$$\Psi_{ij} = \frac{\sum_{k \in \mathcal{N}} w_k \cdot s_k}{\sum_{k \in \mathcal{N}} w_k}$$

The transition probability for an undeveloped cell to become developed is:

$$P(s_{ij}: 0 \to 1) = f(\Psi_{ij}, \, z_{ij}, \, r_{ij})$$

where:

• $\Psi_{ij}$: neighbourhood influence (fraction of developed neighbours, distance-weighted)
• $z_{ij}$: land suitability (topography, soil, zoning — static or slowly varying)
• $r_{ij}$: stochastic perturbation (captures unmodelled heterogeneity)

A common functional form is the logistic (sigmoid):

$$P(0 \to 1) = \frac{1}{1 + \exp\left(-\beta_0 - \beta_1 \Psi_{ij} - \beta_2 z_{ij} - \beta_3 r_{ij}\right)}$$

Irreversibility Constraint

In most urban CA models, development is irreversible:$P(1 \to 0) = 0$. Once land is urbanised, it stays urbanised. This introduces a fundamental asymmetry that drives the monotonic expansion of the urban boundary.

7.4 Statistical Properties of CA Growth

Stochastic urban CA produce fractal-like growth patterns. Key statistics include:

Urban Patch Compactness

$$C = \frac{4\pi A}{P^2}$$

where $A$ is the urban area and$P$ is the perimeter.$C = 1$ for a circle;$C \to 0$ for fractal sprawl.

Fractal Dimension

$$N(\epsilon) \sim \epsilon^{-D_f} \qquad \Rightarrow \qquad D_f = -\lim_{\epsilon \to 0} \frac{\ln N(\epsilon)}{\ln \epsilon}$$

Real cities typically have $D_f \in [1.5, 1.9]$, indicating significant fractal structure at the urban-rural boundary.

7.5 Python: Stochastic CA Urban Growth

A full implementation of a stochastic CA with configurable Moore neighbourhood, land suitability, and sigmoid transition probability.

Stochastic CA: Urban Growth Simulation

Python

script.py107 lines

import numpy as np

np.random.seed(42)

# ── Stochastic Urban CA ──
N_grid = 60       # grid size
n_steps = 40      # time steps
beta_0 = -2.0     # base log-odds (controls overall growth rate)
beta_1 = 6.0      # neighbourhood influence weight
beta_2 = 1.5      # suitability weight
beta_3 = 0.5      # stochastic weight

# Suitability map (higher near center and along an "axis" = road)
x = np.arange(N_grid)
y = np.arange(N_grid)
xx, yy = np.meshgrid(x, y)
center = N_grid // 2
dist_center = np.sqrt((xx - center)**2 + (yy - center)**2)
suitability = np.exp(-dist_center / (N_grid * 0.4))
# Add road effect (horizontal axis)
suitability += 0.3 * np.exp(-(yy - center)**2 / 8)
suitability /= np.max(suitability)

# Initial seed: small cluster at center
grid = np.zeros((N_grid, N_grid), dtype=int)
grid[center-2:center+3, center-2:center+3] = 1

def moore_influence(grid, i, j, radius=1):
    """Compute distance-weighted Moore neighbourhood influence."""
    total_w = 0.0
    total_ws = 0.0
    for di in range(-radius, radius+1):
        for dj in range(-radius, radius+1):
            if di == 0 and dj == 0:
                continue
            ni, nj = i + di, j + dj
            if 0 <= ni < N_grid and 0 <= nj < N_grid:
                d = max(abs(di), abs(dj))  # Chebyshev distance
                w = 1.0 / d
                total_w += w
                total_ws += w * grid[ni, nj]
    return total_ws / total_w if total_w > 0 else 0.0

print("=" * 65)
print("  Stochastic Cellular Automata: Urban Growth Simulation")
print("=" * 65)
print(f"  Grid: {N_grid}x{N_grid}, Steps: {n_steps}")
print(f"  Coefficients: beta_0={beta_0}, beta_1={beta_1}, beta_2={beta_2}")
print()

stats = []
for step in range(n_steps):
    new_cells = 0
    grid_new = grid.copy()
    for i in range(N_grid):
        for j in range(N_grid):
            if grid[i, j] == 1:
                continue  # already developed
            psi = moore_influence(grid, i, j, radius=2)
            if psi == 0:
                continue  # no developed neighbours in range
            z = suitability[i, j]
            r = np.random.normal(0, 1)
            logit = beta_0 + beta_1 * psi + beta_2 * z + beta_3 * r
            prob = 1.0 / (1.0 + np.exp(-logit))
            if np.random.random() < prob:
                grid_new[i, j] = 1
                new_cells += 1
    grid = grid_new
    total_urban = np.sum(grid)
    frac = total_urban / (N_grid * N_grid)
    stats.append((step+1, total_urban, new_cells, frac))

print(f"  {'Step':>6} {'Urban cells':>13} {'New cells':>11} {'Urban frac':>12}")
print("-" * 46)
for s in stats[::4]:
    print(f"  {s[0]:6d} {s[1]:13d} {s[2]:11d} {s[3]:12.4f}")

# Final statistics
total = np.sum(grid)
# Compute perimeter (cells with at least one non-urban neighbour)
perimeter = 0
for i in range(N_grid):
    for j in range(N_grid):
        if grid[i, j] == 1:
            for di, dj in [(-1,0),(1,0),(0,-1),(0,1)]:
                ni, nj = i+di, j+dj
                if ni < 0 or ni >= N_grid or nj < 0 or nj >= N_grid or grid[ni,nj] == 0:
                    perimeter += 1
                    break

compactness = 4 * np.pi * total / (perimeter**2) if perimeter > 0 else 0

print()
print(f"  Final urban area:  {total} cells")
print(f"  Perimeter:         {perimeter} cells")
print(f"  Compactness (C):   {compactness:.4f}  (1 = circle, 0 = fractal)")
print()

# Display a slice of the final grid
print("  Final grid (central 30x30 region, '#' = urban):")
offset = 15
for i in range(center - offset, center + offset):
    row = ""
    for j in range(center - offset, center + offset):
        row += "#" if grid[i, j] == 1 else "."
    print(f"  {row}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7.6 Fortran: High-Performance CA Core

A Fortran implementation of the CA core with periodic boundary conditions for performance-critical large-grid simulations.

Fortran: High-Performance Stochastic CA

Fortran

program.f9084 lines

program stochastic_ca
  implicit none
  integer, parameter :: N = 40, n_steps = 30
  integer :: grid(N, N), grid_new(N, N)
  real(8) :: suit(N, N), psi, z_val, prob, logit, rnd
  integer :: i, j, di, dj, ni, nj, step, new_cells, total
  integer :: center, d_cheb
  real(8) :: w, total_w, total_ws
  real(8), parameter :: beta0 = -2.0d0, beta1 = 6.0d0, beta2 = 1.5d0

center = N / 2

! Initialize suitability (radial decay from center)
  do i = 1, N
    do j = 1, N
      suit(i,j) = exp(-sqrt(dble((i-center)**2 + (j-center)**2)) / (0.4d0*N))
    end do
  end do

! Initialize grid with central seed
  grid = 0
  do i = center-1, center+1
    do j = center-1, center+1
      grid(i,j) = 1
    end do
  end do

write(*,'(A)')    '================================================='
  write(*,'(A)')    '  Fortran Stochastic CA (Periodic Boundaries)'
  write(*,'(A)')    '================================================='
  write(*,'(A)')    '  Step   Urban   NewCells   Fraction'
  write(*,'(A)')    '-------------------------------------------------'

do step = 1, n_steps
    grid_new = grid
    new_cells = 0

do i = 1, N
      do j = 1, N
        if (grid(i,j) == 1) cycle

! Moore neighbourhood (radius 2) with periodic BC
        total_w  = 0.0d0
        total_ws = 0.0d0
        do di = -2, 2
          do dj = -2, 2
            if (di == 0 .and. dj == 0) cycle
            ni = modulo(i + di - 1, N) + 1
            nj = modulo(j + dj - 1, N) + 1
            d_cheb = max(abs(di), abs(dj))
            w = 1.0d0 / dble(d_cheb)
            total_w  = total_w  + w
            total_ws = total_ws + w * dble(grid(ni, nj))
          end do
        end do

psi = total_ws / total_w
        if (psi < 1.0d-8) cycle

z_val = suit(i,j)
        logit = beta0 + beta1 * psi + beta2 * z_val
        prob = 1.0d0 / (1.0d0 + exp(-logit))

call random_number(rnd)
        if (rnd < prob) then
          grid_new(i,j) = 1
          new_cells = new_cells + 1
        end if
      end do
    end do

grid = grid_new
    total = sum(grid)

if (mod(step, 5) == 0 .or. step == 1) then
      write(*,'(I6, I8, I10, F12.4)') step, total, new_cells, dble(total)/dble(N*N)
    end if
  end do

write(*,'(A)')    '-------------------------------------------------'
  write(*,'(A,I6)') '  Final urban cells: ', sum(grid)
  write(*,'(A)')    '================================================='

end program stochastic_ca

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

7.7 Summary & Key Takeaways

• Stochastic CA model urban growth on discrete lattices with probabilistic transition rules
• Neighbourhood influence: $\Psi_{ij} = \sum w_k s_k / \sum w_k$
• Transition probability: sigmoid function of neighbourhood influence, suitability, and noise
• Moore (8-connected) vs von Neumann (4-connected) neighbourhoods produce different growth morphologies
• CA naturally produce fractal urban boundaries with $D_f \in [1.5, 1.9]$

← Instanton Formalism Next: SLEUTH Model →