add urban_growth example

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`.

examples/urban_growth/run.py

@@ -0,0 +1,3 @@
+from urban_growth.server import server
+
+server.launch()

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

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,
+)

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
+ )