diff --git a/src/temporal/t.rast.exporttree/Makefile b/src/temporal/t.rast.exporttree/Makefile new file mode 100644 index 0000000..5758aeb --- /dev/null +++ b/src/temporal/t.rast.exporttree/Makefile @@ -0,0 +1,7 @@ +MODULE_TOPDIR = ../../ + +PGM = t.rast.exporttree + +include $(MODULE_TOPDIR)/include/Make/Script.make + +default: script $(TEST_DST) diff --git a/src/temporal/t.rast.exporttree/t.rast.exporttree.html b/src/temporal/t.rast.exporttree/t.rast.exporttree.html new file mode 100644 index 0000000..5e5f654 --- /dev/null +++ b/src/temporal/t.rast.exporttree/t.rast.exporttree.html @@ -0,0 +1,48 @@ +

DESCRIPTION

+ +t.rast.exporttree is a wrapper module around r.out.gdalto +export raster maps from a STRDS to GeoTiffs in a temporal directory tree. The +tree structure can be defined with the temporal_tree option and the +s-flag (to include the semantic_label in the directory +structure). The output directory needs to exist, but the directory structure +below will be created if necessary. Parallel export of maps is supported +with the nprocs option. + +

EXAMPLES

+ +
+temp_dir=$(g.tempfile -d pid=1)
+mkdir $temp_dir
+target_dir=$(g.tempfile -d pid=1)
+mkdir $target_dir
+
+g.region -ag s=0 n=80 w=0 e=120 res=1
+r.external.out format="GTiff" directory=$temp_dir extension="tif" options="COMPRESS=LZW"
+for rmap_idx in 1 2 3
+do
+  for prefix in a b
+  do
+    r.mapcalc expression="${prefix}_${rmap_idx} = ${rmap_idx}00 --overwrite
+    r.support map="${prefix}_${rmap_idx} semantic_label=$prefix
+  done
+done
+t.create type="strds" temporaltype="absolute" output="A" \
+    title="A test" description="A test" --overwrite
+t.register -i type="raster" input="A" maps="a_1,a_2,a_3" \
+    start="2001-01-01" increment="3 months" --overwrite
+t.create type="strds" temporaltype="absolute" output="B" \
+    title="B test" description="B test" --overwrite
+t.register -i type="raster" input="B" maps="b_1,b_2,b_3" \
+    start="2001-01-01" increment="1 day" --overwrite
+
+t.rast.exporttree -f input="A" temporal_tree="%Y/%m" nprocs=2 \
+    output_directory=$target_dir
+
+t.rast.copytree -s input="B" temporal_tree="%Y/%m/%d" nprocs=2 \
+output_directory=$target_dir
+
+ + +

AUTHOR

+ +Stefan Blumentrath diff --git a/src/temporal/t.rast.exporttree/t.rast.exporttree.py b/src/temporal/t.rast.exporttree/t.rast.exporttree.py new file mode 100644 index 0000000..24cd7f6 --- /dev/null +++ b/src/temporal/t.rast.exporttree/t.rast.exporttree.py @@ -0,0 +1,322 @@ +#! /usr/bin/python3 +"""MODULE: t.rast.exporttree +AUTHOR(S): Stefan Blumentrath +PURPOSE: Transfer raster map files from STRDS in external GDAL format to target directory +COPYRIGHT: (C) 2024 by Stefan Blumentrath and the GRASS Development Team + +This program is free software under the GNU General +Public License (>=v2). Read the file COPYING that +comes with GRASS for details. +""" + +# %module +# % description: Export raster maps from a STRDS to GeoTiff format in a temporal directory structure. +# % keyword: temporal +# % keyword: move +# % keyword: export +# % keyword: GDAL +# % keyword: directory +# %end + +# %option G_OPT_STRDS_INPUT +# %end + +# %option G_OPT_T_WHERE +# %end + +# %option G_OPT_M_DIR +# % key: output_directory +# % description: Path to the output / destination directory +# % required: yes +# %end + +# %option +# % key: format +# % type: string +# % description: Geotiff format flavor to export +# % options: COG,GTiff +# % answer: COG +# % required: no +# %end + +# %option +# % key: resampling +# % type: string +# % description: Resampling to be used for overviews +# % options: NEAREST,​AVERAGE,​BILINEAR​,​CUBIC​,​CUBICSPLINE​,​LANCZOS​,​MODE,​RMS +# % required: no +# %end + +# %option +# % key: level +# % type: integer +# % description: Compression level for the output files (0-22), GDAL default for ZSTD is 9, for DEFLATE is 6 +# % required: no +# %end + +# %option +# % key: over +# % type: integer +# % description: Number of overviews to create (4 is often a good choice) +# % required: no +# %end + +# %option +# % key: temporal_tree +# % type: string +# % description: Strftime format to create temporal directory name or tree (e.g. "%Y/%m/%d") +# % required: no +# %end + +# %option G_OPT_F_SEP +# %end + +# %option G_OPT_M_NPROCS +# %end + +# %flag +# % key: f +# % label: Force export of floating point maps as Float32 +# % description: Force export of floating point maps as Float32 +# %end + +# %flag +# % key: s +# % label: Use semantic label in directory structure +# % description: Use semantic label in directory structure +# %end + +import sys +from functools import partial +from multiprocessing import Pool +from pathlib import Path +from subprocess import PIPE + +import grass.script as gs +import grass.temporal as tgis +from grass.pygrass.modules.interface import Module + +OVERWRITE = False + + +def get_target_directory( + map_row: dict, + output_directory: Path | None = None, + temporal_tree: str = "%Y/%m/%d", + use_semantic_label: bool = False, +) -> Path: + """Get target directory tree for raster map.""" + if use_semantic_label: + output_directory /= map_row["semantic_label"] + return output_directory / map_row["start_time"].strftime(temporal_tree) + + +def check_datatype( + map_info: dict, + force_float32: bool = False, +) -> tuple[str, int | None]: + """Check the datatype of the map and return the smallest appropriate GDAL type. + + :param map_info: Dictionary containing map information + :param force_float32: Force the use of Float32 type for all floating point maps + """ + int_ranges = { + "Byte": (0, 255), + "UInt16": (0, 65535), + "Int16": (-32768, 32767), + "UInt32": (0, 4294967295), + "Int32": (-2147483648, 2147483647), + } + # Check for integer types + if map_info["datatype"] == "CELL": + # Check for integer types + # Int8 is not yet supported by r.out.gdal + for dtype, int_range in int_ranges.items(): + if map_info["min"] >= int_range[0] and map_info["max"] < int_range[1]: + return dtype, int_range[1] + if map_info["min"] > int_range[0] and map_info["max"] <= int_range[1]: + return dtype, int_range[0] + + # Check for floating point types + elif map_info["datatype"] == "FCELL" or force_float32: + return "Float32", None + return "Float64", None + + +def export_map_row( + map_row: dict, + output_directory: str = ".", + flags: str = "cm", + raster_format: str = "COG", + compression: str = "LZW", + blocksize: int = 256, # must be multiple of 16 + resampling: str | None = None, + level: int | None = None, + overviews: int | None = None, + separator: str = "|", + temporal_tree: str = "%Y/%m/%d", + truncate_float: bool = False, + use_semantic_label: bool = False, +) -> str: + """Export raster map using r.out.gdal. + + The function applies the most appropriate create options + based on the raster map type, range and GDAL format driver. + See: + https://gdal.org/en/stable/drivers/raster/cog.html + https://gdal.org/en/stable/drivers/raster/gtiff.html + """ + # Check the Raster map type and range + raster_info = gs.raster_info(map_row["name"]) + data_type, no_data = check_datatype(raster_info, force_float32=truncate_float) + target_directory = get_target_directory( + map_row, + Path(output_directory), + temporal_tree=temporal_tree, + use_semantic_label=use_semantic_label, + ) + + # Allways create BIGTIFFs + createopt = "BIGTIFF=YES" + # Set the resampling method + if resampling: + createopt += f",RESAMPLING={resampling}" + # Default to NEAREST for Byte, AVERAGE for others + elif data_type == "Byte": + createopt += ",RESAMPLING=NEAREST" + else: + createopt += ",RESAMPLING=CUBIC" + # Set the compression type and make sure a predictor is used + createopt += f",COMPRESS={compression}" + + # Set GTiff specific create options + if raster_format == "GTiff": + createopt += ",TILED=YES" + if compression in {"DEFLATE", "LERC_DEFLATE"} and level: + createopt += f",ZLEVEL={level}" + elif compression in {"ZSTD", "LERC_ZSTD"} and level: + createopt += f",ZSTD_LEVEL={level}" + if data_type in {"Float32", "Float64"}: + createopt += ",PREDICTOR=3" + else: + createopt += ",PREDICTOR=2" + createopt += f",BLOCKXSIZE={blocksize},BLOCKYSIZE={blocksize}" + else: + createopt += f",BLOCKSIZE={blocksize}" + # Statistics are useful for registering the map in an STRDS + if ( + compression in {"DEFLATE", "ZSTD", "LERC_DEFLATE", "LERC_ZSTD", "LZMA"} + and level + ): + createopt += f",LEVEL={level}" + createopt += ",PREDICTOR=YES,STATISTICS=YES" + + if data_type == "Float32" and truncate_float and "f" not in flags: + flags += "f" + + output_file = target_directory / f"{map_row['name']}.tif" + if output_file.exists() and not OVERWRITE: + gs.warning( + _( + "Output file {} already exists and overwrite is false, skipping...", + ).format(str(output_file)), + ) + else: + export_module = Module( + "r.out.gdal", + flags=flags, + input=map_row["name"], + output=str(output_file), + format=raster_format, + nodata=no_data, + type=data_type, + overwrite=OVERWRITE, + createopt=createopt, + overviews=overviews, + quiet=True, + stderr_=PIPE, + ) + if export_module.returncode != 0: + gs.fatal( + _( + "Exporting map {} failed with error:\n{}", + ).format(map_row["name"], export_module.stderr), + ) + return separator.join( + [ + map_row["name"], + map_row["start_time"].strftime("%Y-%m-%d %H:%M:%S"), + ( + map_row["end_time"].strftime("%Y-%m-%d %H:%M:%S") + if map_row["end_time"] + else "" + ), + map_row["semantic_label"] or "", + str(output_file), + ], + ) + + +def main() -> None: + """Do the main work.""" + options, flags = gs.parser() + global OVERWRITE + OVERWRITE = gs.overwrite() + nprocs = int(options["nprocs"]) + output_directory = Path(options["output_directory"]) + # Check if maps are exported to GDAL formats + temporal_tree = options["temporal_tree"] or "%Y/%m/%d" + + tgis.init() + input_strds = tgis.open_old_stds(options["input"], "strds") + input_strds_maps = input_strds.get_registered_maps( + columns="name,start_time,end_time,semantic_label", + where=options["where"], + ) + if not input_strds_maps: + gs.warning(_("No maps found in the space-time raster dataset.")) + return + + # Create the output directory structure if it does not exist + # Doing this sequentialy to avoid race condisions + output_directories = { + get_target_directory( + map_row, + output_directory, + temporal_tree=temporal_tree, + use_semantic_label=flags["s"], + ) + for map_row in input_strds_maps + } + for output_dir in output_directories: + output_dir.mkdir(parents=True, exist_ok=True) + export_map_row_tif = partial( + export_map_row, + output_directory=str(output_directory), + flags="cm", + raster_format=options["format"], + compression=options["compression"], + level=int(options["level"]) if options["level"] else None, + overviews=int(options["overviews"]) if options["overviews"] else None, + truncate_float=flags["f"], + separator=options["separator"], + temporal_tree=options["temporal_tree"] or "%Y/%m/%d", + use_semantic_label=flags["s"], + ) + + # Export maps to GeoTIFF / COG + if nprocs > 1: + with Pool(nprocs) as pool: + register_strings = pool.map( + export_map_row_tif, + (dict(input_map) for input_map in input_strds_maps), + ) + else: + register_strings = [export_map_row_tif(map_row) for map_row in input_strds_maps] + + # Print register information + print("\n".join(register_strings)) + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/temporal/t.rast.exporttree/testsuite/test_t_rast_exporttree.py b/src/temporal/t.rast.exporttree/testsuite/test_t_rast_exporttree.py new file mode 100644 index 0000000..4bc8f57 --- /dev/null +++ b/src/temporal/t.rast.exporttree/testsuite/test_t_rast_exporttree.py @@ -0,0 +1,150 @@ +"""Test t.rast.copytree + +(C) 2024 by NVE, Stefan Blumentrath and the GRASS GIS Development Team +This program is free software under the GNU General Public +License (>=v2). Read the file COPYING that comes with GRASS +for details. + +:authors: Stefan Blumentrath +""" + +# ruff: noqa: RUF012 + +from pathlib import Path + +import grass.script as gs +from grass.gunittest.case import TestCase + + +class TestTRastCopytree(TestCase): + """Test case for moving files from STRDS to directory trees""" + + default_region = {"s": 0, "n": 80, "w": 0, "e": 120, "b": 0, "t": 50} + + @classmethod + def setUpClass(cls): + """Initiate the temporal GIS and set the region""" + cls.use_temp_region() + cls.tempdir_target = Path(gs.tempdir()) + cls.runModule("g.region", **cls.default_region, res=1, res3=1) + + for rmap_idx in range(1, 4): + for prefix in ("a", "b", "c"): + cls.runModule( + "r.mapcalc", + expression=f"{prefix}_{rmap_idx} = {rmap_idx}00", + overwrite=True, + ) + cls.runModule( + "r.support", map=f"{prefix}_{rmap_idx}", semantic_label=prefix + ) + + cls.runModule( + "t.create", + type="strds", + temporaltype="absolute", + output="A", + title="A test", + description="A test", + overwrite=True, + ) + cls.runModule( + "t.register", + flags="i", + type="raster", + input="A", + maps="a_1,a_2,a_3", + start="2001-01-01", + increment="3 months", + overwrite=True, + ) + + cls.runModule( + "t.create", + type="strds", + temporaltype="absolute", + output="B", + title="B test", + description="B test", + overwrite=True, + ) + cls.runModule( + "t.register", + flags="i", + type="raster", + input="B", + maps="b_1,b_2,b_3", + start="2001-01-01", + increment="1 day", + overwrite=True, + ) + + cls.runModule( + "t.create", + type="strds", + temporaltype="absolute", + output="C", + title="C test", + description="C test", + overwrite=True, + ) + cls.runModule( + "t.register", + flags="i", + type="raster", + input="C", + maps="c_1,c_2,c_3", + start="2001-01-01", + increment="1 day", + overwrite=True, + ) + + @classmethod + def tearDownClass(cls): + """Remove the temporary region""" + cls.del_temp_region() + cls.runModule("t.remove", flags="df", type="strds", inputs="A") + gs.utils.try_rmdir(str(cls.tempdir_target)) + + def test_t_rast_copytree_move(self): + """Test moving files into directory tree""" + # Check that t.rast.copytree runs successfully + self.assertModule( + "t.rast.exporttree", + flags="f", + input="A", + temporal_tree="%Y/%m", + format="COG", + compression="ZSTD", + level=12, + output_directory=str(self.tempdir_target), + nprocs="2", + ) + + self.assertFileExists(str(self.tempdir_target / "2001/01/a_1.tif")) + self.assertFileExists(str(self.tempdir_target / "2001/04/a_2.tif")) + self.assertFileExists(str(self.tempdir_target / "2001/07/a_3.tif")) + + def test_t_rast_copytree_copy_semantic_label(self): + """Test moving files into directory tree""" + # Check that t.rast.copytree runs successfully + self.assertModule( + "t.rast.exporttree", + flags="s", + input="B", + format="GTiff", + compression="LZW", + temporal_tree="%Y/%m/%d", + output_directory=str(self.tempdir_target), + nprocs="2", + ) + + self.assertFileExists(str(self.tempdir_target / "b/2001/01/01/b_1.tif")) + self.assertFileExists(str(self.tempdir_target / "b/2001/01/02/b_2.tif")) + self.assertFileExists(str(self.tempdir_target / "b/2001/01/03/b_3.tif")) + + +if __name__ == "__main__": + from grass.gunittest.main import test + + test()