Merge pull request #342 from ungarj/fix_325_cog_overviews

Fix single GTiff overview creation

mapchete/_processing.py

@@ -359,15 +359,16 @@ def __init__(
         self.start_method = start_method
         self.max_workers = max_workers or os.cpu_count()
         self.multiprocessing_module = multiprocessing_module
-        logger.debug(
-            "init %s Executor with start_method %s and %s workers",
-            self.multiprocessing_module,
-            self.start_method,
-            self.max_workers,
-        )
-        self._pool = self.multiprocessing_module.get_context(self.start_method).Pool(
-            self.max_workers
-        )
+        if self.max_workers != 1:
+            logger.debug(
+                "init %s Executor with start_method %s and %s workers",
+                self.multiprocessing_module,
+                self.start_method,
+                self.max_workers,
+            )
+            self._pool = self.multiprocessing_module.get_context(
+                self.start_method
+            ).Pool(self.max_workers)
 
     def as_completed(
         self, func=None, iterable=None, fargs=None, fkwargs=None, chunksize=1
@@ -393,14 +394,16 @@ def as_completed(
 
     def __enter__(self):
         """Enter context manager."""
-        self._pool.__enter__()
+        if self.max_workers != 1:
+            self._pool.__enter__()
         return self
 
     def __exit__(self, *args):
         """Exit context manager."""
-        logger.debug("closing %s and workers", self._pool)
-        self._pool.__exit__(*args)
-        logger.debug("%s closed", self._pool)
+        if self.max_workers != 1:
+            logger.debug("closing %s and workers", self._pool)
+            self._pool.__exit__(*args)
+            logger.debug("%s closed", self._pool)
 
 
 class FinishedTask:

mapchete/formats/default/gtiff.py

@@ -41,6 +41,7 @@
 from rasterio.enums import Resampling
 from rasterio.io import MemoryFile
 from rasterio.profiles import Profile
+from rasterio.rio.overview import get_maximum_overview_level
 from rasterio.shutil import copy
 from rasterio.windows import from_bounds
 from shapely.geometry import box
@@ -336,7 +337,7 @@ def profile(self, tile=None):
                 k: v
                 for k, v in self.output_params.items()
                 if k not in _OUTPUT_PARAMETERS
-            }
+            },
         )
         dst_metadata.pop("transform", None)
         if tile is not None:
@@ -422,16 +423,6 @@ def __init__(self, output_params, **kwargs):
             raise ValueError("single file output only works with one zoom level")
         self.zoom = output_params["delimiters"]["zoom"][0]
         self.cog = output_params.get("cog", False)
-        if self.cog or "overviews" in output_params:
-            self.overviews = True
-            self.overviews_resampling = output_params.get(
-                "overviews_resampling", "nearest"
-            )
-            self.overviews_levels = output_params.get(
-                "overviews_levels", [2 ** i for i in range(1, self.zoom + 1)]
-            )
-        else:
-            self.overviews = False
         self.in_memory = output_params.get("in_memory", True)
         _bucket = self.path.split("/")[2] if self.path.startswith("s3://") else None
         self._bucket_resource = get_boto3_bucket(_bucket) if _bucket else None
@@ -480,16 +471,43 @@ def prepare(self, process_area=None, **kwargs):
                     k: self.output_params.get(k, GTIFF_DEFAULT_PROFILE[k])
                     for k in GTIFF_DEFAULT_PROFILE.keys()
                 },
-                **creation_options
+                **creation_options,
             ),
-            bigtiff=self.output_params.get("bigtiff", "NO")
+            bigtiff=self.output_params.get("bigtiff", "NO"),
         )
         logger.debug("single GTiff profile: %s", self._profile)
