diff --git a/examples/urban_growth/README.md b/examples/urban_growth/README.md
new file mode 100644
index 00000000..d1516965
--- /dev/null
+++ b/examples/urban_growth/README.md
@@ -0,0 +1,14 @@
+UrbanGrowth Model
+==============
+
+An implementation of the [UrbanGrowth Model](https://github.com/abmgis/abmgis/tree/master/Chapter06-IntegratingABMandGIS/Models/UrbanGrowth) in Python, using [Mesa](https://github.com/projectmesa/mesa) and [mesa-geo](https://github.com/projectmesa/mesa-geo).
+
+## How to run
+
+To run the model interactively, run `mesa runserver` in this directory. e.g.
+
+```bash
+mesa runserver
+```
+
+Then open your browser to [http://127.0.0.1:8521/](http://127.0.0.1:8521/) and press `Start`.
diff --git a/examples/urban_growth/data/excluded_santafe.asc.gz b/examples/urban_growth/data/excluded_santafe.asc.gz
new file mode 100755
index 00000000..43e7e2cb
Binary files /dev/null and b/examples/urban_growth/data/excluded_santafe.asc.gz differ
diff --git a/examples/urban_growth/data/landuse_santafe.asc.gz b/examples/urban_growth/data/landuse_santafe.asc.gz
new file mode 100755
index 00000000..5de1a72b
Binary files /dev/null and b/examples/urban_growth/data/landuse_santafe.asc.gz differ
diff --git a/examples/urban_growth/data/road1_santafe.asc.gz b/examples/urban_growth/data/road1_santafe.asc.gz
new file mode 100755
index 00000000..6971459f
Binary files /dev/null and b/examples/urban_growth/data/road1_santafe.asc.gz differ
diff --git a/examples/urban_growth/data/slope_santafe.asc.gz b/examples/urban_growth/data/slope_santafe.asc.gz
new file mode 100755
index 00000000..4697b913
Binary files /dev/null and b/examples/urban_growth/data/slope_santafe.asc.gz differ
diff --git a/examples/urban_growth/data/urban_santafe.asc.gz b/examples/urban_growth/data/urban_santafe.asc.gz
new file mode 100755
index 00000000..85f72aed
Binary files /dev/null and b/examples/urban_growth/data/urban_santafe.asc.gz differ
diff --git a/examples/urban_growth/run.py b/examples/urban_growth/run.py
new file mode 100644
index 00000000..7291302e
--- /dev/null
+++ b/examples/urban_growth/run.py
@@ -0,0 +1,3 @@
+from urban_growth.server import server
+
+server.launch()
diff --git a/examples/urban_growth/urban_growth/model.py b/examples/urban_growth/urban_growth/model.py
new file mode 100644
index 00000000..63bbf465
--- /dev/null
+++ b/examples/urban_growth/urban_growth/model.py
@@ -0,0 +1,160 @@
+import random
+
+import mesa
+import numpy as np
+
+from .space import City
+
+
+class UrbanGrowth(mesa.Model):
+    def __init__(
+        self,
+        world_width=531,
+        world_height=394,
+        max_coefficient=100,
+        dispersion_coefficient=20,
+        spread_coefficient=27,
+        breed_coefficient=5,
+        rg_coefficient=10,
+        slope_coefficient=50,
+        critical_slope=25,
+        road_influence=False,
+    ):
+        super().__init__()
+        self.world_width = world_width
+        self.world_height = world_height
+        self.max_coefficient = max_coefficient
+        self.dispersion_coefficient = dispersion_coefficient
+        self.spread_coefficient = spread_coefficient
+        self.breed_coefficient = breed_coefficient
+        self.rg_coefficient = rg_coefficient
+        self.slope_coefficient = slope_coefficient
+        self.critical_slope = critical_slope
+        self.road_influence = road_influence
+
+        self.dispersion_value = (dispersion_coefficient * 0.005) * (
+            world_width**2 + world_height**2
+        ) ** 0.5
+        self.rg_value = (
+            rg_coefficient / max_coefficient * ((world_width + world_height) / 16.0)
+        )
+        self.max_search = 4 * (self.rg_value * (1 + self.rg_value))
+
+        self._load_data()
+        self._check_suitability()
+        self.space.check_road()
+        for cell in self.space.raster_layer:
+            if cell.road:
+                cell.run_value = cell.road_1 / 4 * self.dispersion_coefficient
+            cell.road_found = False
+            cell.road_pixel = None
+
+        self.schedule = mesa.time.RandomActivation(self)
+        self.datacollector = mesa.DataCollector(
+            {
+                "Percentage Urbanized": "pct_urbanized",
+            }
+        )
+        self.datacollector.collect(self)
+
+    @property
+    def pct_urbanized(self) -> float:
+        total_urban = sum([cell.urban for cell in self.space.raster_layer])
+        return total_urban / (self.world_width * self.world_height) * 100.0
+
+    def _load_data(self) -> None:
+        self.space = City(
+            width=self.world_width,
+            height=self.world_height,
+            crs="epsg:3857",
+            total_bounds=[-901575.0, 1442925.0, -885645.0, 1454745.0],
+        )
+        self.space.load_datasets(
+            urban_data="data/urban_santafe.asc.gz",
+            slope_data="data/slope_santafe.asc.gz",
+            road_data="data/road1_santafe.asc.gz",
+            excluded_data="data/excluded_santafe.asc.gz",
+            land_use_data="data/landuse_santafe.asc.gz",
+        )
+
+    def _check_suitability(self) -> None:
+        max_slope = max([cell.slope for cell in self.space.raster_layer])
+
+        prob_to_build = []
+
+        i = 0
+        while i <= self.critical_slope:
+            i += 1
+            val = (self.critical_slope - i) / self.critical_slope
+            prob_to_build.append((val ** (self.slope_coefficient / 200)))
+
+        j = 1
+        while j <= (max_slope - self.critical_slope):
+            prob_to_build.append(0)
+            j += 1
+
+        pct_suitable = 0
+
+        for cell in self.space.raster_layer:
+            if (
+                int(cell.slope) != -9999
+                and np.random.uniform() >= prob_to_build[int(cell.slope)]
+            ) or (cell.excluded == 0.0):
+                cell.suitable = False
+            else:
+                cell.suitable = True
+                pct_suitable += 1
+
+        pct_suitable = pct_suitable / (self.world_width * self.world_height) * 100
+
+    def step(self):
+        self._spontaneous_growth()
+        self._new_spreading_center_growth()
+        self._edge_growth()
+        if self.road_influence:
+            self._road_influenced_growth()
+        self.datacollector.collect(self)
+
+    def _spontaneous_growth(self) -> None:
+        i = 0
+        while i < self.dispersion_value:
+            i += 1
+            w = np.random.randint(self.world_width)
+            h = np.random.randint(self.world_height)
+            random_cell = self.space.raster_layer.cells[w][h]
+            if (not random_cell.urban) and random_cell.suitable:
+                random_cell.urban = True
+                random_cell.new_urbanized = True
+
+    def _new_spreading_center_growth(self) -> None:
+        for cell in self.space.raster_layer:
+            if cell.new_urbanized:
+                x = np.random.randint(self.max_coefficient)
+                if x < self.breed_coefficient:
+                    neighbors = self.space.raster_layer.get_neighboring_cells(
+                        cell.pos, moore=True
+                    )
+                    for random_neighbor in random.choices(neighbors, k=2):
+                        if (not random_neighbor.urban) and random_neighbor.suitable:
+                            random_neighbor.urban = True
+                            random_neighbor.new_urbanized = True
+                cell.new_urbanized = False
+
+    def _edge_growth(self) -> None:
+        for cell in self.space.raster_layer:
+            if cell.urban:
+                x = np.random.randint(self.max_coefficient)
+                if x < self.spread_coefficient:
+                    neighbors = self.space.raster_layer.get_neighboring_cells(
+                        cell.pos, moore=True
+                    )
+                    urban_neighbors = [c for c in neighbors if c.urban]
+                    non_urban_neighbors = [c for c in neighbors if not c.urban]
+                    if len(urban_neighbors) > 1 and len(non_urban_neighbors) > 0:
+                        random_non_urban_neighbor = random.choice(non_urban_neighbors)
+                        if random_non_urban_neighbor.suitable:
+                            random_non_urban_neighbor.urban = True
+                            random_non_urban_neighbor.new_urbanized = True
+
+    def _road_influenced_growth(self) -> None:
+        pass
diff --git a/examples/urban_growth/urban_growth/server.py b/examples/urban_growth/urban_growth/server.py
new file mode 100644
index 00000000..d0e447fa
--- /dev/null
+++ b/examples/urban_growth/urban_growth/server.py
@@ -0,0 +1,63 @@
+from typing import Tuple
+
+import mesa
+
+from mesa_geo.visualization.ModularVisualization import ModularServer
+from mesa_geo.visualization.modules import MapModule
+from .model import UrbanGrowth
+from .space import UrbanCell
+
+
+def cell_portrayal(cell: UrbanCell) -> Tuple[float, float, float, float]:
+    if cell.urban:
+        if cell.new_urbanized:
+            return 255, 0, 0, 1
+        else:
+            return 0, 0, 255, 1
+    else:
+        return 0, 0, 0, 0
+
+
+class UrbanizedText(mesa.visualization.TextElement):
+    def render(self, model):
+        return f"Percentage Urbanized: {model.pct_urbanized:.2f}%"
+
+
+model_params = {
+    "max_coefficient": mesa.visualization.NumberInput("max_coefficient", 100),
+    "dispersion_coefficient": mesa.visualization.Slider(
+        "dispersion_coefficient", 20, 0, 100, 1
+    ),
+    "spread_coefficient": mesa.visualization.Slider(
+        "spread_coefficient", 27, 0, 100, 1
+    ),
+    "breed_coefficient": mesa.visualization.Slider("breed_coefficient", 5, 0, 100, 1),
+    "rg_coefficient": mesa.visualization.Slider("rg_coefficient", 10, 0, 100, 1),
+    "slope_coefficient": mesa.visualization.Slider("slope_coefficient", 50, 0, 100, 1),
+    "critical_slope": mesa.visualization.Slider("critical_slope", 25, 0, 100, 1),
+    "road_influence": mesa.visualization.Choice(
+        "road_influence", False, choices=[True, False]
+    ),
+}
+
+
+map_module = MapModule(
+    portrayal_method=cell_portrayal,
+    view=[12.904598815296707, -8.027435210420451],
+    zoom=12.1,
+    map_height=394,
+    map_width=531,
+)
+urbanized_text = UrbanizedText()
+urbanized_chart = mesa.visualization.ChartModule(
+    [
+        {"Label": "Percentage Urbanized", "Color": "Black"},
+    ]
+)
+
+server = ModularServer(
+    UrbanGrowth,
+    [map_module, urbanized_text, urbanized_chart],
+    "UrbanGrowth Model",
+    model_params,
+)
diff --git a/examples/urban_growth/urban_growth/space.py b/examples/urban_growth/urban_growth/space.py
new file mode 100644
index 00000000..b1bd356d
--- /dev/null
+++ b/examples/urban_growth/urban_growth/space.py
@@ -0,0 +1,89 @@
+from __future__ import annotations
+import gzip
+
+import mesa
+import rasterio as rio
+
+from mesa_geo import Cell, RasterLayer
+from mesa_geo.geospace import GeoSpace
+
+
+class UrbanCell(Cell):
+    urban: bool | None
+    slope: int | None
+    road_1: int | None
+    excluded: int | None
+    land_use: int | None
+
+    suitable: bool | None
+    road: bool | None
+    run_value: int | None
+    road_found: bool | None
+    road_pixel: UrbanCell | None
+
+    new_urbanized: bool | None
+
+    def __init__(
+        self,
+        pos: mesa.space.Coordinate | None = None,
+        indices: mesa.space.Coordinate | None = None,
+    ):
+        super().__init__(pos, indices)
+        self.urban = None
+        self.slope = None
+        self.road_1 = None
+        self.excluded = None
+        self.land_use = None
+
+        self.suitable = None
+        self.road = None
+        self.run_value = None
+        self.road_found = None
+        self.road_pixel = None
+        self.new_urbanized = None
+
+    def step(self):
+        pass
+
+
+class City(GeoSpace):
+    def __init__(self, width, height, crs, total_bounds):
+        super().__init__(crs=crs)
+        self.add_layer(
+            RasterLayer(width, height, crs, total_bounds, cell_cls=UrbanCell)
+        )
+
+    def load_datasets(
+        self, urban_data, slope_data, road_data, excluded_data, land_use_data
+    ):
+        data = {
+            "urban": urban_data,
+            "slope": slope_data,
+            "road_1": road_data,
+            "excluded": excluded_data,
+            "land_use": land_use_data,
+        }
+        for attribute_name, data_file in data.items():
+            with gzip.open(data_file, "rb") as f:
+                with rio.open(f, "r") as dataset:
+                    values = dataset.read()
+            self.raster_layer.apply_raster(values, name=attribute_name)
+
+        for cell in self.raster_layer:
+            cell.urban = True if cell.urban == 2 else False
+
+    def check_road(self):
+        for cell in self.raster_layer:
+            cell.road = True if cell.road_1 > 0 else False
+
+    @property
+    def raster_layer(self):
+        return self.layers[0]
+
+    def is_at_boundary(self, row_idx, col_idx):
+        return (
+            row_idx == 0
+            or row_idx == self.raster_layer.height
+            or col_idx == 0
+            or col_idx == self.raster_layer.width
+        )