Merge pull request #207 from Vizzuality/data/nigeria_fcs

SofiaAldabet · web-flow · commit 1c2cc5d87c2c · 2026-04-16T16:23:13.000+02:00
Nigeria FCS — Sentinel-1 MTC layer processing
diff --git a/data-processing/notebooks/01_raster_layers.ipynb b/data-processing/notebooks/01_raster_layers.ipynb
@@ -35,7 +35,7 @@
     "\n",
     "sys.path.append(\"../src\")\n",
     "\n",
-    "from data_processing.process_layers import process_raster_layers\n",
+    "from data_processing.process_layers import process_prestyled_raster, process_raster_layers\n",
     "from data_processing.utils import clip_rasters_by_vector, merge_tifs, resample_raster\n",
     "\n",
     "project_root = Path(\"/Users/sofia/Documents/Repos/esa/data-processing\")\n",
@@ -1036,6 +1036,29 @@
     "    upload_override=True\n",
     ")"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "nigeria-fcs-mtc-header",
+   "metadata": {},
+   "source": [
+    "#### Nigeria - FCS"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "nigeria-fcs-mtc-cell",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "process_prestyled_raster(\n",
+    "    input_file=project_root / \"data/raw/FCS/Nigeria/Rasters/MTC/SENTINEL-1_MTC_IW1-IW2-IW3_VV_RIGHT_ASCENDING_20150523174456_20150604174457.tif\",\n",
+    "    output_file=project_root / \"data/processed/FCS/Nigeria/Rasters/MTC/sentinel1_mtc.tif\",\n",
+    "    layer_name=\"Nigeria_Sentinel1_MTC\",\n",
+    "    upload=False,\n",
+    ")"
+   ]
   }
  ],
  "metadata": {
@@ -1059,4 +1082,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 0
-}
+}
diff --git a/data-processing/src/data_processing/process_layers.py b/data-processing/src/data_processing/process_layers.py
@@ -8,8 +8,13 @@
 
 import yaml
 
+import os
+
 from .raster_processor import RasterProcessor
 from .vector_processor import VectorProcessor
+from helpers.cog_converter import COGConverter
+from helpers.mapbox_uploader import upload_to_mapbox
+from helpers.mbtiles_converter import MBTilesConverterFactory
 
 # ── Shared helpers ────────────────────────────────────────────────────────────
 