+
+        logger.debug(
+            get_maximum_overview_level(
+                width, height, minsize=self._profile["blockxsize"]
+            )
+        )
+        if self.cog or "overviews" in self.output_params:
+            self.overviews = True
+            self.overviews_resampling = self.output_params.get(
+                "overviews_resampling", "nearest"
+            )
+            self.overviews_levels = self.output_params.get(
+                "overviews_levels",
+                [
+                    2 ** i
+                    for i in range(
+                        1,
+                        get_maximum_overview_level(
+                            width, height, minsize=self._profile["blockxsize"]
+                        ),
+                    )
+                ],
+            )
+        else:
+            self.overviews = False
+
         self.in_memory = (
             self.in_memory
             if self.in_memory is False
             else height * width < IN_MEMORY_THRESHOLD
         )
+
         # set up rasterio
         if path_exists(self.path):
             if self.output_params["mode"] != "overwrite":
@@ -506,11 +524,13 @@ def prepare(self, process_area=None, **kwargs):
         # (1) use memfile if output is remote or COG
         if self.cog or path_is_remote(self.path):
             if self.in_memory:
+                logger.debug("create MemoryFile")
                 self._memfile = self._ctx.enter_context(MemoryFile())
                 self.dst = self._ctx.enter_context(self._memfile.open(**self._profile))
             else:
                 # in case output raster is too big, use tempfile on disk
                 self._tempfile = self._ctx.enter_context(NamedTemporaryFile())
+                logger.debug(f"create tempfile in {self._tempfile.name}")
                 self.dst = self._ctx.enter_context(
                     rasterio.open(self._tempfile.name, "w+", **self._profile)
                 )
@@ -601,7 +621,7 @@ def write(self, process_tile, data):
                         *out_tile.bounds,
                         transform=self.dst.transform,
                         height=self.dst.height,
-                        width=self.dst.width
+                        width=self.dst.width,
                     )
                     .round_lengths(pixel_precision=0)
                     .round_offsets(pixel_precision=0)
@@ -646,7 +666,9 @@ def close(self, exc_type=None, exc_value=None, exc_traceback=None):
                         self.overviews_levels, Resampling[self.overviews_resampling]
                     )
                     self.dst.update_tags(
-                        ns="rio_overview", resampling=self.overviews_resampling
+                        OVR_RESAMPLING_ALG=Resampling[
+                            self.overviews_resampling
+                        ].name.upper()
                     )
                 # write
                 if self.cog:
@@ -663,7 +685,7 @@ def close(self, exc_type=None, exc_value=None, exc_traceback=None):
                                 self.dst,
                                 tmp_dst.name,
                                 copy_src_overviews=True,
-                                **self._profile
+                                **self._profile,
                             )
                             self._bucket_resource.upload_file(
                                 Filename=tmp_dst.name,
@@ -676,7 +698,7 @@ def close(self, exc_type=None, exc_value=None, exc_traceback=None):
                             self.dst,
                             self.path,
                             copy_src_overviews=True,
-                            **self._profile
+                            **self._profile,
                         )
                 else:
                     if path_is_remote(self.path):

test/test_cli.py

@@ -400,7 +400,7 @@ def test_convert_single_gtiff_overviews(cleantopo_br_tif, mp_tmpdir):
             "--output-pyramid",
             "geodetic",
             "-z",
-            "3",
+            "7",
             "--overviews",
             "--overviews-resampling-method",
             "bilinear",

test/test_formats_geotiff.py

@@ -314,7 +314,12 @@ def test_output_single_gtiff_overviews(output_single_gtiff):
 
     with rasterio.open(mp.config.output.path) as src:
         assert src.overviews(1)
-        assert src.tags(ns="rio_overview").get("resampling") == "bilinear"
+        assert src.tags().get("OVR_RESAMPLING_ALG").lower() == "bilinear"
+        for o in [1, 2, 4, 8]:
+            a = src.read(
+                masked=True, out_shape=(1, int(src.height / o), int(src.width / o))
+            )
+            assert not a.mask.all()
 
 
 @pytest.mark.remote