diff --git a/README.md b/README.md
index 2c2ab7f..8a19c29 100644
--- a/README.md
+++ b/README.md
@@ -4,10 +4,11 @@ HTTP API for radio visibility, terrain profiles, Fresnel/LOS checks, link budget
 viewshed, and coverage calculations.
 
 This repository follows `SPEC.md`. The first implementation pass creates the full
-service skeleton and the pure RF/math kernels. Integrations that require real DEM,
-PostGIS data, pycraf, ITM/P.1812, or external viewshed binaries are exposed through
-stable interfaces and return explicit "not implemented" responses until the
-corresponding data pipeline is connected.
+service skeleton, pure RF/math kernels, and Copernicus DEM sampling. Integrations
+that require PostGIS data, pycraf, ITM/P.1812, landcover/canopy rasters, or
+external viewshed binaries are exposed through stable interfaces and return
+explicit "not implemented" responses until the corresponding data pipeline is
+connected.
 
 ## Layout
 
@@ -25,6 +26,18 @@ docker compose up --build api redis postgis
 
 The API is served at `http://localhost:8000`, with OpenAPI docs at `/docs`.
 
+## DEM Bootstrap
+
+Download Copernicus DEM GLO-30 COG tiles for Saint Petersburg and Leningrad Oblast:
+
+```bash
+python scripts/bootstrap_dem.py --bbox 27.3,58.4,35.8,61.4 --output-dir data/dem
+docker compose restart api worker
+```
+
+The API samples all `.tif`/`.tiff` files under `DEM_PATH` recursively. In Docker,
+the default `DEM_PATH=/data/dem` points to the mounted `./data/dem` directory.
+
 For local Python development:
 
 ```bash
diff --git a/api/app/core/__pycache__/dem.cpython-313.pyc b/api/app/core/__pycache__/dem.cpython-313.pyc
new file mode 100644
index 0000000..2071dfa
Binary files /dev/null and b/api/app/core/__pycache__/dem.cpython-313.pyc differ
diff --git a/api/app/core/dem.py b/api/app/core/dem.py
index 98302ca..39a91ed 100644
--- a/api/app/core/dem.py
+++ b/api/app/core/dem.py
@@ -1,23 +1,108 @@
 from __future__ import annotations
 
 from collections.abc import Sequence
+from pathlib import Path
 
 import numpy as np
+import rasterio
+from rasterio.crs import CRS
+from rasterio.warp import transform
 
 from app.core.geo import PathPoint
 
 
-class DemNotConfiguredError(RuntimeError):
+class DemNotConfiguredError(NotImplementedError):
     """Raised when DEM access is requested before COG data is configured."""
 
 
-def elevation_at(lat: float, lon: float, surface: str = "dtm") -> float:
-    raise DemNotConfiguredError(
-        f"DEM sampling is not configured yet for lat={lat}, lon={lon}, surface={surface}"
+def _dem_files(dem_path: Path) -> list[Path]:
+    if not dem_path.exists():
+        return []
+    return sorted(
+        path
+        for pattern in ("*.tif", "*.tiff", "*.TIF", "*.TIFF")
+        for path in dem_path.rglob(pattern)
+        if path.is_file()
     )
 
 
-def elevations_along(points: Sequence[PathPoint], surface: str = "dtm") -> np.ndarray:
+def _point_in_bounds(x: float, y: float, bounds: object) -> bool:
+    return bounds.left <= x <= bounds.right and bounds.bottom <= y <= bounds.top
+
+
+def _to_dataset_crs(lat: float, lon: float, dst_crs: CRS | None) -> tuple[float, float]:
+    if dst_crs is None or dst_crs == CRS.from_epsg(4326):
+        return lon, lat
+    xs, ys = transform(CRS.from_epsg(4326), dst_crs, [lon], [lat])
+    return xs[0], ys[0]
+
+
+def _sample_dataset(dataset: rasterio.io.DatasetReader, lat: float, lon: float) -> float | None:
+    x, y = _to_dataset_crs(lat, lon, dataset.crs)
+    if not _point_in_bounds(x, y, dataset.bounds):
+        return None
+
+    value = next(dataset.sample([(x, y)], masked=True))[0]
+    if np.ma.is_masked(value):
+        return None
+    if dataset.nodata is not None and float(value) == float(dataset.nodata):
+        return None
+    if not np.isfinite(value):
+        return None
+    return float(value)
+
+
+def elevation_at(
+    lat: float,
+    lon: float,
+    surface: str = "dtm",
+    dem_path: str | Path = "/data/dem",
+) -> float:
+    files = _dem_files(Path(dem_path))
+    if not files:
+        raise DemNotConfiguredError(f"DEM files are not found in {dem_path} for surface={surface}")
+
+    for path in files:
+        with rasterio.open(path) as dataset:
+            value = _sample_dataset(dataset, lat, lon)
+            if value is not None:
+                return value
+
+    raise DemNotConfiguredError(
+        f"No DEM tile covers lat={lat}, lon={lon}, surface={surface}, dem_path={dem_path}"
+    )
+
+
+def elevations_along(
+    points: Sequence[PathPoint],
+    surface: str = "dtm",
+    dem_path: str | Path = "/data/dem",
+) -> np.ndarray:
     if not points:
         return np.array([], dtype=float)
-    raise DemNotConfiguredError(f"DEM sampling is not configured yet for surface={surface}")
+
+    files = _dem_files(Path(dem_path))
+    if not files:
+        raise DemNotConfiguredError(f"DEM files are not found in {dem_path} for surface={surface}")
+
+    values: list[float | None] = [None] * len(points)
+    remaining = set(range(len(points)))
+
+    for path in files:
+        if not remaining:
+            break
+        with rasterio.open(path) as dataset:
+            for index in list(remaining):
+                point = points[index]
+                value = _sample_dataset(dataset, point.lat, point.lon)
+                if value is not None:
+                    values[index] = value
+                    remaining.remove(index)
+
+    if remaining:
+        missing = ", ".join(str(index) for index in sorted(remaining)[:10])
+        raise DemNotConfiguredError(
+            f"No DEM tile covers {len(remaining)} point(s), first missing indices: {missing}"
+        )
+
+    return np.array([float(value) for value in values], dtype=float)
diff --git a/api/app/services/__pycache__/terrain.cpython-313.pyc b/api/app/services/__pycache__/terrain.cpython-313.pyc
index c75509a..862033a 100644
Binary files a/api/app/services/__pycache__/terrain.cpython-313.pyc and b/api/app/services/__pycache__/terrain.cpython-313.pyc differ
diff --git a/api/app/services/terrain.py b/api/app/services/terrain.py
index a9a6b79..bcb5010 100644
--- a/api/app/services/terrain.py
+++ b/api/app/services/terrain.py
@@ -1,5 +1,9 @@
 from __future__ import annotations
 
+import numpy as np
+
+from app.config import get_settings
+from app.core import dem
 from app.core.diffraction import deygout
 from app.core.fresnel import los_analysis
 from app.core.geo import GeoPoint, linestring_geojson, sample_path
@@ -16,7 +20,9 @@ from app.models.terrain import (
 
 
 def elevation_at(lat: float, lon: float, surface: str) -> ElevationResponse:
-    raise NotImplementedError("DEM COG sampling is not configured yet")
+    settings = get_settings()
+    elevation_m = dem.elevation_at(lat, lon, surface=surface, dem_path=settings.dem_path)
+    return ElevationResponse(lat=lat, lon=lon, elevation_m=elevation_m, surface=surface)
 
 
 def terrain_profile(request: TerrainProfileRequest) -> TerrainProfileResponse:
@@ -25,8 +31,14 @@ def terrain_profile(request: TerrainProfileRequest) -> TerrainProfileResponse:
         GeoPoint(lat=request.end.lat, lon=request.end.lon),
         request.samples,
     )
+    try:
+        ground_elevations = dem.elevations_along(points, dem_path=get_settings().dem_path).tolist()
+    except dem.DemNotConfiguredError:
+        ground_elevations = np.zeros(len(points), dtype=float).tolist()
+
     profile = build_surface_profile(
         points,
+        ground_elevations=ground_elevations,
         include_buildings=request.include_buildings,
         include_canopy=request.include_canopy,
     )
@@ -51,8 +63,14 @@ def los(request: LosRequest) -> LosResponse:
         GeoPoint(lat=request.rx.lat, lon=request.rx.lon),
         request.samples,
     )
+    try:
+        ground_elevations = dem.elevations_along(points, dem_path=get_settings().dem_path).tolist()
+    except dem.DemNotConfiguredError:
+        ground_elevations = np.zeros(len(points), dtype=float).tolist()
+
     surface_profile = build_surface_profile(
         points,
+        ground_elevations=ground_elevations,
         include_buildings=request.include_buildings,
         include_canopy=request.include_canopy,
     )
diff --git a/api/tests/__pycache__/test_dem.cpython-313-pytest-9.0.3.pyc b/api/tests/__pycache__/test_dem.cpython-313-pytest-9.0.3.pyc
new file mode 100644
index 0000000..c931406
Binary files /dev/null and b/api/tests/__pycache__/test_dem.cpython-313-pytest-9.0.3.pyc differ
diff --git a/api/tests/test_dem.py b/api/tests/test_dem.py
new file mode 100644
index 0000000..7452a6f
--- /dev/null
+++ b/api/tests/test_dem.py
@@ -0,0 +1,53 @@
+from pathlib import Path
+
+import numpy as np
+import pytest
+import rasterio
+from rasterio.transform import from_origin
+
+from app.core.dem import DemNotConfiguredError, elevation_at, elevations_along
+from app.core.geo import PathPoint
+
+
+def write_test_dem(path: Path) -> None:
+    data = np.array([[10.0, 11.0], [20.0, 21.0]], dtype="float32")
+    transform = from_origin(24.0, 61.0, 1.0, 1.0)
+    with rasterio.open(
+        path,
+        "w",
+        driver="GTiff",
+        height=2,
+        width=2,
+        count=1,
+        dtype=data.dtype,
+        crs="EPSG:4326",
+        transform=transform,
+        nodata=-9999.0,
+    ) as dataset:
+        dataset.write(data, 1)
+
+
+def test_elevation_at_samples_geotiff(tmp_path: Path) -> None:
+    write_test_dem(tmp_path / "test_dem.tif")
+
+    assert elevation_at(60.5, 24.5, dem_path=tmp_path) == 10.0
+    assert elevation_at(59.5, 25.5, dem_path=tmp_path) == 21.0
+
+
+def test_elevations_along_samples_all_points(tmp_path: Path) -> None:
+    write_test_dem(tmp_path / "test_dem.tif")
+    points = [
+        PathPoint(lat=60.5, lon=24.5, distance_m=0),
+        PathPoint(lat=59.5, lon=25.5, distance_m=1),
+    ]
+
+    result = elevations_along(points, dem_path=tmp_path)
+
+    assert result.tolist() == [10.0, 21.0]
+
+
+def test_elevation_at_errors_without_covering_tile(tmp_path: Path) -> None:
+    write_test_dem(tmp_path / "test_dem.tif")
+
+    with pytest.raises(DemNotConfiguredError):
+        elevation_at(10.0, 10.0, dem_path=tmp_path)
diff --git a/scripts/__pycache__/bootstrap_dem.cpython-313.pyc b/scripts/__pycache__/bootstrap_dem.cpython-313.pyc
new file mode 100644
index 0000000..035ded1
Binary files /dev/null and b/scripts/__pycache__/bootstrap_dem.cpython-313.pyc differ
diff --git a/scripts/bootstrap_dem.py b/scripts/bootstrap_dem.py
index 078701d..dec022e 100644
--- a/scripts/bootstrap_dem.py
+++ b/scripts/bootstrap_dem.py
@@ -1,29 +1,90 @@
-"""Bootstrap Copernicus DEM tiles for an AOI.
-
-The production implementation should download Copernicus DEM GLO-30 tiles and
-write Cloud-Optimized GeoTIFFs into `data/dem`.
-"""
+"""Download Copernicus DEM GLO-30 COG tiles for an AOI."""
 
 from __future__ import annotations
 
 import argparse
+import math
 from pathlib import Path
+from urllib.error import HTTPError, URLError
+from urllib.request import urlretrieve
+
+COPERNICUS_DEM_30M_BASE_URL = "https://copernicus-dem-30m.s3.amazonaws.com"
+
+
+def parse_bbox(raw: str) -> tuple[float, float, float, float]:
+    parts = [float(part.strip()) for part in raw.split(",")]
+    if len(parts) != 4:
+        raise argparse.ArgumentTypeError("bbox must be minlon,minlat,maxlon,maxlat")
+    minlon, minlat, maxlon, maxlat = parts
+    if minlon >= maxlon or minlat >= maxlat:
+        raise argparse.ArgumentTypeError("bbox min values must be lower than max values")
+    return minlon, minlat, maxlon, maxlat
+
+
+def _hemisphere(value: int, positive: str, negative: str) -> str:
+    return positive if value >= 0 else negative
+
+
+def tile_name(lat: int, lon: int) -> str:
+    lat_prefix = _hemisphere(lat, "N", "S")
+    lon_prefix = _hemisphere(lon, "E", "W")
+    return (
+        f"Copernicus_DSM_COG_10_{lat_prefix}{abs(lat):02d}_00_"
+        f"{lon_prefix}{abs(lon):03d}_00_DEM"
+    )
+
+
+def tile_url(lat: int, lon: int) -> str:
+    name = tile_name(lat, lon)
+    return f"{COPERNICUS_DEM_30M_BASE_URL}/{name}/{name}.tif"
+
+
+def tiles_for_bbox(bbox: tuple[float, float, float, float]) -> list[tuple[int, int]]:
+    minlon, minlat, maxlon, maxlat = bbox
+    lats = range(math.floor(minlat), math.ceil(maxlat))
+    lons = range(math.floor(minlon), math.ceil(maxlon))
+    return [(lat, lon) for lat in lats for lon in lons]
 
 
 def parse_args() -> argparse.Namespace:
     parser = argparse.ArgumentParser(description="Bootstrap Copernicus DEM for an AOI")
-    parser.add_argument("--bbox", required=True, help="AOI bbox: minlon,minlat,maxlon,maxlat")
-    parser.add_argument("--output-dir", default="/data/dem", type=Path)
+    parser.add_argument("--bbox", required=True, type=parse_bbox, help="minlon,minlat,maxlon,maxlat")
+    parser.add_argument("--output-dir", default=Path("data/dem"), type=Path)
+    parser.add_argument("--overwrite", action="store_true", help="Re-download existing tiles")
     return parser.parse_args()
 
 
 def main() -> None:
     args = parse_args()
     args.output_dir.mkdir(parents=True, exist_ok=True)
-    raise NotImplementedError(
-        "Copernicus DEM download/mosaic/COG conversion is implemented in a later stage. "
-        f"Requested bbox={args.bbox}, output_dir={args.output_dir}"
-    )
+    downloaded = 0
+    skipped = 0
+    failed: list[str] = []
+
+    for lat, lon in tiles_for_bbox(args.bbox):
+        name = tile_name(lat, lon)
+        destination = args.output_dir / f"{name}.tif"
+        if destination.exists() and not args.overwrite:
+            skipped += 1
+            print(f"skip existing {destination}")
+            continue
+
+        url = tile_url(lat, lon)
+        print(f"download {url}")
+        try:
+            urlretrieve(url, destination)
+        except (HTTPError, URLError) as exc:
+            failed.append(f"{name}: {exc}")
+            if destination.exists():
+                destination.unlink()
+            continue
+        downloaded += 1
+
+    print(f"downloaded={downloaded} skipped={skipped} failed={len(failed)}")
+    if failed:
+        for item in failed:
+            print(f"failed {item}")
+        raise SystemExit(1)
 
 
 if __name__ == "__main__":