@@ -32,6 +37,140 @@ def _resolve_upload_flag(
     return layer_config.get('upload', False)
 
 
+# ── Pre-styled raster (no QML) ───────────────────────────────────────────────
+
+def _native_max_zoom(file: Path) -> int:
+    """Return the native max zoom of a raster based on its resolution."""
+    import math
+    import rasterio
+    from rasterio.warp import transform_bounds
+    with rasterio.open(file) as src:
+        bounds = transform_bounds(src.crs, "EPSG:4326", *src.bounds)
+        res_x = abs(src.transform.a)
+        # Convert native resolution to degrees, then to web mercator metres
+        lng_span = bounds[2] - bounds[0]
+        px_deg = lng_span / src.width if src.crs.is_geographic else None
+        if px_deg is None:
+            # Projected: convert pixel size to approximate degrees at scene centre
+            centre_lat = (bounds[1] + bounds[3]) / 2
+            metres_per_deg = 111320 * math.cos(math.radians(centre_lat))
+            px_deg = res_x / metres_per_deg
+    # zoom = log2(360 / (px_deg * 256))
+    zoom = math.log2(360 / (px_deg * 256))
+    return max(1, round(zoom))
+
+
+def _normalize_to_uint8(input_file: Path, output_file: Path, percentile: float = 2.0) -> None:
+    """
+    Stretch each band to uint8 using percentile clipping.
+
+    Values below the low percentile → 0, above the high percentile → 255.
+    This avoids black output when the source is float32 with a narrow value range.
+
+    Args:
+        input_file: Source float raster.
+        output_file: Destination uint8 GeoTIFF.
+        percentile: Low/high percentile used for clipping (default 2 / 98).
+    """
+    import numpy as np
+    import rasterio
+
+    with rasterio.open(input_file) as src:
+        meta = src.meta.copy()
+        data = src.read()
+        nodata = src.nodata
+
+    bands, height, width = data.shape
+    out = np.zeros((bands, height, width), dtype=np.uint8)
+
+    for i in range(bands):
+        band = data[i].astype(np.float64)
+        valid_mask = band != nodata if nodata is not None else np.ones_like(band, dtype=bool)
+        valid = band[valid_mask]
+        if valid.size == 0:
+            continue
+        lo = np.percentile(valid, percentile)
+        hi = np.percentile(valid, 100 - percentile)
+        if hi == lo:
+            continue
+        stretched = np.clip((band - lo) / (hi - lo) * 255, 0, 255)
+        out[i] = stretched.astype(np.uint8)
+
+    meta.update({"dtype": "uint8", "nodata": None})
+    with rasterio.open(output_file, "w", **meta) as dst:
+        dst.write(out)
+
+
+def process_prestyled_raster(
+    input_file: Path,
+    output_file: Path,
+    layer_name: str,
+    max_zoom: int = None,
+    percentile: float = 2.0,
+    upload: bool = False,
+):
+    """
+    Convert an already-styled (RGB/RGBA) raster to COG + MBTiles and optionally
+    upload to Mapbox. Use this when the source file already has colour bands and
+    no QML styling is needed.
+
+    Float32 inputs are automatically normalized to uint8 via percentile stretching
+    before COG/MBTiles conversion (otherwise tiles render as all black).
+
+    Args:
+        input_file: Path to the input GeoTIFF.
+        output_file: Path for the output GeoTIFF / MBTiles (same stem, different suffix).
+        layer_name: Display name used when uploading to Mapbox.
+        max_zoom: Maximum zoom level for COG conversion. If None, auto-detected
+                  from the file's native resolution.
+        percentile: Percentile used for float→uint8 contrast stretch (default 2/98).
+        upload: Whether to upload the MBTiles to Mapbox.
+    """
+    import rasterio
+
+    input_file = Path(input_file)
+    output_file = Path(output_file)
+    output_file.parent.mkdir(parents=True, exist_ok=True)
+
+    if max_zoom is None:
+        max_zoom = _native_max_zoom(input_file)
+        print(f"🔍 Auto-detected max zoom: {max_zoom}")
+
+    # Normalize float data to uint8 so PNG tiles render correctly
+    with rasterio.open(input_file) as src:
+        needs_normalization = src.dtypes[0] not in ("uint8", "byte")
+
+    cog_input = input_file
+    if needs_normalization:
+        normalized_path = output_file.parent / f"normalized_{output_file.name}"
+        print(f"⚖️  Normalizing float32 → uint8 (percentile clip {percentile}/{100-percentile})")
+        _normalize_to_uint8(input_file, normalized_path, percentile=percentile)
+        cog_input = normalized_path
+
+    print(f"☁️  Converting to COG: {cog_input.name}")
+    COGConverter.convert(cog_input, output_file, max_zoom=max_zoom)
+
+    if needs_normalization:
+        normalized_path.unlink()
+
+    tile_format = "JPEG" if needs_normalization else "PNG"
+    mbtiles_path = output_file.with_suffix(".mbtiles")
+    print(f"📦 Converting to MBTiles ({tile_format}): {mbtiles_path.name}")
+    MBTilesConverterFactory.convert(output_file, mbtiles_path, tile_format=tile_format)
+
+    if upload:
+        print(f"☁️  Uploading to Mapbox: {layer_name}")
+        upload_to_mapbox(
+            source=mbtiles_path,
+            display_name=layer_name,
+            username=os.getenv("MAPBOX_USER"),
+            token=os.getenv("MAPBOX_TOKEN"),
+        )
+
+    print(f"✅ Done: {mbtiles_path}")
+    return mbtiles_path
+
+
 # ── Raster ────────────────────────────────────────────────────────────────────
 
 def _process_rasters(
diff --git a/data-processing/src/helpers/mbtiles_converter.py b/data-processing/src/helpers/mbtiles_converter.py
@@ -25,23 +25,36 @@ class GeoTIFFConverter(MBTilesConverter):
     """Converter for GeoTIFF files."""
 
     def convert(
-        self, input_path: Path, output_path: Path, overview_levels: list[int] = None
+        self,
+        input_path: Path,
+        output_path: Path,
+        overview_levels: list[int] = None,
+        tile_format: str = "PNG",
     ) -> None:
-        """Convert a Cloud-Optimized GeoTIFF to MBTiles with overviews."""
+        """Convert a Cloud-Optimized GeoTIFF to MBTiles with overviews.
+
+        Args:
+            tile_format: ``"PNG"`` (lossless, best for classified maps) or
+                         ``"JPEG"`` (lossy, much smaller for continuous imagery).
+        """
         if overview_levels is None:
             overview_levels = [2, 4, 8, 16, 32, 64]
-        self._convert_to_mbtiles(input_path, output_path)
+        self._convert_to_mbtiles(input_path, output_path, tile_format=tile_format)
+        self._add_overviews(output_path, overview_levels)
 
-    def _convert_to_mbtiles(self, input_path: Path, output_path: Path) -> None:
-        # Ensure output directory exists
+    def _convert_to_mbtiles(
+        self, input_path: Path, output_path: Path, tile_format: str = "PNG"
+    ) -> None:
         output_path.parent.mkdir(parents=True, exist_ok=True)
 
         command = [
             "gdal_translate",
             "-of",
             "MBTiles",
+            "-a_nodata",
+            "none",
             "-co",
-            "TILE_FORMAT=PNG",
+            f"TILE_FORMAT={tile_format}",
             "-co",
             "BLOCKSIZE=512",
             "--config",
@@ -142,10 +155,12 @@ def create_converter(cls, file_path: Path) -> MBTilesConverter:
         return converter_class()
 
     @classmethod
-    def convert(cls, input_path: Path, output_path: Path, **kwargs) -> None:
+    def convert(
+        cls, input_path: Path, output_path: Path, tile_format: str = "PNG", **kwargs
+    ) -> None:
         """Convenience method to convert file using appropriate converter."""
         converter = cls.create_converter(input_path)
-        converter.convert(input_path, output_path, **kwargs)
+        converter.convert(input_path, output_path, tile_format=tile_format, **kwargs)
 
     @classmethod
     def register_converter(cls, extension: str, converter_class: Type[MBTilesConverter]) -> None:
