Repository: ungarj/mapchete
Pull request #461 (2022-04-28 11:04:43)

add cli command to write STACTA prototype datasets required by GDAL


Patch:
diff --git a/mapchete/cli/default/stac.py b/mapchete/cli/default/stac.py
index 402ef5c..00e8a83 100644
--- a/mapchete/cli/default/stac.py
+++ b/mapchete/cli/default/stac.py
@@ -5,10 +5,11 @@ import json
 import os
 import oyaml as yaml
 
+import mapchete
 from mapchete.cli import options
 from mapchete.config import raw_conf, raw_conf_output_pyramid
 from mapchete.formats import read_output_metadata
-from mapchete.stac import tile_directory_stac_item
+from mapchete.stac import tile_directory_stac_item, create_prototype_files
 from mapchete.validate import validate_zooms
 
 
@@ -128,3 +129,16 @@ def output_info(inp):
         None,
         None,
     )
+
+
+@stac.command(name="create-prototype-files", help="Create STAC item prototype files.")
+@options.arg_input
+@options.opt_force
+@options.opt_debug
+def prototype_files(
+    input_,
+    force=None,
+    **kwargs,
+):
+    with mapchete.open(input_, mode="readonly") as mp:
+        create_prototype_files(mp)
diff --git a/mapchete/io/raster.py b/mapchete/io/raster.py
index 5ab37bc..ca0f9c2 100644
--- a/mapchete/io/raster.py
+++ b/mapchete/io/raster.py
@@ -458,6 +458,7 @@ def write_raster_window(
     out_path=None,
     tags=None,
     fs=None,
+    write_empty=False,
     **kwargs,
 ):
     """
@@ -501,7 +502,7 @@ def write_raster_window(
         out_profile["transform"] = out_profile.pop("affine")
 
     # write if there is any band with non-masked data
-    if window_data.all() is not ma.masked:
+    if write_empty or (window_data.all() is not ma.masked):
 
         try:
             with rasterio_write(out_path, "w", fs=fs, **out_profile) as dst:
diff --git a/mapchete/stac.py b/mapchete/stac.py
index 462f49b..0242f81 100644
--- a/mapchete/stac.py
+++ b/mapchete/stac.py
@@ -1,16 +1,17 @@
 from collections import OrderedDict
 import datetime
 import logging
-from nis import match
+import numpy as np
+import numpy.ma as ma
 import os
 from pyproj import CRS
-from typing import Type
 from shapely.geometry import box, mapping
 
 from tilematrix._funcs import Bounds
 
 from mapchete.errors import ReprojectionFailed
-from mapchete.io import fs_from_path
+from mapchete.io import makedirs
+from mapchete.io.raster import write_raster_window
 from mapchete.io.vector import reproject_geometry
 from mapchete.tile import BufferedTilePyramid
 
@@ -439,3 +440,34 @@ def zoom_levels_from_item(item):
         return zoom_levels
     else:
         raise AttributeError("STAC item does have tile matrix sets defined")
+
+
+def create_prototype_files(mp):
+    # for each zoom level get tile output for 0/0
+    for zoom in mp.config.init_zoom_levels:
+        prototype_tile = mp.config.output_pyramid.tile(zoom, 0, 0)
+        tile_path = mp.config.output.get_path(prototype_tile)
+        # if tile exists, skip
+        if mp.config.output.tiles_exist(output_tile=prototype_tile):
+            logger.debug("prototype tile %s already exists", tile_path)
+            pass
+        # if not, write empty tile
+        else:
+            logger.debug("creating prototype tile %s", tile_path)
+            out_profile = mp.config.output.profile(prototype_tile)
+            makedirs(os.path.dirname(tile_path))
+            write_raster_window(
+                in_tile=prototype_tile,
+                in_data=ma.masked_array(
+                    data=np.full(
+                        (out_profile["count"],) + prototype_tile.shape,
+                        out_profile["nodata"],
+                        dtype=out_profile["dtype"],
+                    ),
+                    mask=True,
+                ),
+                out_profile=out_profile,
+                out_tile=prototype_tile,
+                out_path=tile_path,
+                write_empty=True,
+            )


Test patch:
diff --git a/test/test_cli.py b/test/test_cli.py
index c1978bf..8c686f5 100644
--- a/test/test_cli.py
+++ b/test/test_cli.py
@@ -4,10 +4,10 @@ from click.testing import CliRunner
 import fiona
 import geobuf
 import os
+from packaging import version
 import pytest
 from shapely import wkt
 from shapely.geometry import shape
-import shapely.geos
 import rasterio
 from rasterio.io import MemoryFile
 from rio_cogeo.cogeo import cog_validate
@@ -17,7 +17,6 @@ import yaml
 import mapchete
 from mapchete.cli.main import main as mapchete_cli
 from mapchete.cli import options
-from mapchete.errors import MapcheteProcessOutputError
 
 
 SCRIPTDIR = os.path.dirname(os.path.realpath(__file__))
@@ -1499,6 +1498,7 @@ def test_fs_opt_extractor():
 
 
 def test_stac_mapchete_file(cleantopo_br):
+    run_cli(["execute", cleantopo_br.path])
     run_cli(["stac", "create-item", cleantopo_br.path, "-z", "5", "--force"])
 
 
@@ -1515,3 +1515,13 @@ def test_stac_tiledir(http_tiledir, mp_tmpdir):
             f"{mp_tmpdir}/stac_example.json",
         ]
     )
+
+
+@pytest.mark.skipif(
+    version.parse(rasterio.__gdal_version__) < version.parse("3.3.0"),
+    reason="required STACTA driver is only available in GDAL>=3.3.0",
+)
+def test_stac_prototype_files(cleantopo_br):
+    run_cli(["execute", cleantopo_br.path])
+    run_cli(["stac", "create-prototype-files", cleantopo_br.path])
+    rasterio.open(cleantopo_br.mp().config.output.stac_path)
diff --git a/test/test_stac.py b/test/test_stac.py
index 6316613..77c8179 100644
--- a/test/test_stac.py
+++ b/test/test_stac.py
@@ -1,13 +1,19 @@
 import json
+from packaging import version
 import pytest
+import rasterio
+from rasterio.errors import RasterioIOError
 from shapely.geometry import box, shape
 
+from mapchete.commands import execute
+from mapchete.io import fs_from_path
 from mapchete.io.vector import reproject_geometry
 from mapchete.stac import (
     tile_directory_stac_item,
     update_tile_directory_stac_item,
     tile_pyramid_from_item,
     zoom_levels_from_item,
+    create_prototype_files,
 )
 from mapchete.tile import BufferedTilePyramid
 
@@ -244,3 +250,40 @@ def test_zoom_levels_from_item_errors():
     item.properties = {}
     with pytest.raises(AttributeError):
         zoom_levels_from_item(item)
+
+
+@pytest.mark.skipif(
+    version.parse(rasterio.__gdal_version__) < version.parse("3.3.0"),
+    reason="required STACTA driver is only available in GDAL>=3.3.0",
+)
+def test_create_prototype_file(example_mapchete):
+    # create sparse tiledirectory with no tiles at row/col 0/0
+    execute(example_mapchete.dict, zoom=[10, 11])
+
+    # read STACTA with rasterio and expect an exception
+    stac_path = example_mapchete.mp().config.output.stac_path
+    assert fs_from_path(stac_path).exists(stac_path)
+
+    with pytest.raises(RasterioIOError):
+        rasterio.open(stac_path)
+
+    # create prototype file and assert reading is possible
+    create_prototype_files(example_mapchete.mp())
+    rasterio.open(stac_path)
+
+
+@pytest.mark.skipif(
+    version.parse(rasterio.__gdal_version__) < version.parse("3.3.0"),
+    reason="required STACTA driver is only available in GDAL>=3.3.0",
+)
+def test_create_prototype_file_exists(cleantopo_tl):
+    # create sparse tiledirectory with no tiles at row/col 0/0
+    execute(cleantopo_tl.dict)
+
+    # read STACTA with rasterio and expect an exception
+    stac_path = cleantopo_tl.mp().config.output.stac_path
+    assert fs_from_path(stac_path).exists(stac_path)
+
+    # create prototype file and assert reading is possible
+    create_prototype_files(cleantopo_tl.mp())
+    rasterio.open(stac_path)