diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3de96c1..aadc511 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,15 +1,107 @@
 # Changelog
 
+All notable changes to this project will be documented in this file.
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+---
+
+## [0.2.1] - 2026-03
+
+### Changed
+- `parse_idx` now returns `GridPos(row, col)` namedtuple instead of a plain tuple,
+  eliminating the `(col, row)` vs `(row, col)` ambiguity. Tuple unpacking remains
+  fully compatible — no breaking change for existing callers.
+- `regular_grid()` renamed to `vector_grid()`. Old name kept as a `DeprecationWarning`
+  alias and will be removed in v0.3.0.
+
+### Added
+- `GridPos` namedtuple exported from `dissmodel.geo.vector.regular_grid`.
+
+### Fixed
+- `RasterModel.setup()` and `RasterCellularAutomaton.setup()` no longer accept
+  `**kwargs`, resolving a salabim incompatibility that raised
+  `TypeError: parameter 'kwargs' not allowed`.
+
+---
+
 ## [0.2.0] - 2026-03
 
 ### Added
-- Raster backend architecture
-- Examples for raster simulations
-- CLI examples
 
-### Improved
-- Package organization
-- Documentation
+#### Raster substrate
+- `RasterBackend` — named NumPy array store with vectorized spatial operations:
+  `shift2d`, `focal_sum`, `focal_sum_mask`, `neighbor_contact`, `snapshot`.
+- `RasterModel` — base class for raster push models, providing `self.backend`,
+  `self.shape`, `self.shift`, and `self.dirs` (Moore neighbourhood).
+- `RasterCellularAutomaton` — vectorized CA base class; `rule(arrays) → dict`
+  replaces per-cell iteration.
+- `raster_grid()` / `make_raster_grid()` — `RasterBackend` factory.
+- `DIRS_MOORE` and `DIRS_VON_NEUMANN` — neighbourhood direction constants.
+- `RasterMap` — visualization component supporting categorical (`color_map`) and
+  continuous (`cmap`) modes; renders to Streamlit, Jupyter, interactive window,
+  or headless PNG frames.
+
+#### Vector substrate
+- `SpatialModel` — GeoDataFrame-based push model with `create_neighborhood()`,
+  `neighs_id()`, `neighs()`, and `neighbor_values()`.
+- `vector_grid()` — replaces `regular_grid()` as the canonical grid factory name.
+
+#### Examples and benchmarks
+- `GameOfLifeRaster` — Conway's Game of Life on raster substrate.
+- `FireModelRaster` — forest fire spread on raster substrate.
+- `benchmark_game_of_life.py` — vector vs raster benchmark with exact cell-by-cell
+  validation; supports `--steps`, `--sizes`, `--no-validation` flags.
+- `benchmark_raster_vs_vector.py` — flood model benchmark at realistic workload.
+- CLI examples updated.
+
+#### Tests
+- Full test suite: `tests/vector/`, `tests/raster/`, `tests/integration/`.
+- `tests/integration/test_game_of_life.py` — exact cell-by-cell validation across
+  5×5, 10×10, 20×20 grids and 5 seeds.
+- `tests/integration/test_flood_model.py` — cross-substrate equivalence with 95%
+  match threshold.
+
+### Changed
+- Package reorganized: `dissmodel.geo.vector.*` and `dissmodel.geo.raster.*`
+  replace the flat `dissmodel.geo.*` layout.
+- Documentation migrated to MkDocs Material with separate API reference pages for
+  vector and raster substrates.
 
 ### Fixed
-- Import paths in examples
\ No newline at end of file
+- Import paths corrected across all examples after package reorganization.
+- `celullar_automaton.py` filename typo fixed → `cellular_automaton.py`.
+
+### Performance
+Benchmarks on Conway's Game of Life (10 steps, Python 3.12, NumPy):
+
+| Grid | Cells | Raster (ms/step) | Vector (ms/step) | Speedup |
+|-----:|------:|-----------------:|-----------------:|--------:|
+| 10×10 | 100 | 0.15 | 30.11 | 206× |
+| 50×50 | 2,500 | 0.20 | 647.22 | 3,164× |
+| 100×100 | 10,000 | 0.60 | 2,715.16 | 4,491× |
+| 1,000×1,000 | 1,000,000 | 25.85 | — | — |
+
+---
+
+## [0.1.5] - 2026-02
+
+### Added
+- JOSS submission at this version.
+
+### Fixed
+- Minor documentation and packaging fixes.
+
+---
+
+## [0.1.0] - 2025
+
+### Added
+- Initial release.
+- `CellularAutomaton` and `SpatialModel` on GeoDataFrame substrate.
+- `regular_grid()`, `fill()`, `FillStrategy`.
+- System Dynamics models: `SIR`, `PredatorPrey`, `PopulationGrowth`, `Lorenz`, `Coffee`.
+- Cellular Automata models: `GameOfLife`, `FireModel`, `FireModelProb`, `Snow`, `Growth`,
+  `Propagation`, `Anneal`.
+- Salabim integration via `dissmodel.core.Environment` and `Model`.
+- Streamlit examples.
diff --git a/dissmodel/visualization/__init__.py b/dissmodel/visualization/__init__.py
index 60e190d..2aaab5b 100644
--- a/dissmodel/visualization/__init__.py
+++ b/dissmodel/visualization/__init__.py
@@ -1,3 +1,4 @@
 from dissmodel.visualization.map import Map
