diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9989028..b451e26 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,13 @@ All notable changes to `disslucc` are documented here.
 
 ## [Unreleased]
 
+### Added
+- `DemandComputeTwoDates` and `DemandComputeThreeDates`: demand computed from the land-use
+  layers (start-year areas and one or two later reference years, interpolated linearly, rounded
+  as LuccME does). They are what LuccME labs 04, 05, 16 and 17 need on top of components that
+  already existed; with them those four labs reproduce the TerraME goldens
+  ([disslucc-benchmark](https://github.com/LambdaGeo/disslucc-benchmark)).
+
 ### Changed
 - `benchmark/goldens/` updated to terrame-docker v0.1.1: same results (CSVs within
   1e-12), with `terrame.log` and `manifest.json` in English; `benchmark/README.md`
diff --git a/docs/api.md b/docs/api.md
index eee8d44..7d2fe91 100644
--- a/docs/api.md
+++ b/docs/api.md
@@ -50,6 +50,45 @@ scripts, no external file needed.
 DemandInline(values: list[list[float]], land_use_types: list[str])
 ```
 
+### `DemandComputeTwoDates` and `DemandComputeThreeDates`
+
+Demand computed from the land-use layers instead of a table (LuccME's
+components of the same names). The start-year area of each class is read
+from the layers at construction; the later area(s) from other columns of the
+layer (for example `f2014`); the years in between are interpolated linearly.
+Both subclass `DemandPreComputedValues`, so they implement the same
+`DemandProtocol`. Create them **before** `env.run()`.
+
+```python
+DemandComputeTwoDates(
+    backend,                          # layers: land_use_types + final_land_use_types (+ optional "mask")
+    land_use_types: list[str],        # e.g. ["f", "d", "outros"], layers at start_year
+    final_land_use_types: list[str],  # e.g. ["f2014", "d2014", "outros"], same order
+    start_year: int, end_year: int,   # rows = end_year - start_year + 1
+    final_year: int,                  # year the final layers describe
+    cell_area: float,                 # unit of the demand
+    direction_for_interpolation: list[str] | None = None,  # "Increase" / "Decrease" per class
+)
+
+DemandComputeThreeDates(
+    backend, land_use_types,
+    middle_land_use_types,            # e.g. ["f2011", "d2011", "outros"]
+    final_land_use_types,
+    start_year, end_year, middle_year, final_year,   # start < middle < final
+    cell_area, direction_for_interpolation=None,
+)
+```
+
+Rounding follows LuccME: each year adds the step to the **already rounded**
+previous demand and rounds half up (`floor(x + 0.5)`); the reference years are
+set to their exact area, rounded. So the demand is integer-valued and can
+differ by a few units from exact linear interpolation. The final year may lie
+after `end_year` (lab17); beyond it the same step keeps being applied. With
+float32 layers the areas carry float32 rounding: use float64 layers when the
+demand must match a reference to the unit. The pure functions
+`interpolate_two_dates` / `interpolate_three_dates` (areas in, rows out) are
+exported from `disslucc.components.demand` too.
+
 ### `load_demand_csv(raw: str, land_use_types: list[str]) -> list[list[float]]`
 
 Parses a CSV with one column per class, one row per step. Column order
diff --git a/src/disslucc/__init__.py b/src/disslucc/__init__.py
index 07d30ee..6fa1296 100644
--- a/src/disslucc/__init__.py
+++ b/src/disslucc/__init__.py
@@ -4,7 +4,13 @@
 (the CLUE algorithm itself is faithful to the original).
 """
 from .components.allocation import AllocationClueLike, AllocationClueLikeSaturation, AllocationDClueSLike
-from .components.demand import DemandInline, DemandPreComputedValues, load_demand_csv
+from .components.demand import (
+    DemandComputeThreeDates,
+    DemandComputeTwoDates,
+    DemandInline,
+    DemandPreComputedValues,
+    load_demand_csv,
+)
 from .components.potential import (
     PotentialDLogisticRegression,
     PotentialLinearRegression,
@@ -28,6 +34,8 @@
     "AllocationClueLikeSaturation",
     "AllocationDClueSLike",
     "AllocationSpec",
+    "DemandComputeThreeDates",
+    "DemandComputeTwoDates",
     "DemandInline",
     "DemandPreComputedValues",
     "DemandProtocol",
diff --git a/src/disslucc/components/demand/__init__.py b/src/disslucc/components/demand/__init__.py
index d7b878e..dcbe7ef 100644
--- a/src/disslucc/components/demand/__init__.py
+++ b/src/disslucc/components/demand/__init__.py
@@ -1,4 +1,18 @@
+from .dates import (
+    DemandComputeThreeDates,
+    DemandComputeTwoDates,
+    interpolate_three_dates,
+    interpolate_two_dates,
+)
 from .inline import DemandInline
 from .precomputed import DemandPreComputedValues, load_demand_csv
 
-__all__ = ["DemandInline", "DemandPreComputedValues", "load_demand_csv"]
+__all__ = [
+    "DemandComputeThreeDates",
+    "DemandComputeTwoDates",
+    "DemandInline",
+    "DemandPreComputedValues",
+    "interpolate_three_dates",
+    "interpolate_two_dates",
+    "load_demand_csv",
+]
diff --git a/src/disslucc/components/demand/dates.py b/src/disslucc/components/demand/dates.py
new file mode 100644
index 0000000..2c43599
--- /dev/null
+++ b/src/disslucc/components/demand/dates.py
@@ -0,0 +1,202 @@
+"""
+disslucc.demand.dates
+----------------------
+Port of LuccME's DemandComputeTwoDates and DemandComputeThreeDates.
+
+Both build the annual demand from the land-use maps themselves instead of a
+table: the initial demand is the area of each class in the layer at the start
+year, the later demands (one or two reference years) are the areas read from
+other columns of the layer (e.g. ``f2011``, ``f2014``), and the years in
+between are interpolated linearly.
+
+The result is handed to :class:`DemandPreComputedValues`, so both classes are
+substrate-neutral from there on and implement the same ``DemandProtocol``.
+
+Rounding follows LuccME exactly (and it matters): every year adds the step to
+the **already rounded** demand of the previous year and rounds again with
+``floor(x + 0.5)``; the reference years are set to their exact value, rounded.
+That is why the demand is integer-valued and can differ by a few units from a
+table computed with exact arithmetic.
+"""
+from __future__ import annotations
+
+import math
+
+import numpy as np
+
+from .precomputed import DemandPreComputedValues
+
+_DIRECTIONS = {"Increase": 1, "Decrease": -1}
+
+
+def _round(x: float) -> int:
+    """LuccME's `math.floor(x + 0.5)` (rounds half up, unlike Python's round())."""
+    return math.floor(x + 0.5)
+
+
+def _directions(
+    direction: list[str] | None, start: list[float], end: list[float]
+) -> list[int]:
+    """+1 / -1 / 0 per class: the given `direction`, or the sign of end - start."""
+    if direction is not None:
+        # anything other than "Increase" / "Decrease" means stable, as in LuccME
+        return [_DIRECTIONS.get(d, 0) for d in direction]
+    return [(e > s) - (e < s) for s, e in zip(start, end)]
+
+
+def interpolate_two_dates(
+    initial: list[float],
+    final: list[float],
+    final_offset: int,
+    n_years: int,
+    direction: list[str] | None = None,
+) -> list[list[int]]:
+    """
+    Annual demand from two dates (the algorithm of ``DemandComputeTwoDates``).
+
+    initial, final : total area per class at the start year and at the final year.
+    final_offset   : final year - start year (> 0).
+    n_years        : number of rows to generate (end year - start year + 1).
+                     Beyond the final year the same step keeps being applied.
+    direction      : optional "Increase" / "Decrease" per class; default is the
+                     sign of ``final - initial``.
+    """
+    if final_offset <= 0:
+        raise ValueError("final_year must be after start_year")
+    sign = _directions(direction, initial, final)
+    step = [abs(f - i) / final_offset for i, f in zip(initial, final)]
+    rows = [[_round(v) for v in initial]]
+    for year in range(1, n_years):                      # year = offset from the start
+        if year == final_offset:
+            rows.append([_round(v) for v in final])
+        else:
+            prev = rows[-1]
+            rows.append([_round(prev[j] + step[j] * sign[j]) for j in range(len(initial))])
+    return rows[:n_years]
+
+
+def interpolate_three_dates(
+    initial: list[float],
+    middle: list[float],
+    final: list[float],
+    middle_offset: int,
+    final_offset: int,
+    n_years: int,
+    direction: list[str] | None = None,
+) -> list[list[int]]:
+    """
+    Annual demand from three dates (``DemandComputeThreeDates``): linear from the
+    start to the middle year, then linear from the middle to the final year, each
+    segment with its own step and direction. The final year may lie after the end
+    of the simulation (then only the first segment, or part of the second, is used).
+    """
+    if not 0 < middle_offset < final_offset:
+        raise ValueError("start_year < middle_year < final_year is required")
+    seg1, seg2 = middle_offset, final_offset - middle_offset
+    sign1, sign2 = _directions(direction, initial, middle), _directions(direction, middle, final)
+    step1 = [abs(m - i) / seg1 for i, m in zip(initial, middle)]
+    step2 = [abs(f - m) / seg2 for m, f in zip(middle, final)]
+    rows = [[_round(v) for v in initial]]
+    for year in range(1, n_years):
+        if year == middle_offset:
+            rows.append([_round(v) for v in middle])
+        elif year == final_offset:
+            rows.append([_round(v) for v in final])
+        else:
+            step, sign = (step1, sign1) if year < middle_offset else (step2, sign2)
+            prev = rows[-1]
+            rows.append([_round(prev[j] + step[j] * sign[j]) for j in range(len(initial))])
+    return rows[:n_years]
+
+
+def _areas(backend, names: list[str], cell_area: float) -> list[float]:
+    """Total area of each named layer over the valid cells (mask), cell by cell as LuccME does."""
+    first = np.asarray(backend.get(names[0]))
+    mask = np.asarray(backend.arrays.get("mask", np.ones(first.shape, dtype=bool))).astype(bool)
+    return [float((np.asarray(backend.get(n), dtype=np.float64)[mask] * cell_area).sum()) for n in names]
+
+
+class DemandComputeTwoDates(DemandPreComputedValues):
+    """
+    Demand interpolated between the start year and one final year.
+
+    Parameters (setup)
+    -------------------
+    backend                  : substrate with the class layers and the final-year layers.
+    land_use_types           : class layers at the start year, e.g. ["f", "d", "outros"].
+    final_land_use_types     : layers with the same classes at the final year, in the
+                               same order, e.g. ["f2014", "d2014", "outros"].
+    start_year, end_year     : first and last simulated year (rows = end - start + 1).
+    final_year               : year the ``final_land_use_types`` layers describe.
+    cell_area                : area of one cell (the unit of the demand).
+    direction_for_interpolation : optional "Increase" / "Decrease" per class.
+
+    Create it before ``env.run()``: the start-year areas are read at construction.
+    With float32 layers the areas carry float32 rounding; use float64 layers when
+    the demand must match a reference to the unit.
+    """
+
+    def setup(
+        self,
+        backend,
+        land_use_types: list[str],
+        final_land_use_types: list[str],
+        start_year: int,
+        end_year: int,
+        final_year: int,
+        cell_area: float,
+        direction_for_interpolation: list[str] | None = None,
+    ) -> None:
+        if len(final_land_use_types) != len(land_use_types):
+            raise ValueError("final_land_use_types must have one layer per land use type")
+        if direction_for_interpolation is not None and len(direction_for_interpolation) != len(land_use_types):
+            raise ValueError("direction_for_interpolation must have one entry per land use type")
+        rows = interpolate_two_dates(
+            _areas(backend, land_use_types, cell_area),
+            _areas(backend, final_land_use_types, cell_area),
+            final_year - start_year,
+            end_year - start_year + 1,
+            direction_for_interpolation,
+        )
+        super().setup(annual_demand=rows, land_use_types=land_use_types)
+
+
+class DemandComputeThreeDates(DemandPreComputedValues):
+    """
+    Demand interpolated through a middle year to a final year (see
+    :class:`DemandComputeTwoDates` for the shared parameters).
+
+    Extra parameters (setup)
+    -------------------------
+    middle_land_use_types : layers at the middle year, e.g. ["f2011", "d2011", "outros"].
+    middle_year           : year they describe; start_year < middle_year < final_year.
+    """
+
+    def setup(
+        self,
+        backend,
+        land_use_types: list[str],
+        middle_land_use_types: list[str],
+        final_land_use_types: list[str],
+        start_year: int,
+        end_year: int,
+        middle_year: int,
+        final_year: int,
+        cell_area: float,
+        direction_for_interpolation: list[str] | None = None,
+    ) -> None:
+        n = len(land_use_types)
+        if len(middle_land_use_types) != n or len(final_land_use_types) != n:
+            raise ValueError("middle and final land use types must have one layer per land use type")
+        if direction_for_interpolation is not None and len(direction_for_interpolation) != n:
+            raise ValueError("direction_for_interpolation must have one entry per land use type")
+        rows = interpolate_three_dates(
+            _areas(backend, land_use_types, cell_area),
+            _areas(backend, middle_land_use_types, cell_area),
+            _areas(backend, final_land_use_types, cell_area),
+            middle_year - start_year,
+            final_year - start_year,
+            end_year - start_year + 1,
+            direction_for_interpolation,
+        )
+        super().setup(annual_demand=rows, land_use_types=land_use_types)
diff --git a/tests/test_demand_dates.py b/tests/test_demand_dates.py
new file mode 100644
index 0000000..dc64190
--- /dev/null
+++ b/tests/test_demand_dates.py
@@ -0,0 +1,203 @@
+"""
+tests/test_demand_dates.py
+---------------------------
+DemandComputeTwoDates / DemandComputeThreeDates (LuccME demand interpolated
+from the land-use layers).
+
+The expected tables below are what these components produce on the layers in
+data/input; with them lab04, lab05, lab16 and lab17 reproduce the TerraME
+goldens (disslucc-benchmark), so they are the demand LuccME computes.
+"""
+from pathlib import Path
+
+import geopandas as gpd
+import numpy as np
+import pytest
+from dissmodel.core import Environment
+from dissmodel.geo.raster.backend import RasterBackend
+
+from disslucc import DemandComputeThreeDates, DemandComputeTwoDates
+from disslucc.components.demand.dates import interpolate_three_dates, interpolate_two_dates
+
+DATA = Path(__file__).resolve().parent.parent / "data" / "input"
+
+
+# ── the interpolation itself ─────────────────────────────────────────────────────────────────
+
+def test_rounding_accumulates_on_the_rounded_previous_year():
+    # LuccME adds the step to the already rounded demand: 10 + 0.4 -> 10 every year,
+    # then the final year is set exactly. Exact linear interpolation would give 10, 11, 11, 12.
+    assert interpolate_two_dates([10], [12], final_offset=5, n_years=6) == [[10], [10], [10], [10], [10], [12]]
+
+
+def test_rounds_half_up():
+    # step 0.5: 0 -> floor(0.5 + 0.5) = 1 (Python's round() would give 0)
+    assert interpolate_two_dates([0], [1], final_offset=2, n_years=3) == [[0], [1], [1]]
+
+
+def test_direction_from_the_data_and_static_class():
+    rows = interpolate_two_dates([100, 50, 7], [110, 40, 7], final_offset=2, n_years=3)
+    assert rows == [[100, 50, 7], [105, 45, 7], [110, 40, 7]]
+
+
+def test_given_direction_overrides_the_sign():
+    # "Decrease" on a growing class steps down, then the final year is set exactly (as in LuccME)
+    rows = interpolate_two_dates([100], [110], final_offset=2, n_years=3, direction=["Decrease"])
+    assert rows == [[100], [95], [110]]
+    # anything else than Increase / Decrease means stable
+    rows = interpolate_two_dates([100], [110], final_offset=2, n_years=3, direction=["Stable"])
+    assert rows == [[100], [100], [110]]
+
+
+def test_keeps_the_same_step_after_the_final_year():
+    assert interpolate_two_dates([0], [10], final_offset=2, n_years=5) == [[0], [5], [10], [15], [20]]
+
+
+def test_stops_before_the_final_year():
+    assert interpolate_two_dates([0], [10], final_offset=5, n_years=3) == [[0], [2], [4]]
+
+
+def test_three_dates_two_segments():
+    rows = interpolate_three_dates([0], [10], [14], middle_offset=2, final_offset=4, n_years=5)
+    assert rows == [[0], [5], [10], [12], [14]]
+
+
+def test_three_dates_final_year_after_the_end_of_the_run():
+    # lab17: the final year (2007) is later than the last simulated year (2004)
+    rows = interpolate_three_dates([0], [10], [14], middle_offset=2, final_offset=4, n_years=3)
+    assert rows == [[0], [5], [10]]
+
+
+@pytest.mark.parametrize("call", [
+    lambda: interpolate_two_dates([1], [2], final_offset=0, n_years=3),
+    lambda: interpolate_three_dates([1], [2], [3], middle_offset=3, final_offset=3, n_years=5),
+    lambda: interpolate_three_dates([1], [2], [3], middle_offset=0, final_offset=3, n_years=5),
+])
+def test_invalid_years_are_rejected(call):
+    with pytest.raises(ValueError):
+        call()
+
+
+# ── the components: areas read from the layers ───────────────────────────────────────────────
+
+def _backend(layers: dict, mask=None) -> RasterBackend:
+    shape = np.asarray(next(iter(layers.values()))).shape
+    backend = RasterBackend(shape=shape)
+    backend.set("mask", np.ones(shape, bool) if mask is None else mask)
+    for name, arr in layers.items():
+        backend.set(name, np.asarray(arr, dtype=np.float64))
+    return backend
+
+
+def test_two_dates_sums_cell_area_over_valid_cells_only():
+    mask = np.array([[True, True], [True, False]])
+    backend = _backend({
+        "f": [[1, 1], [0, 1]], "d": [[0, 0], [1, 0]],     # the unmasked cell is ignored
+        "f9": [[1, 0], [0, 1]], "d9": [[0, 1], [1, 0]],
+    }, mask)
+    Environment(end_time=2)
+    demand = DemandComputeTwoDates(
+        backend=backend, land_use_types=["f", "d"], final_land_use_types=["f9", "d9"],
+        start_year=2000, end_year=2002, final_year=2002, cell_area=25.0,
+    )
+    # start: f = 2 cells (50), d = 1 cell (25); final: f = 1 cell (25), d = 2 cells (50); step 12.5,
+    # applied to the rounded previous year: 50 - 12.5 = 37.5 -> 38 and 25 + 12.5 = 37.5 -> 38
+    assert demand.annual_demand == [[50, 25], [38, 38], [25, 50]]
+
+
+def test_the_mask_matters():
+    layers = {"f": [[1, 1], [0, 1]], "f9": [[1, 0], [0, 1]]}
+    Environment(end_time=1)
+    full = DemandComputeTwoDates(backend=_backend(layers), land_use_types=["f"], final_land_use_types=["f9"],
+                                 start_year=0, end_year=1, final_year=1, cell_area=1.0)
+    assert full.annual_demand == [[3], [2]]
+
+
+def test_component_exposes_the_demand_protocol():
+    backend = _backend({"f": [[4]], "d": [[0]], "f9": [[2]], "d9": [[2]]})
+    env = Environment(end_time=1)
+    demand = DemandComputeTwoDates(backend=backend, land_use_types=["f", "d"], final_land_use_types=["f9", "d9"],
+                                   start_year=0, end_year=1, final_year=1, cell_area=1.0)
+    env._now = 1
+    demand.execute()
+    assert demand.get_current_lu_demand(0) == 2 and demand.get_previous_lu_demand(0) == 4
+    assert demand.get_current_lu_direction(0) == -1 and demand.get_current_lu_direction(1) == 1
+
+
+def test_mismatched_layer_lists_are_rejected():
+    backend = _backend({"f": [[1]], "d": [[1]], "f9": [[1]]})
+    Environment(end_time=1)
+    with pytest.raises(ValueError):
+        DemandComputeTwoDates(backend=backend, land_use_types=["f", "d"], final_land_use_types=["f9"],
+                              start_year=0, end_year=1, final_year=1, cell_area=1.0)
+    with pytest.raises(ValueError):
+        DemandComputeThreeDates(backend=backend, land_use_types=["f"], middle_land_use_types=["f9"],
+                                final_land_use_types=["f9"], start_year=0, end_year=2, middle_year=1, final_year=2,
+                                cell_area=1.0, direction_for_interpolation=["Increase", "Decrease"])
+
+
+# ── LuccME labs 04, 05, 16, 17 on the layers in data/input ───────────────────────────────────
+
+def _layer_backend(name: str, columns: list[str]) -> RasterBackend:
+    cells = gpd.read_file(DATA / f"{name}.zip")
+    rows = cells["row" if "row" in cells else "lin"].astype(int).values
+    cols = cells["col"].astype(int).values
+    shape = (rows.max() + 1, cols.max() + 1)
+    backend = RasterBackend(shape=shape)
+    mask = np.zeros(shape, bool)
+    mask[rows, cols] = True
+    backend.set("mask", mask)
+    for c in columns:
+        arr = np.zeros(shape)
+        arr[rows, cols] = cells[c].astype(float).values
+        backend.set(c, arr)
+    return backend
+
+
+CSAC = ["f", "d", "outros", "f2011", "d2011", "f2014", "d2014"]
+MOJU = ["f", "d", "o", "f04", "d04", "f07", "d07"]
+
+
+def test_lab04_two_dates_csac():
+    backend = _layer_backend("csAC", CSAC)
+    Environment(end_time=6)
+    demand = DemandComputeTwoDates(
+        backend=backend, land_use_types=["f", "d", "outros"], final_land_use_types=["f2014", "d2014", "outros"],
+        start_year=2008, end_year=2014, final_year=2014, cell_area=25.0)
+    assert demand.annual_demand == [
+        [137878, 19983, 6489], [137607, 20254, 6489], [137336, 20525, 6489], [137065, 20796, 6489],
+        [136794, 21067, 6489], [136523, 21338, 6489], [136253, 21607, 6489]]
+
+
+def test_lab05_three_dates_csac():
+    backend = _layer_backend("csAC", CSAC)
+    Environment(end_time=6)
+    demand = DemandComputeThreeDates(
+        backend=backend, land_use_types=["f", "d", "outros"], middle_land_use_types=["f2011", "d2011", "outros"],
+        final_land_use_types=["f2014", "d2014", "outros"], start_year=2008, end_year=2014, middle_year=2011,
+        final_year=2014, cell_area=25.0)
+    assert demand.annual_demand == [
+        [137878, 19983, 6489], [137622, 20239, 6489], [137366, 20495, 6489], [137110, 20750, 6489],
+        [136824, 21036, 6489], [136538, 21322, 6489], [136253, 21607, 6489]]
+
+
+MOJU_TABLE = [[5706, 205, 3], [5658, 253, 3], [5610, 301, 3], [5562, 349, 3], [5514, 397, 3], [5468, 443, 3]]
+
+
+def test_lab16_two_dates_moju():
+    backend = _layer_backend("cs_moju", MOJU)
+    Environment(end_time=5)
+    demand = DemandComputeTwoDates(
+        backend=backend, land_use_types=["f", "d", "o"], final_land_use_types=["f04", "d04", "o"],
+        start_year=1999, end_year=2004, final_year=2004, cell_area=1.0)
+    assert demand.annual_demand == MOJU_TABLE
+
+
+def test_lab17_three_dates_moju_final_year_after_the_run():
+    backend = _layer_backend("cs_moju", MOJU)
+    Environment(end_time=5)
+    demand = DemandComputeThreeDates(
+        backend=backend, land_use_types=["f", "d", "o"], middle_land_use_types=["f04", "d04", "o"],
+        final_land_use_types=["f07", "d07", "o"], start_year=1999, end_year=2004, middle_year=2004,
+        final_year=2007, cell_area=1.0)
+    assert demand.annual_demand == MOJU_TABLE