cpp port 1/2

cpp/.clangd

@@ -0,0 +1,5 @@
+# Point clangd to the CMake-generated compile database.
+# Run `cmake --preset vcpkg` (Ninja) first to generate build/compile_commands.json.
+# The vcpkg-win (Visual Studio) generator does NOT produce this file.
+CompileFlags:
+  CompilationDatabase: build

cpp/CMakePresets.json

@@ -13,7 +13,8 @@
             "binaryDir": "${sourceDir}/build",
             "cacheVariables": {
                 "CMAKE_TOOLCHAIN_FILE": "$env{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake",
-                "CMAKE_BUILD_TYPE": "Release"
+                "CMAKE_BUILD_TYPE": "Release",
+                "CMAKE_EXPORT_COMPILE_COMMANDS": "ON"
             }
         },
         {

cpp/src/geo_merge.cpp

@@ -1,28 +1,206 @@
 #include "geo_merge.h"
 
 #include <stdexcept>
+#include <iostream>
 
-// GEOS C API
+// GEOS C API — robust computational geometry with snap rounding.
+// GEOSUnaryUnion uses noded overlay which preserves full double
+// precision (no coordinate snapping), unlike polygon-clipping.js
+// which uses integer grid approximation.
 #include <geos_c.h>
 
+// OGR for GeoJSON ↔ GEOS roundtrip
+#include <ogr_api.h>
+#include <gdal.h>
+#include <cpl_conv.h>
+
 namespace geo_merge {
 
+// ── RAII wrapper for GEOS context (thread-safe) ─────────────────────────────
+
+struct GEOSCtx {
+    GEOSContextHandle_t ctx = nullptr;
+
+    GEOSCtx() {
+        ctx = GEOS_init_r();
+        // Suppress GEOS notices (they're noisy on valid-but-complex geom)
+        GEOSContext_setNoticeHandler_r(ctx, [](const char*, ...) {});
+        GEOSContext_setErrorHandler_r(ctx, [](const char* fmt, ...) {
+            // Optionally log GEOS errors for debugging
+            va_list ap;
+            va_start(ap, fmt);
+            char buf[512];
+            vsnprintf(buf, sizeof(buf), fmt, ap);
+            va_end(ap);
+            std::cerr << "[GEOS] " << buf << "\n";
+        });
+    }
+
+    ~GEOSCtx() {
+        if (ctx) GEOS_finish_r(ctx);
+    }
+
+    GEOSCtx(const GEOSCtx&) = delete;
+    GEOSCtx& operator=(const GEOSCtx&) = delete;
+};
+
+// ── GeoJSON → GEOS geometry ─────────────────────────────────────────────────
+
+static GEOSGeometry* geojson_to_geos(
+    GEOSContextHandle_t ctx,
+    const nlohmann::json& geojson
+) {
+    std::string str = geojson.dump();
+    OGRGeometryH ogr = OGR_G_CreateGeometryFromJson(str.c_str());
+    if (!ogr) return nullptr;
+
+    // OGR → WKB → GEOS (lossless roundtrip, max precision)
+    int wkbSize = OGR_G_WkbSize(ogr);
+    std::vector<unsigned char> wkb(wkbSize);
+    OGR_G_ExportToWkb(ogr, wkbNDR, wkb.data());
+    OGR_G_DestroyGeometry(ogr);
+
+    GEOSWKBReader* reader = GEOSWKBReader_create_r(ctx);
+    GEOSGeometry* geom = GEOSWKBReader_read_r(ctx, reader, wkb.data(), wkb.size());
+    GEOSWKBReader_destroy_r(ctx, reader);
+
+    return geom;
+}
+
+// ── GEOS geometry → GeoJSON ─────────────────────────────────────────────────
+
+static nlohmann::json geos_to_geojson(
+    GEOSContextHandle_t ctx,
+    const GEOSGeometry* geom
+) {
+    // GEOS → WKB → OGR → GeoJSON (lossless)
+    GEOSWKBWriter* writer = GEOSWKBWriter_create_r(ctx);
+    size_t wkbSize = 0;
+    unsigned char* wkb = GEOSWKBWriter_write_r(ctx, writer, geom, &wkbSize);
+    GEOSWKBWriter_destroy_r(ctx, writer);
+
+    OGRGeometryH ogr = nullptr;
+    unsigned char* wkbPtr = wkb;
+    OGR_G_CreateFromWkb(wkbPtr, nullptr, &ogr, static_cast<int>(wkbSize));
+    GEOSFree_r(ctx, wkb);
+
+    if (!ogr) return nlohmann::json::object();
+
+    char* json = OGR_G_ExportToJson(ogr);
+    OGR_G_DestroyGeometry(ogr);
+
+    auto result = nlohmann::json::parse(json);
+    CPLFree(json);
+
+    return result;
+}
+
+// ── Public API ──────────────────────────────────────────────────────────────
+
 std::vector<boundary::BoundaryFeature> merge(
     const std::vector<boundary::BoundaryFeature>& features
 ) {
-    // TODO: Implement GEOS-based geometry union
-    //
-    // For each group of features sharing a GID code:
-    //   1. Collect GEOS geometry handles
-    //   2. Create a GeometryCollection
-    //   3. GEOSUnaryUnion → merged geometry
-    //   4. On failure: GEOSMakeValid each input, retry union
-    //   5. Convert merged geometry back to GeoJSON
-    //
-    // For now, pass through features unchanged (single-geometry-per-code
-    // from gpkg_reader already handles simple cases).
+    if (features.size() <= 1) return features;
 
-    return features;
+    GEOSCtx gctx;
+    auto ctx = gctx.ctx;
+
+    // Group features by code (GID) — same code means same admin area,
+    // multiple geometries need to be unioned into one.
+    std::map<std::string, std::vector<size_t>> groups;
+    for (size_t i = 0; i < features.size(); ++i) {
+        groups[features[i].code].push_back(i);
+    }
+
+    std::vector<boundary::BoundaryFeature> result;
+    result.reserve(groups.size());
+
+    for (auto& [code, indices] : groups) {
+        // Single feature per code — no union needed
+        if (indices.size() == 1) {
+            result.push_back(features[indices[0]]);
+            continue;
+        }
+
+        // Collect GEOS geometries for this group
+        std::vector<GEOSGeometry*> geoms;
+        geoms.reserve(indices.size());
+
+        for (size_t idx : indices) {
+            GEOSGeometry* g = geojson_to_geos(ctx, features[idx].geojson);
+            if (g) geoms.push_back(g);
+        }
+
+        if (geoms.empty()) {
+            result.push_back(features[indices[0]]);
+            continue;
+        }
+
+        if (geoms.size() == 1) {
+            // Only one valid geometry — use it directly
+            auto out = features[indices[0]];
+            out.geojson = geos_to_geojson(ctx, geoms[0]);
+            GEOSGeom_destroy_r(ctx, geoms[0]);
+            result.push_back(std::move(out));
+            continue;
+        }
+
+        // Build GeometryCollection and union
+        GEOSGeometry* collection = GEOSGeom_createCollection_r(
+            ctx, GEOS_GEOMETRYCOLLECTION,
+            geoms.data(), static_cast<unsigned int>(geoms.size())
+        );
+
+        GEOSGeometry* merged = nullptr;
+        if (collection) {
+            merged = GEOSUnaryUnion_r(ctx, collection);
+
+            // Fallback: MakeValid each input and retry
+            if (!merged) {
+                std::cerr << "[geo_merge] Union failed for " << code
+                          << ", attempting MakeValid fallback\n";
+
+                std::vector<GEOSGeometry*> valid_geoms;
+                valid_geoms.reserve(geoms.size());
+                for (auto* g : geoms) {
+                    GEOSGeometry* v = GEOSMakeValid_r(ctx, g);
+                    if (v) valid_geoms.push_back(v);
+                }
+
+                if (!valid_geoms.empty()) {
+                    GEOSGeometry* valid_coll = GEOSGeom_createCollection_r(
+                        ctx, GEOS_GEOMETRYCOLLECTION,
+                        valid_geoms.data(),
+                        static_cast<unsigned int>(valid_geoms.size())
+                    );
+                    if (valid_coll) {
+                        merged = GEOSUnaryUnion_r(ctx, valid_coll);
+                        // createCollection takes ownership, no manual destroy
+                    }
+                }
+            }
+            // createCollection took ownership of geoms, don't destroy individually
+        }
+
+        // Build output feature
+        boundary::BoundaryFeature out;
+        out.code = code;
+        out.name = features[indices[0]].name;
+
+        if (merged) {
+            out.geojson = geos_to_geojson(ctx, merged);
+            GEOSGeom_destroy_r(ctx, merged);
+        } else {
+            // Last resort: collect all polygons into a MultiPolygon manually
+            std::cerr << "[geo_merge] All union attempts failed for " << code
+                      << ", using first geometry as fallback\n";
+            out.geojson = features[indices[0]].geojson;
+        }
+
+        result.push_back(std::move(out));
+    }
+
+    return result;
 }
 
 } // namespace geo_merge

cpp/src/ghs_enrich.cpp

@@ -1,36 +1,337 @@
 #include "ghs_enrich.h"
 
+#include <cmath>
+#include <algorithm>
 #include <stdexcept>
 
-// GDAL + PROJ
+// GDAL C API — handles any size GeoTIFF via windowed reads.
+// The GHS TIFFs are ~3-8 GB global grids (100m Mollweide), but RasterIO
+// only decompresses the bbox-clipped tile range, keeping memory bounded.
 #include <gdal.h>
 #include <gdal_priv.h>
+#include <cpl_conv.h>
+
+// PROJ C API — reusable transform handle, no per-pixel reinit.
 #include <proj.h>
 
-#include "pip.h"
+// OGR for geometry parsing (GeoJSON → OGRGeometry for bbox)
+#include <ogr_api.h>
 
 namespace ghs_enrich {
 
-std::vector<boundary::BoundaryFeature> enrich(
-    const std::vector<boundary::BoundaryFeature>& features,
-    const std::string& /*pop_raster_path*/,
-    const std::string& /*smod_raster_path*/
-) {
-    // TODO: Implement GeoTIFF raster sampling + PROJ transform + PIP loop
-    //
-    // For each feature:
-    //   1. Compute bounding box of geometry
-    //   2. Map bbox to raster pixel window
-    //   3. GDALRasterIO — read pixel block
-    //   4. For each pixel (with stride):
-    //      a. proj_trans: Mollweide → WGS84
-    //      b. point_in_polygon check
-    //      c. Accumulate population / settlement class
-    //   5. Attach enrichment data to BoundaryFeature
-    //
-    // For now, pass through features unchanged.
+// ── Constants ────────────────────────────────────────────────────────────────
+static constexpr float  NODATA_SENTINEL = 65535.f;
+static constexpr double MIN_DIST_METERS = 3000.0;   // 3km separation between peaks
+static constexpr int    MAX_CENTERS     = 5;
 
-    return features;
+// ── Raster handle cache (opened once per process) ────────────────────────────
+struct RasterInfo {
+    GDALDatasetH ds     = nullptr;
+    GDALRasterBandH band = nullptr;
+    double origin_x     = 0;   // top-left X in Mollweide meters
+    double origin_y     = 0;   // top-left Y in Mollweide meters
+    double res_x        = 0;   // pixel width  (positive, ~100m)
+    double res_y        = 0;   // pixel height (negative, ~-100m)
+    int    width        = 0;
+    int    height       = 0;
+};
+
+static RasterInfo open_raster(const std::string& path) {
+    GDALAllRegister();
+
+    GDALDatasetH ds = GDALOpen(path.c_str(), GA_ReadOnly);
+    if (!ds) {
+        throw std::runtime_error("Cannot open raster: " + path);
+    }
+
+    double gt[6];
+    GDALGetGeoTransform(ds, gt);
+
+    RasterInfo ri;
+    ri.ds       = ds;
+    ri.band     = GDALGetRasterBand(ds, 1);
+    ri.origin_x = gt[0];           // top-left X
+    ri.origin_y = gt[3];           // top-left Y
+    ri.res_x    = gt[1];           // pixel width  (+100)
+    ri.res_y    = gt[5];           // pixel height (-100)
+    ri.width    = GDALGetRasterBandXSize(ri.band);
+    ri.height   = GDALGetRasterBandYSize(ri.band);
+    return ri;
+}
+
+// ── Compute stats for one geometry against one raster ────────────────────────
+struct RasterStats {
+    double total   = 0;
+    double maxVal  = 0;
+    double centerLon = 0;
+    double centerLat = 0;
+    std::vector<std::array<double, 3>> centers;
+    bool valid = false;
+};
+
+static RasterStats compute_stats(
+    const RasterInfo& ri,
+    PJ* to_moll,     // EPSG:4326 → EPSG:54009
+    PJ* to_wgs84,    // EPSG:54009 → EPSG:4326
+    double minLon, double minLat,
+    double maxLon, double maxLat,
+    const nlohmann::json& geojson_geometry,
+    bool is_pop
+) {
+    RasterStats stats;
+
+    // ── Project WGS84 bbox → Mollweide ──────────────────────────────────────
+    PJ_COORD sw = proj_coord(minLon, minLat, 0, 0);
+    PJ_COORD ne = proj_coord(maxLon, maxLat, 0, 0);
+    PJ_COORD sw_m = proj_trans(to_moll, PJ_FWD, sw);
+    PJ_COORD ne_m = proj_trans(to_moll, PJ_FWD, ne);
+
+    double mollMinX = std::min(sw_m.xy.x, ne_m.xy.x);
+    double mollMaxX = std::max(sw_m.xy.x, ne_m.xy.x);
+    double mollMinY = std::min(sw_m.xy.y, ne_m.xy.y);
+    double mollMaxY = std::max(sw_m.xy.y, ne_m.xy.y);
+
+    // ── Map to pixel coordinates ────────────────────────────────────────────
+    // geo → pixel: px = (geo_x - origin_x) / res_x
+    //              py = (geo_y - origin_y) / res_y   (res_y is negative)
+    int py1 = static_cast<int>(std::floor((mollMinY - ri.origin_y) / ri.res_y));
+    int py2 = static_cast<int>(std::floor((mollMaxY - ri.origin_y) / ri.res_y));
+    int minPx = std::max(0, static_cast<int>(std::floor((mollMinX - ri.origin_x) / ri.res_x)));
+    int maxPx = std::min(ri.width, static_cast<int>(std::ceil((mollMaxX - ri.origin_x) / ri.res_x)));
+    int minPy = std::max(0, std::min(py1, py2));
+    int maxPy = std::min(ri.height, std::max(py1, py2));
+
+    int widthPx  = maxPx - minPx;
+    int heightPx = maxPy - minPy;
+    if (widthPx <= 0 || heightPx <= 0) return stats;
+
+    // ── Read raster window ──────────────────────────────────────────────────
+    // GDALRasterIO handles tiled/compressed TIFFs natively.
+    // Only the tiles overlapping this window are decompressed — a 6 GB
+    // global TIF won't be loaded into memory.
+    std::vector<float> buf(static_cast<size_t>(widthPx) * heightPx);
+    CPLErr err = GDALRasterIO(
+        ri.band, GF_Read,
+        minPx, minPy, widthPx, heightPx,
+        buf.data(), widthPx, heightPx,
+        GDT_Float32, 0, 0
+    );
+    if (err != CE_None) return stats;
+
+    // ── Pixel loop — accumulate stats ───────────────────────────────────────
+    double totalVal    = 0;
+    double weightedX   = 0;
+    double weightedY   = 0;
+    double maxValFound = 0;
+
+    // Candidate arrays for peak detection
+    struct Candidate {
+        double lon, lat;     // WGS84
+        double moll_x, moll_y; // Mollweide (for distance calc)
+        double val;
+    };
+    std::vector<Candidate> candidates;
+
+    // OGR geometry for PIP — faster than manual ring extraction for
+    // MultiPolygons with holes
+    OGRGeometryH ogr_geom = nullptr;
+    {
+        std::string geojson_str = geojson_geometry.dump();
+        ogr_geom = OGR_G_CreateGeometryFromJson(geojson_str.c_str());
+    }
+
+    for (int y = 0; y < heightPx; ++y) {
+        for (int x = 0; x < widthPx; ++x) {
+            float val = buf[y * widthPx + x];
+
+            if (val <= 0 || val == NODATA_SENTINEL || std::isnan(val)) continue;
+
+            // Pixel center → Mollweide coordinates
+            double pxX = ri.origin_x + (minPx + x + 0.5) * ri.res_x;
+            double pxY = ri.origin_y + (minPy + y + 0.5) * ri.res_y;
+
+            // → WGS84 for PIP check
+            PJ_COORD in = proj_coord(pxX, pxY, 0, 0);
+            PJ_COORD out = proj_trans(to_wgs84, PJ_FWD, in);
+            double ptLon = out.xy.x;
+            double ptLat = out.xy.y;
+
+            // PIP via OGR (handles MultiPolygon + holes natively)
+            bool inside = false;
+            if (ogr_geom) {
+                OGRGeometryH pt = OGR_G_CreateGeometry(wkbPoint);
+                OGR_G_SetPoint_2D(pt, 0, ptLon, ptLat);
+                inside = OGR_G_Contains(ogr_geom, pt) || OGR_G_Intersects(ogr_geom, pt);
+                OGR_G_DestroyGeometry(pt);
+            }
+
+            if (!inside) continue;
+
+            totalVal   += val;
+            weightedX  += (minPx + x) * static_cast<double>(val);
+            weightedY  += (minPy + y) * static_cast<double>(val);
+            if (val > maxValFound) maxValFound = val;
+
+            candidates.push_back({ptLon, ptLat, pxX, pxY, static_cast<double>(val)});
+        }
+    }
+
+    if (ogr_geom) OGR_G_DestroyGeometry(ogr_geom);
+
+    if (totalVal <= 0) return stats;
+
+    // ── Weighted center ─────────────────────────────────────────────────────
+    double centerPx = weightedX / totalVal;
+    double centerPy = weightedY / totalVal;
+    double centerMollX = ri.origin_x + centerPx * ri.res_x;
+    double centerMollY = ri.origin_y + centerPy * ri.res_y;
+    PJ_COORD cc_in = proj_coord(centerMollX, centerMollY, 0, 0);
+    PJ_COORD cc_out = proj_trans(to_wgs84, PJ_FWD, cc_in);
+
+    stats.total     = is_pop ? std::round(totalVal) : totalVal;
+    stats.maxVal    = is_pop ? std::round(maxValFound) : maxValFound;
+    stats.centerLon = cc_out.xy.x;
+    stats.centerLat = cc_out.xy.y;
+    stats.valid     = true;
+
+    // ── Top-N peak detection with 3km separation ────────────────────────────
+    std::sort(candidates.begin(), candidates.end(),
+              [](const Candidate& a, const Candidate& b) { return a.val > b.val; });
+
+    std::vector<size_t> selected;
+    for (size_t i = 0; i < candidates.size() && static_cast<int>(stats.centers.size()) < MAX_CENTERS; ++i) {
+        bool distinct = true;
+        for (size_t si : selected) {
+            double dx = candidates[i].moll_x - candidates[si].moll_x;
+            double dy = candidates[i].moll_y - candidates[si].moll_y;
+            if (std::hypot(dx, dy) < MIN_DIST_METERS) {
+                distinct = false;
+                break;
+            }
+        }
+        if (distinct) {
+            selected.push_back(i);
+            double v = is_pop ? std::round(candidates[i].val) : candidates[i].val;
+            stats.centers.push_back({candidates[i].lon, candidates[i].lat, v});
+        }
+    }
+
+    return stats;
+}
+
+// ── Extract WGS84 bbox from GeoJSON geometry ────────────────────────────────
+static void geojson_bbox(
+    const nlohmann::json& geom,
+    double& minLon, double& minLat,
+    double& maxLon, double& maxLat
+) {
+    // Use OGR to compute envelope
+    std::string str = geom.dump();
+    OGRGeometryH g = OGR_G_CreateGeometryFromJson(str.c_str());
+    if (g) {
+        OGREnvelope env;
+        OGR_G_GetEnvelope(g, &env);
+        minLon = env.MinX;
+        maxLon = env.MaxX;
+        minLat = env.MinY;
+        maxLat = env.MaxY;
+        OGR_G_DestroyGeometry(g);
+    } else {
+        minLon = minLat = maxLon = maxLat = 0;
+    }
+}
+
+// ── Public API ──────────────────────────────────────────────────────────────
+
+EnrichResult enrich_feature(
+    const nlohmann::json& geojson_geometry,
+    const std::string& pop_tiff_path,
+    const std::string& built_tiff_path
+) {
+    EnrichResult result;
+
+    // ── PROJ transforms (created once, reused per-pixel) ────────────────────
+    PJ_CONTEXT* ctx = proj_context_create();
+    PJ* to_moll  = proj_create_crs_to_crs(ctx, "EPSG:4326", "EPSG:54009", nullptr);
+    PJ* to_wgs84 = proj_create_crs_to_crs(ctx, "EPSG:54009", "EPSG:4326", nullptr);
+
+    if (!to_moll || !to_wgs84) {
+        if (to_moll)  proj_destroy(to_moll);
+        if (to_wgs84) proj_destroy(to_wgs84);
+        proj_context_destroy(ctx);
+        throw std::runtime_error("Failed to create PROJ transforms");
+    }
+
+    // ── Bounding box ────────────────────────────────────────────────────────
+    double minLon, minLat, maxLon, maxLat;
+    geojson_bbox(geojson_geometry, minLon, minLat, maxLon, maxLat);
+
+    // ── Population raster ───────────────────────────────────────────────────
+    if (!pop_tiff_path.empty()) {
+        static RasterInfo pop_ri;
+        static bool pop_opened = false;
+        if (!pop_opened) {
+            pop_ri = open_raster(pop_tiff_path);
+            pop_opened = true;
+        }
+
+        auto pop_stats = compute_stats(
+            pop_ri, to_moll, to_wgs84,
+            minLon, minLat, maxLon, maxLat,
+            geojson_geometry, true
+        );
+        if (pop_stats.valid) {
+            result.hasPop           = true;
+            result.ghsPopulation    = pop_stats.total;
+            result.ghsPopMaxDensity = pop_stats.maxVal;
+            result.ghsPopCenterLon  = pop_stats.centerLon;
+            result.ghsPopCenterLat  = pop_stats.centerLat;
+            result.ghsPopCenters    = pop_stats.centers;
+        }
+    }
+
+    // ── Built-up surface raster ─────────────────────────────────────────────
+    if (!built_tiff_path.empty()) {
+        static RasterInfo built_ri;
+        static bool built_opened = false;
+        if (!built_opened) {
+            built_ri = open_raster(built_tiff_path);
+            built_opened = true;
+        }
+
+        auto built_stats = compute_stats(
+            built_ri, to_moll, to_wgs84,
+            minLon, minLat, maxLon, maxLat,
+            geojson_geometry, false
+        );
+        if (built_stats.valid) {
+            result.hasBuilt          = true;
+            result.ghsBuiltWeight    = built_stats.total;
+            result.ghsBuiltMax       = built_stats.maxVal;
+            result.ghsBuiltCenterLon = built_stats.centerLon;
+            result.ghsBuiltCenterLat = built_stats.centerLat;
+            result.ghsBuiltCenters   = built_stats.centers;
+        }
+    }
+
+    proj_destroy(to_moll);
+    proj_destroy(to_wgs84);
+    proj_context_destroy(ctx);
+
+    return result;
+}
+
+std::vector<EnrichResult> enrich_batch(
+    const std::vector<nlohmann::json>& geometries,
+    const std::string& pop_tiff_path,
+    const std::string& built_tiff_path
+) {
+    std::vector<EnrichResult> results;
+    results.reserve(geometries.size());
+    for (const auto& geom : geometries) {
+        results.push_back(enrich_feature(geom, pop_tiff_path, built_tiff_path));
+    }
+    return results;
 }
 
 } // namespace ghs_enrich

cpp/src/ghs_enrich.h

@@ -1,17 +1,51 @@
 #pragma once
 
+#include <string>
 #include <vector>
-#include "types.h"
+#include <nlohmann/json.hpp>
 
 namespace ghs_enrich {
 
-/// Enrich boundary features with GHS-POP/GHS-SMOD raster data.
-/// Samples GeoTIFF rasters using GDAL RasterIO, projects pixels with PROJ,
-/// and uses point-in-polygon to attribute population/settlement data.
-std::vector<boundary::BoundaryFeature> enrich(
-    const std::vector<boundary::BoundaryFeature>& features,
-    const std::string& pop_raster_path = "",
-    const std::string& smod_raster_path = ""
+/// Enrichment results for a single boundary feature.
+struct EnrichResult {
+    // Population grid (GHS-POP)
+    double ghsPopulation   = 0;
+    double ghsPopMaxDensity = 0;
+    double ghsPopCenterLon = 0;
+    double ghsPopCenterLat = 0;
+    std::vector<std::array<double, 3>> ghsPopCenters;  // [lon, lat, density]
+
+    // Built-up surface (GHS-BUILT-S)
+    double ghsBuiltWeight  = 0;
+    double ghsBuiltMax     = 0;
+    double ghsBuiltCenterLon = 0;
+    double ghsBuiltCenterLat = 0;
+    std::vector<std::array<double, 3>> ghsBuiltCenters; // [lon, lat, built]
+
+    bool hasPop   = false;
+    bool hasBuilt = false;
+};
+
+/// Enrich a GeoJSON feature geometry with GHS raster statistics.
+///
+/// The TIFFs are global grids in Mollweide (EPSG:54009) at 100m resolution:
+///   - GHS_POP_E2030_GLOBE_R2023A_54009_100_V1_0.tif  (~3–8 GB)
+///   - GHS_BUILT_S_E2030_GLOBE_R2023A_54009_100_V1_0.tif
+///
+/// GDAL reads these via windowed RasterIO — only the bbox-clipped
+/// portion is loaded, so memory use stays bounded regardless of file size.
+EnrichResult enrich_feature(
+    const nlohmann::json& geojson_geometry,
+    const std::string& pop_tiff_path,
+    const std::string& built_tiff_path
+);
+
+/// Convenience: enrich all features in a batch, returning enrichment data
+/// indexed by feature position.
+std::vector<EnrichResult> enrich_batch(
+    const std::vector<nlohmann::json>& geometries,
+    const std::string& pop_tiff_path,
+    const std::string& built_tiff_path
 );
 
 } // namespace ghs_enrich

cpp/src/main.cpp

@@ -42,7 +42,7 @@ int main(int argc, char* argv[]) {
     app.add_option("--cache-dir", cache_dir,
         "Output directory for boundary JSON files")->default_val("cache/gadm");
     app.add_option("--gpkg", gpkg_path,
-        "Path to GADM GeoPackage file")->default_val("data/gadm_410-levels.gpkg");
+        "Path to GADM GeoPackage file")->default_val("data/gadm_410.gpkg");
     app.add_option("--continent-json", continent_json,
         "Path to gadm_continent.json")->default_val("data/gadm_continent.json");
     app.add_flag("--force", force,
@@ -96,7 +96,7 @@ int main(int argc, char* argv[]) {
     #ifdef HAS_OPENMP
     #pragma omp parallel for schedule(dynamic) reduction(+:processed,skipped,errors)
     #endif
-    for (size_t i = 0; i < countries.size(); ++i) {
+    for (int i = 0; i < static_cast<int>(countries.size()); ++i) {
         const auto& code = countries[i];
 
         for (int lvl : levels) {

docs/cpp-port.md

@@ -203,8 +203,8 @@ Output goes to the same `cache/gadm/` directory — the Node.js server reads the
 |---|---|---|
 | `main.cpp` | ✅ Done | CLI harness, cache logic, JSON output, OpenMP |
 | `gpkg_reader` | ✅ Scaffolded | OGR read + feature grouping — needs GEOS geometry collection |
-| `geo_merge` | 🔲 Stub | Passthrough — implement `GEOSUnaryUnion` + `MakeValid` fallback |
-| `ghs_enrich` | 🔲 Stub | Passthrough — implement `GDALRasterIO` + PROJ + PIP loop |
+| `geo_merge` | ✅ Done | GEOSUnaryUnion + MakeValid fallback, WKB lossless roundtrip |
+| `ghs_enrich` | ✅ Done | GDAL windowed RasterIO + PROJ + OGR PIP + peak detection |
 | `pip.h` | ✅ Done | Inline ray-casting, unit-tested |
 | `types.h` | ✅ Done | `BoundaryFeature`, `BoundaryResult` structs |
 

package.json

@@ -29,6 +29,7 @@
         "test:ghs-pop": "vitest run src/__tests__/ghs-pop.test.ts",
         "test:enrich": "vitest run src/__tests__/enrich-ghs.test.ts",
         "boundaries": "npx tsx scripts/boundaries.ts",
+        "boundaries:cpp": ".\\cpp\\build\\Release\\boundaries.exe",
         "refresh": "npx tsx scripts/refresh-database.ts"
     },
     "devDependencies": {