+from dissmodel.visualization.raster_map import RasterMap
 from dissmodel.visualization.chart import Chart, track_plot
 from dissmodel.visualization.widgets import display_inputs
\ No newline at end of file
diff --git a/docs/api/core.md b/docs/api/core.md
index 728f816..9dc57e5 100644
--- a/docs/api/core.md
+++ b/docs/api/core.md
@@ -1,4 +1,145 @@
 # Core
 
+The `dissmodel.core` module provides the simulation clock and execution lifecycle,
+built on top of [Salabim](https://www.salabim.org/)'s discrete event engine.
+
+All models and visualization components must be instantiated **after** the
+`Environment` — they register themselves automatically on creation.
+
+```
+Environment  →  Model  →  Visualization  →  env.run()
+     ↑             ↑            ↑                ↑
+  first         second        third           fourth
+```
+
+## Usage
+
+```python
+from dissmodel.core import Environment, Model
+
+env = Environment(start_time=1, end_time=10)
+
+class MyModel(Model):
+    def setup(self):
+        pass
+
+    def execute(self):
+        print(f"step {self.env.now()}")
+
+MyModel()
+env.run()
+```
+
+## Object-Oriented Modeling
+
+Object-oriented modeling is a core feature of DisSModel, inherited directly from
+Python's class system. Just as TerraME defines agents as objects with encapsulated
+attributes and behaviours, DisSModel uses class inheritance to build structured,
+reusable, and modular models.
+
+Every model is a subclass of `Model`, which guarantees automatic registration with
+the active `Environment`. This means the simulation clock, the execution lifecycle,
+and any visualization components are wired together without any boilerplate.
+
+```python
+from dissmodel.core import Model, Environment
+
+class SIR(Model):
+
+    def setup(self, susceptible=9998, infected=2, recovered=0,
+              duration=2, contacts=6, probability=0.25):
+        self.susceptible  = susceptible
+        self.infected     = infected
+        self.recovered    = recovered
+        self.duration     = duration
+        self.contacts     = contacts
+        self.probability  = probability
+
+    def execute(self):
+        total       = self.susceptible + self.infected + self.recovered
+        alpha       = self.contacts * self.probability
+        new_inf     = self.infected * alpha * (self.susceptible / total)
+        new_rec     = self.infected / self.duration
+        self.susceptible -= new_inf
+        self.infected    += new_inf - new_rec
+        self.recovered   += new_rec
+```
+
+Instantiation is clean and parametric:
+
+```python
+env = Environment(end_time=30)
+SIR(susceptible=9998, infected=2, recovered=0,
+    duration=2, contacts=6, probability=0.25)
+env.run()
+```
+
+!!! tip "Why subclass Model?"
+    - **Automatic clock integration** — `self.env.now()` is always available inside `execute()`.
+    - **Encapsulation** — each model owns its state; multiple instances can run in the same environment independently.
+    - **Extensibility** — override `setup()` to add parameters, `execute()` to define the transition rule. Nothing else is required.
+    - **Composability** — models can read each other's state, enabling coupled CA + SysDyn simulations within a single `env.run()`.
+
+
+
+Each model can define its own `start_time` and `end_time`, independent of the
+environment interval. This allows different parts of a simulation to be active
+at different periods within the same run.
+
+```python
+from dissmodel.core import Model, Environment
+
+class ModelA(Model):
+    def execute(self):
+        print(f"[A] t={self.env.now()}")
+
+class ModelB(Model):
+    def execute(self):
+        print(f"[B] t={self.env.now()}")
+
+class ModelC(Model):
+    def execute(self):
+        print(f"[C] t={self.env.now()}")
+
+env = Environment(start_time=2010, end_time=2016)
+
+ModelA(start_time=2012)        # active from 2012 to end
+ModelB(end_time=2013)          # active from start to 2013
+ModelC()                       # active throughout
+
+env.run()
+```
+
+Expected output:
+
+```
+Running from 2010 to 2016 (duration: 6)
+[B] t=2010.0
+[C] t=2010.0
+[B] t=2011.0
+[C] t=2011.0
+[A] t=2012.0
+[B] t=2012.0
+[C] t=2012.0
+[A] t=2013.0
+[C] t=2013.0
+[A] t=2014.0
+[C] t=2014.0
+[A] t=2015.0
+[C] t=2015.0
+[A] t=2016.0
+[C] t=2016.0
+```
+
+!!! note
+    Models with no `start_time` / `end_time` inherit the environment's interval.
+    Models are synchronised — all active models execute at each time step before
+    the clock advances.
+
+---
+
+## API Reference
+
 ::: dissmodel.core.Environment
-::: dissmodel.core.Model
+
+::: dissmodel.core.Model
\ No newline at end of file
diff --git a/docs/api/geo/index.md b/docs/api/geo/index.md
new file mode 100644
index 0000000..f27a0dd
--- /dev/null
+++ b/docs/api/geo/index.md
@@ -0,0 +1,206 @@
+# Geo
+
+The `dissmodel.geo` module provides the spatial infrastructure for building
+simulation models. It handles grid generation, neighbourhood computation, and
+attribute initialization — without imposing any domain logic.
+
+```python
+from dissmodel.geo import vector_grid, fill, FillStrategy
+from dissmodel.geo.vector.neighborhood import attach_neighbors
+from dissmodel.geo.raster.backend import RasterBackend
+from dissmodel.geo.raster.regular_grid import raster_grid
+```
+
+---
+
+## Dual-substrate design
+
+The module provides two independent spatial substrates. Both share the same
+salabim `Environment` and clock — a vector model and a raster model can run
+side by side in the same `env.run()`.
+
+| | Vector | Raster |
+|---|---|---|
+| **Module** | `dissmodel.geo.vector` | `dissmodel.geo.raster` |
+| **Data structure** | `GeoDataFrame` (GeoPandas) | `RasterBackend` (NumPy 2D arrays) |
+| **Grid factory** | `vector_grid()` | `raster_grid()` |
+| **Neighbourhood** | Queen / Rook (libpysal) | Moore / Von Neumann (`shift2d`) |
+| **Rule pattern** | `rule(idx)` per cell | `rule(arrays) → dict` vectorized |
+| **GIS integration** | CRS, projections, spatial joins | rasterio I/O via `load_geotiff` |
+| **Best for** | Irregular grids, real-world data | Large grids, performance-critical models |
+
+---
+
+## Vector substrate
+
+The vector substrate uses a `GeoDataFrame` as the spatial grid. Any model can
+operate directly on real geographic data — shapefiles, GeoJSON, real CRS — with
+no conversion step.
+
+```python
+import geopandas as gpd
+from dissmodel.core import Model, Environment
+from dissmodel.visualization.map import Map
+
+gdf = gpd.read_file("area.shp")
+gdf.set_index("object_id", inplace=True)
+
+env = Environment(start_time=1, end_time=20)
+
+class ElevationModel(Model):
+    def setup(self, gdf, rate=0.01):
+        self.gdf  = gdf
+        self.rate = rate
+
+    def execute(self):
+        self.gdf["alt"] += self.rate
+
+ElevationModel(gdf=gdf, rate=0.01)
+Map(gdf=gdf, plot_params={"column": "alt", "cmap": "Blues", "legend": True})
+env.run()
+```
+
+For abstract (non-georeferenced) grids, use `vector_grid()`:
+
+```python
+from dissmodel.geo import vector_grid
+
+# from dimension + resolution
+gdf = vector_grid(dimension=(10, 10), resolution=1)
+
+# from bounding box + resolution
+gdf = vector_grid(bounds=(0, 0, 1000, 1000), resolution=100)
+
+# from an existing GeoDataFrame
+gdf = vector_grid(gdf=base_gdf, resolution=50)
+```
+
+---
+
+## Raster substrate
+
+The raster substrate stores named NumPy arrays in a `RasterBackend`. All
+operations (`shift2d`, `focal_sum`, `neighbor_contact`) are fully vectorized —
+no Python loops over cells.
+
+```python
+from dissmodel.geo.raster.regular_grid import raster_grid
+from dissmodel.geo.raster.backend import RasterBackend
+import numpy as np
+
+backend = raster_grid(rows=100, cols=100, attrs={"state": 0, "alt": 0.0})
+
+# read / write arrays
+state = backend.get("state").copy()    # snapshot — equivalent to .past in TerraME
+backend.arrays["state"] = new_state
+
+# vectorized neighbourhood operations
+shifted    = RasterBackend.shift2d(state, -1, 0)   # northern neighbour of each cell
+n_active   = backend.focal_sum_mask(state == 1)    # count active Moore neighbours
+has_active = backend.neighbor_contact(state == 1)  # bool mask: any active neighbour?
+```
+
+---
+
+## Filling grid attributes
+
+The `fill()` function initialises GeoDataFrame columns from spatial data sources,
+avoiding manual cell-by-cell loops.
+
+```python
+from dissmodel.geo import fill, FillStrategy
+```
+
+### Zonal statistics from a raster
+
+```python
+import rasterio
+
+with rasterio.open("altitude.tif") as src:
+    raster  = src.read(1)
+    affine  = src.transform
+
+fill(
+    FillStrategy.ZONAL_STATS,
+    vectors=gdf, raster_data=raster, affine=affine,
+    stats=["mean", "min", "max"], prefix="alt_",
+)
+# → adds columns alt_mean, alt_min, alt_max to gdf
+```
+
+### Minimum distance to features
+
+```python
+rivers = gpd.read_file("rivers.shp")
+
+fill(FillStrategy.MIN_DISTANCE, from_gdf=gdf, to_gdf=rivers, attr_name="dist_river")
+```
+
+### Random sampling
+
+```python
+fill(
+    FillStrategy.RANDOM_SAMPLE,
+    gdf=gdf, attr="land_use",
+    data={0: 0.7, 1: 0.3},   # 70% class 0, 30% class 1
+    seed=42,
+)
+```
+
+### Fixed pattern (useful for tests)
+
+```python
+pattern = [[1, 0, 0],
+           [0, 1, 0],
+           [0, 0, 1]]
+
+fill(FillStrategy.PATTERN, gdf=gdf, attr="zone", pattern=pattern)
+```
+
+Custom strategies can be registered:
+
+```python
+from dissmodel.geo.vector.fill import register_strategy
+
+@register_strategy("my_strategy")
+def fill_my_strategy(gdf, attr, **kwargs):
+    ...
+```
+
+---
+
+## Neighbourhood
+
+Spatial neighbourhoods are built via `attach_neighbors()` or directly through
+`create_neighborhood()` on any `CellularAutomaton` or `SpatialModel`.
+
+```python
+from libpysal.weights import Queen, Rook, KNN
+from dissmodel.geo.vector.neighborhood import attach_neighbors
+
+# topological (Queen — edge or vertex contact)
+gdf = attach_neighbors(gdf, strategy=Queen)
+
+# topological (Rook — edge contact only)
+gdf = attach_neighbors(gdf, strategy=Rook)
+
+# distance-based (k nearest neighbours)
+gdf = attach_neighbors(gdf, strategy=KNN, k=4)
+
+# precomputed (from dict or JSON file — faster for large grids)
+gdf = attach_neighbors(gdf, neighbors_dict="neighborhood.json")
+```
+
+| Strategy | Use case |
+|----------|----------|
+| `Queen` | Standard CA — cells share an edge or vertex |
+| `Rook` | Von Neumann-style — edge contact only |
+| `KNN` | Point data, non-contiguous polygons |
+| `neighbors_dict` | Precomputed — skip recomputation on repeated runs |
+
+---
+
+## See also
+
+- [Vector API Reference](vector.md)
+- [Raster API Reference](raster.md)
diff --git a/docs/api/visualization.md b/docs/api/visualization.md
index c0e1798..fe1a7da 100644
--- a/docs/api/visualization.md
+++ b/docs/api/visualization.md
@@ -1,7 +1,182 @@
 # Visualization
 
+The `dissmodel.visualization` module provides graphical and interactive representations
+of running simulations. All visualization components inherit from `Model` and are
+therefore integrated into the simulation clock — they update automatically at each step.
+
+Three main components are available:
+
+| Component | Substrate | Description |
+|-----------|-----------|-------------|
+| `Chart` | Any | Time-series plots from tracked model variables |
+| `Map` | Vector (GeoDataFrame) | Dynamic spatial maps updated each step |
+| `RasterMap` | Raster (NumPy) | Raster array rendering — categorical or continuous |
+
+All three support **three output targets**: local `matplotlib` window,
+Jupyter inline display, and Streamlit `st.empty()` placeholder.
+
+---
+
+## `@track_plot`
+
+The `track_plot` decorator marks model attributes to be collected and plotted
+by `Chart`. Each call defines the variable label, colour, and plot type.
+
+```python
+from dissmodel.core import Model
+from dissmodel.visualization import track_plot
+
+@track_plot("Susceptible", "green")
+@track_plot("Infected",    "red")
+@track_plot("Recovered",   "blue")
+class SIR(Model):
+
+    def setup(self, susceptible=9998, infected=2, recovered=0,
+              duration=2, contacts=6, probability=0.25):
+        self.susceptible = susceptible
+        self.infected    = infected
+        self.recovered   = recovered
+        self.duration    = duration
+        self.contacts    = contacts
+        self.probability = probability
+
+    def execute(self):
+        total   = self.susceptible + self.infected + self.recovered
+        alpha   = self.contacts * self.probability
+        new_inf = self.infected * alpha * (self.susceptible / total)
+        new_rec = self.infected / self.duration
+        self.susceptible -= new_inf
+        self.infected    += new_inf - new_rec
+        self.recovered   += new_rec
+```
+
+---
+
+## `Chart`
+
+Displays time-series data from variables annotated with `@track_plot`.
+
+```python
+from dissmodel.core import Environment
+from dissmodel.models.sysdyn import SIR
+from dissmodel.visualization import Chart
+
+env = Environment(end_time=30)
+SIR()
+Chart(show_legend=True)
+env.run()
+```
+
+**Streamlit:**
+
+```python
+Chart(plot_area=st.empty())
+```
+
+---
+
+## `Map`
+
+Renders spatial data from a GeoDataFrame, updated at every simulation step.
+
+```python
+from dissmodel.visualization.map import Map
+from matplotlib.colors import ListedColormap
+
+Map(
+    gdf=gdf,
+    plot_params={
+        "column": "state",
+        "cmap": ListedColormap(["white", "black"]),
+        "ec": "gray",
+    },
+)
+```
+
+---
+
+## `RasterMap`
+
+Renders a named NumPy array from a `RasterBackend`. Supports categorical
+(value → colour mapping) and continuous (colormap + colorbar) modes.
+
+**Categorical:**
+
+```python
+from dissmodel.visualization.raster_map import RasterMap
+
+RasterMap(
+    backend   = b,
+    band      = "uso",
+    title     = "Land Use",
+    color_map = {1: "#006400", 3: "#00008b", 5: "#d2b48c"},
+    labels    = {1: "Mangrove", 3: "Sea", 5: "Bare soil"},
+)
+```
+
+**Continuous:**
+
+```python
+RasterMap(
+    backend        = b,
+    band           = "alt",
+    title          = "Altimetry",
+    cmap           = "terrain",
+    colorbar_label = "Altitude (m)",
+    mask_band      = "uso",
+    mask_value     = 3,      # mask SEA cells
+)
+```
+
+**Headless** (default when no display is available):
+frames are saved to `raster_map_frames/<band>_step_NNN.png`.
+
+---
+
+## `display_inputs`
+
+Generates Streamlit input widgets automatically from a model's type annotations.
+Integer and float attributes become sliders; booleans become checkboxes.
+
+```python
+from dissmodel.visualization import display_inputs
+
+sir = SIR()
+display_inputs(sir, st.sidebar)
+```
+
+---
+
+## Full Streamlit example
+
+```python
+import streamlit as st
+from dissmodel.core import Environment
+from dissmodel.models.sysdyn import SIR
+from dissmodel.visualization import Chart, display_inputs
+
+st.set_page_config(page_title="SIR Model", layout="centered")
+st.title("SIR Model — DisSModel")
+
+st.sidebar.title("Parameters")
+steps   = st.sidebar.slider("Steps", min_value=1, max_value=50, value=10)
+run_btn = st.button("Run")
+
+env = Environment(end_time=steps, start_time=0)
+sir = SIR()
+display_inputs(sir, st.sidebar)
+Chart(plot_area=st.empty())
+
+if run_btn:
+    env.run()
+```
+
+---
+
+## API Reference
+
 ::: dissmodel.visualization.Chart
+
 ::: dissmodel.visualization.map.Map
-::: dissmodel.visualization.raster_map.RasterMap
-::: dissmodel.visualization.widgets.display_inputs
-::: dissmodel.visualization.track_plot
+
+::: dissmodel.visualization.raster_map.RasterMap
\ No newline at end of file
diff --git a/docs/getting_started.md b/docs/getting_started.md
index 66b61ca..6bc173d 100644
--- a/docs/getting_started.md
+++ b/docs/getting_started.md
@@ -1,5 +1,25 @@
 # Getting Started
 
+## Why Python?
+
+Python has become the standard language for data science and geographic analysis,
+with a rich ecosystem of libraries that complement each other naturally:
+
+- **Greater flexibility** — NumPy, SciPy, Pandas, GeoPandas, rasterio, and
+  machine learning frameworks can all be integrated into the same workflow,
+  allowing models of arbitrary complexity.
+- **Active development** — the Python geospatial community is extremely active,
+  with frequent updates, extensive documentation, and strong support for
+  reproducible research practices.
+- **Accessibility** — Python is widely considered one of the easiest languages
+  to learn, making spatial modeling more approachable for environmental scientists
+  and territorial planners who are not primarily software developers.
+
+DisSModel builds on this ecosystem rather than replacing it. A DisSModel simulation
+is just Python — the full power of the scientific stack is available at every step.
+
+---
+
 ## Installation
 
 ```bash
@@ -14,17 +34,21 @@ cd dissmodel
 pip install -e .
 ```
 
+---
+
 ## Instantiation order
 
 The `Environment` must always be created **before** any model.
 Models connect to the active environment automatically when instantiated.
 
 ```
-Environment  →  Model  →  Visualization
-     ↑             ↑            ↑
-  first         second        third
+Environment  →  Model  →  Visualization  →  env.run()
+     ↑             ↑            ↑                ↑
+  first         second        third           fourth
 ```
 
+---
+
 ## Minimal example
 
 ```python
@@ -33,7 +57,89 @@ from dissmodel.models.sysdyn import SIR
 from dissmodel.visualization import Chart
 
 env = Environment()
-SIR()
+SIR(susceptible=9998, infected=2, recovered=0,
+    duration=2, contacts=6, probability=0.25)
 Chart(show_legend=True)
 env.run(30)
 ```
+
+---
+
+## Execution modes
+
+DisSModel supports three execution strategies. All three follow the same
+`Environment → Model → Visualization → env.run()` pattern — only the
+entry point and how results are displayed differ.
+
+### Command Line (CLI)
+
+Best for automation, batch experiments, and integration with processing pipelines.
+
+```bash
+python examples/cli/sir_model.py
+```
+
+CLI examples are located in `examples/cli/`.
+
+### Jupyter Notebook
+
+Best for incremental exploration, teaching, and visual analysis. DisSModel
+detects the Jupyter environment automatically and renders visualizations inline.
+
+```python
+from dissmodel.core import Environment
+from dissmodel.geo import regular_grid
+from dissmodel.models.ca import GameOfLife
+from dissmodel.visualization.map import Map
+from matplotlib.colors import ListedColormap
+
+gdf = regular_grid(dimension=(30, 30), resolution=1)
+
+env = Environment(end_time=20)
+model = GameOfLife(gdf=gdf)
+model.initialize()
+
+Map(
+    gdf=gdf,
+    plot_params={
+        "column": "state",
+        "cmap": ListedColormap(["white", "black"]),
+        "ec": "gray",
+    },
+)
+
+env.run()
+```
+
+Notebook examples are located in `examples/notebooks/`.
+
+### Streamlit Web Application
+
+Best for interactive demos and non-technical users. Parameters are configured
+via sliders and input fields directly in the browser.
+
+```bash
+streamlit run examples/streamlit/sir_model.py
+
+# or run all models in a single interface:
+streamlit run examples/streamlit/run_all.py
+```
+
+Streamlit examples are located in `examples/streamlit/`.
+
+---
+
+## Choosing a substrate
+
+DisSModel provides two spatial substrates for cellular automata. Choose based
+on your performance and flexibility requirements:
+
+| | Vector | Raster |
+|---|---|---|
+| **Data structure** | GeoDataFrame | NumPy 2D array |
+| **Neighbourhood** | Queen / Rook (libpysal) | Moore / Von Neumann (shift2d) |
+| **Rule pattern** | `rule(idx)` per cell | `rule(arrays) → dict` vectorized |
+| **Performance** | ~2,700 ms/step @ 10k cells | ~0.6 ms/step @ 10k cells |
+| **Best for** | Irregular grids, GIS integration | Large grids, performance-critical models |
+
+See the [API Reference](api/geo/vector.md) for full details on each substrate.
\ No newline at end of file
diff --git a/docs/index.md b/docs/index.md
index 9008ebf..826b1b0 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -7,6 +7,10 @@ Developed by the [LambdaGeo](https://github.com/LambdaGeo) group at the Federal
 it provides a unified environment for building **Cellular Automata (CA)** and **System Dynamics (SysDyn)** models
 on top of the Python geospatial ecosystem.
 
+Inspired by the [TerraME](http://www.terrame.org/) framework, DisSModel brings the same modeling
+expressiveness to Python — replacing the TerraLib/Lua stack with GeoPandas, NumPy, and Salabim,
+while remaining fully interoperable with the broader Python data science ecosystem.
+
 ```bash
 pip install dissmodel
 ```
@@ -15,15 +19,9 @@ pip install dissmodel
 
 ## Why DisSModel?
 
-DisSModel was designed as a modern, Pythonic alternative to the [TerraME](http://www.terrame.org/) framework,
-replacing the TerraLib/Lua stack with the standardized GeoPandas/Python stack.
-This transition gives researchers direct access to the full Python data science ecosystem
-while maintaining the modeling expressiveness required for Land Use and Cover Change (LUCC) applications.
-
-**Core objectives:**
-
 - **Multi-paradigm support** — Cellular Automata, System Dynamics, and Agent-Based Models in a unified environment
-- **Geospatial ecosystem integration** — GeoPandas, libpysal, and rasterstats for advanced spatial operations
+- **Dual-substrate architecture** — vector (GeoDataFrame) for spatial expressiveness, raster (NumPy) for high-performance vectorized computation
+- **Geospatial ecosystem integration** — GeoPandas, libpysal, rasterio, and rasterstats for advanced spatial operations
 - **Open science and reproducibility** — open-source, installable via PyPI, examples included
 - **Standardized implementation** — pure Python lowers the barrier for interdisciplinary collaboration
 
@@ -36,9 +34,30 @@ DisSModel is organized into four modules:
 | Module | Description |
 |:---|:---|
 | `dissmodel.core` | Simulation clock and execution lifecycle powered by [Salabim](https://www.salabim.org/) |
-| `dissmodel.geo` | Spatial data structures — grid generation, fill strategies, neighborhood |
+| `dissmodel.geo` | Spatial data structures — dual-substrate design (vector + raster) |
 | `dissmodel.models` | Ready-to-use CA and SysDyn reference implementations |
-| `dissmodel.visualization` | Observer-based visualization — `Chart`, `Map`, `display_inputs`, `@track_plot` |
+| `dissmodel.visualization` | Observer-based visualization — `Chart`, `Map`, `RasterMap`, `display_inputs` |
+
+### `dissmodel.geo` — Dual Substrate
+
+The `geo` module provides two independent spatial substrates that share the same simulation clock:
+
+**Vector substrate** (`dissmodel.geo.vector`) — backed by GeoDataFrame. Supports irregular
+geometries, direct GIS integration, and libpysal neighbourhoods (Queen, Rook).
+Use for models that require spatial joins, real-world projections, or interoperability
+with existing GIS workflows.
+
+**Raster substrate** (`dissmodel.geo.raster`) — backed by `RasterBackend` (NumPy 2D arrays).
+Replaces cell-by-cell iteration with fully vectorized operations (`shift2d`, `focal_sum`,
+`neighbor_contact`). At 10,000 cells, the raster substrate is **~4,500× faster** than the
+vector substrate.
+
+| | Vector | Raster |
+|---|---|---|
+| Data structure | GeoDataFrame | NumPy 2D array (`RasterBackend`) |
+| Neighbourhood | Queen / Rook (libpysal) | Moore / Von Neumann (shift2d) |
+| Rule pattern | `rule(idx)` per cell | `rule(arrays) → dict` vectorized |
+| Performance @ 10k cells | ~2,700 ms/step | ~0.6 ms/step |
 
 ---
 
@@ -57,7 +76,7 @@ Chart(show_legend=True)
 env.run(30)
 ```
 
-### Cellular Automaton
+### Cellular Automaton — Vector
 
 ```python
 from dissmodel.core import Environment
@@ -72,6 +91,23 @@ fire.initialize()
 env.run()
 ```
 
+### Cellular Automaton — Raster
+
+```python
+from dissmodel.core import Environment
+from dissmodel.geo.raster.regular_grid import raster_grid
+from dissmodel.models.ca import GameOfLifeRaster
+from dissmodel.visualization.raster_map import RasterMap
+
+backend = raster_grid(rows=100, cols=100, attrs={"state": 0})
+
+env = Environment(end_time=20)
+model = GameOfLifeRaster(backend=backend)
+model.initialize()
+RasterMap(backend=backend, band="state")
+env.run()
+```
+
 ### Streamlit
 
 ```bash
@@ -128,7 +164,9 @@ DisSModel builds on well-established, industry-standard libraries:
 |:---|:---|
 | [GeoPandas](https://geopandas.org/) | Vector data and GeoDataFrame operations |
 | [Salabim](https://www.salabim.org/) | Discrete event simulation engine |
+| [NumPy](https://numpy.org/) | Vectorized array operations — raster substrate |
 | [libpysal](https://pysal.org/libpysal/) | Spatial weights and neighborhood analysis |
+| [rasterio](https://rasterio.readthedocs.io/) | GeoTIFF I/O |
 | [rasterstats](https://pythonhosted.org/rasterstats/) | Raster/vector zonal statistics |
 | [Shapely](https://shapely.readthedocs.io/) | Geometric operations |
 | [Matplotlib](https://matplotlib.org/) | Time-series and spatial visualization |
diff --git a/mkdocs.yml b/mkdocs.yml
index 062e922..00d1b0b 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -66,6 +66,7 @@ nav:
   - API Reference:
     - Core: api/core.md
     - Geo:
+      - Overview:  api/geo/index.md   
       - Vector: api/geo/vector.md
       - Raster:  api/geo/raster.md
     - Visualization: api/visualization.md
diff --git a/paper.md b/paper.md
index ebc40ec..a59e53b 100644
--- a/paper.md
+++ b/paper.md
@@ -31,7 +31,7 @@ bibliography: paper.bib
 
 DisSModel (Discrete Spatial Modeling) is a modular Python framework designed for spatially explicit dynamic modeling, specifically targeting the complexities of Land Use and Land Cover Change (LUCC). Developed to bridge the gap between static geospatial analysis and high-level dynamic simulations, DisSModel translates the modeling paradigms of the TerraME framework [@Carneiro2013] into the Python ecosystem. It enables researchers to simulate complex socio-environmental systems—including forest fires, epidemiological spreads, and urban expansion—by integrating the simulation clock of discrete-event engines with the spatial data structures of GeoPandas [@Jordahl2021].
 
-The framework provides a unified workflow for data ingestion, spatial homogenization via cellular grids, and interactive visualization. DisSModel is currently available for the research community through the LambdaGeo GitHub repository and for testing via TestPyPI.
+The framework provides a **dual-substrate architecture**: a vector substrate backed by GeoDataFrame for flexibility and spatial expressiveness, and a raster substrate backed by NumPy 2D arrays for high-performance vectorized computation. DisSModel is available through the LambdaGeo GitHub repository and on PyPI.
 
 ## Statement of Need
 
@@ -47,7 +47,7 @@ DisSModel occupies a unique niche between general-purpose agent-based modeling (
 |--------|---------|--------------|-----------|
 | Language | Lua | Visual/Internal | Python |
 | Simulation Engine | Discrete Event | Cellular Automata | Integrated Salabim (DES) |
-| Spatial Structure | CellularSpace (Fixed) | Cellular Grid | GeoDataFrame (Dynamic) |
+| Spatial Structure | CellularSpace (Fixed) | Cellular Grid | GeoDataFrame + NumPy (Dual) |
 | GIS Integration | TerraLib | Native Raster | GeoPandas / Rasterio |
 | Extensibility | Script-based | Block-based | Class Inheritance |
 | Anisotropy | GPM Support | Limited | GPM Support |
@@ -56,12 +56,28 @@ While frameworks like **NetLogo** and **Mesa** are excellent for ABM, they often
 
 ## Software Design
 
-DisSModel is organized into four modules, following a strict separation of concerns that allows researchers to extend the framework through class inheritance. 
+DisSModel is organized into four modules, following a strict separation of concerns that allows researchers to extend the framework through class inheritance.
 
-- **Core:** The central engine that manages the simulation clock and discrete-event execution via Salabim integration. 
-- **Geo:** Manages spatial representations. The `vector_grid` function is used for spatial homogenization, aligning disparate rasters and vectors into a unified cellular space.
-- **Models:** Provides templates for common paradigms, including Cellular Automata (CA) and System Dynamics. 
-- **Visualization:** Integrates Matplotlib for static outputs and Streamlit for interactive dashboards.
+- **Core:** The central engine that manages the simulation clock and discrete-event execution via Salabim integration.
+- **Geo:** Manages spatial representations through a dual-substrate design. The vector substrate (`vector_grid`, `SpatialModel`, `CellularAutomaton`) operates on GeoDataFrame structures with libpysal neighborhoods [@Rey2021]. The raster substrate (`raster_grid`, `RasterModel`, `RasterCellularAutomaton`) operates on named NumPy arrays via `RasterBackend`, enabling fully vectorized spatial operations — `shift2d`, `focal_sum`, and `neighbor_contact` — that replace cell-by-cell iteration loops.
+- **Models:** Provides templates for common paradigms, including Cellular Automata (CA) and System Dynamics.
+- **Visualization:** Integrates Matplotlib for static outputs, Streamlit for interactive dashboards, and `RasterMap` for step-by-step raster rendering in both headless and interactive modes.
+
+## Performance
+
+The dual-substrate design exposes a fundamental performance trade-off. The vector substrate offers spatial expressiveness and direct integration with GIS workflows, while the raster substrate achieves high throughput through NumPy vectorization.
+
+Benchmarks running Conway's Game of Life under identical conditions on both substrates show the following results:
+
+| Grid | Cells | Raster (ms/step) | Vector (ms/step) | Speedup |
+|-----:|------:|-----------------:|-----------------:|--------:|
+| 10×10 | 100 | 0.15 | 30.11 | 206× |
+| 50×50 | 2,500 | 0.20 | 647.22 | 3,164× |
+| 100×100 | 10,000 | 0.60 | 2,715.16 | 4,491× |
+| 500×500 | 250,000 | 15.36 | — | — |
+| 1,000×1,000 | 1,000,000 | 25.85 | — | — |
+
+Results were validated cell-by-cell: both substrates produce identical output for all comparable grid sizes. The raster substrate processes grids of one million cells in approximately 26 ms per step, making it suitable for large-scale LUCC simulations such as the BR-MANGUE coastal dynamics model (94,704 cells, 88 time steps).
 
 ## Research impact statement