Fix overviews

CHANGELOG.md

@@ -25,6 +25,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ### Fixed
 
+- Colormap must be applied before `gdal_translate` is called. This ensures that the correct settings are used e.g. NEAREST vs CUBIC for generating overviews.[#22](https://github.com/stactools-packages/nrcan-landcover/pull/22)
 - create-item and create-collection shouldn't call save on the STAC object.
 - Metadata URL on STAC Item.
 

src/stactools/nrcan_landcover/cog.py

@@ -1,6 +1,7 @@
 import logging
 import os
 from glob import glob
+from shutil import copyfile
 from subprocess import CalledProcessError, check_output
 from tempfile import TemporaryDirectory
 from zipfile import ZipFile
@@ -11,6 +12,7 @@
 from stactools.nrcan_landcover.constants import (
     COLOUR_MAP,
     JSONLD_HREF,
+    NO_DATA_VALUE,
     TILING_PIXEL_SIZE,
 )
 from stactools.nrcan_landcover.utils import get_metadata
@@ -149,50 +151,59 @@ def create_cog(
         str: The path to the output COG.
     """
 
-    output = None
-    try:
-        if dry_run:
-            logger.info("Would have read TIFF, created COG, and written COG")
-        else:
-            logger.info("Converting TIFF to COG")
-            logger.debug(f"input_path: {input_path}")
-            logger.debug(f"output_path: {output_path}")
-            cmd = [
-                "gdal_translate",
-                "-of",
-                "COG",
-                "-co",
-                "NUM_THREADS=ALL_CPUS",
-                "-co",
-                "BLOCKSIZE=512",
-                "-co",
-                "COMPRESS=DEFLATE",
-                "-co",
-                "LEVEL=9",
-                "-co",
-                "PREDICTOR=YES",
-                "-co",
-                "OVERVIEWS=IGNORE_EXISTING",
-                "-a_nodata",
-                "0",
-                input_path,
-                output_path,
-            ]
-
-            try:
-                output = check_output(cmd)
-            except CalledProcessError as e:
-                output = e.output
-                raise
-            finally:
-                logger.info(f"output: {str(output)}")
-            with rasterio.open(output_path, "r+") as dataset:
-                dataset.write_colormap(1, COLOUR_MAP)
+    with TemporaryDirectory() as tmp_dir:
+        # The colormap must be applied before processing.
+        # This ensures that the correct defaults are used.
+        # e.g. NEAREST rather than CUBIC for calculating overviews.
+        # So we copy the input file to a temp dir and apply it before
+        # running gdal_translate.
+        tmp_path = os.path.join(tmp_dir, os.path.basename(input_path))
+        copyfile(input_path, tmp_path)
+        with rasterio.open(tmp_path, "r+") as dataset:
+            dataset.write_colormap(1, COLOUR_MAP)
+        output = None
+        try:
+            if dry_run:
+                logger.info(
+                    "Would have read TIFF, created COG, and written COG")
+            else:
+                logger.info("Converting TIFF to COG")
+                logger.debug(f"input_path: {input_path}")
+                logger.debug(f"output_path: {output_path}")
+                cmd = [
+                    "gdal_translate",
+                    "-of",
+                    "COG",
+                    "-co",
+                    "NUM_THREADS=ALL_CPUS",
+                    "-co",
+                    "BLOCKSIZE=512",
+                    "-co",
+                    "COMPRESS=DEFLATE",
+                    "-co",
+                    "LEVEL=9",
+                    "-co",
+                    "PREDICTOR=YES",
+                    "-co",
+                    "OVERVIEWS=IGNORE_EXISTING",
+                    "-a_nodata",
+                    str(NO_DATA_VALUE),
+                    tmp_path,
+                    output_path,
+                ]
 
-    except Exception:
-        logger.error("Failed to process {}".format(output_path))
+                try:
+                    output = check_output(cmd)
+                except CalledProcessError as e:
+                    output = e.output
+                    raise
+                finally:
+                    logger.info(f"output: {str(output)}")
 
-        if raise_on_fail:
-            raise
+        except Exception:
+            logger.error("Failed to process {}".format(output_path))
+
+            if raise_on_fail:
+                raise
 
     return output_path

src/stactools/nrcan_landcover/constants.py

@@ -89,6 +89,8 @@
     19: "Snow and Ice",
 }
 
+NO_DATA_VALUE = 0
+
 TILING_PIXEL_SIZE = (10001, 10001)
 
 FULL_DATASET_BBOX = [-2600030.0, -885090.0, 3100000.0, 3914940.0]

src/stactools/nrcan_landcover/stac.py

@@ -43,6 +43,7 @@
     LANDCOVER_TITLE,
     LICENSE,
     LICENSE_LINK,
+    NO_DATA_VALUE,
     NRCAN_PROVIDER,
     THUMBNAIL_HREF,
 )
@@ -247,7 +248,7 @@ def create_item(metadata: Dict[str, Any],
         # Raster Extension
         cog_asset_raster = RasterExtension.ext(cog_asset, add_if_missing=True)
         cog_asset_raster.bands = [
-            RasterBand.create(nodata=0,
+            RasterBand.create(nodata=NO_DATA_VALUE,
                               sampling=Sampling.AREA,
                               data_type=DataType.UINT8,
                               spatial_resolution=30)
@@ -378,7 +379,7 @@ def create_collection(
             "title":
             "Land cover of Canada COG",
             "raster:bands": [
-                RasterBand.create(nodata=0,
+                RasterBand.create(nodata=NO_DATA_VALUE,
                                   sampling=Sampling.AREA,
                                   data_type=DataType.UINT8,
                                   spatial_resolution=30).to_dict()