diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index 8dbea2478027..3b2c309994a4 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -21,6 +21,7 @@ add_library( gdalalg_raster_overview_delete.cpp gdalalg_raster_read.cpp gdalalg_raster_reproject.cpp + gdalalg_raster_stack.cpp gdalalg_raster_write.cpp gdalalg_vector.cpp gdalalg_vector_info.cpp diff --git a/apps/gdalalg_raster.cpp b/apps/gdalalg_raster.cpp index 84808b402f63..ba68d40d55ca 100644 --- a/apps/gdalalg_raster.cpp +++ b/apps/gdalalg_raster.cpp @@ -20,6 +20,7 @@ #include "gdalalg_raster_overview.h" #include "gdalalg_raster_pipeline.h" #include "gdalalg_raster_reproject.h" +#include "gdalalg_raster_stack.h" /************************************************************************/ /* GDALRasterAlgorithm */ @@ -47,6 +48,7 @@ class GDALRasterAlgorithm final : public GDALAlgorithm RegisterSubAlgorithm(); RegisterSubAlgorithm(); RegisterSubAlgorithm(); + RegisterSubAlgorithm(); } private: diff --git a/apps/gdalalg_raster_stack.cpp b/apps/gdalalg_raster_stack.cpp new file mode 100644 index 000000000000..6d8d0552ad6e --- /dev/null +++ b/apps/gdalalg_raster_stack.cpp @@ -0,0 +1,314 @@ +/****************************************************************************** + * + * Project: GDAL + * Purpose: gdal "raster stack" subcommand + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2024, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#include "gdalalg_raster_stack.h" + +#include "cpl_conv.h" +#include "cpl_vsi_virtual.h" + +#include "gdal_priv.h" +#include "gdal_utils.h" + +//! @cond Doxygen_Suppress + +#ifndef _ +#define _(x) (x) +#endif + +/************************************************************************/ +/* GDALRasterStackAlgorithm::GDALRasterStackAlgorithm() */ +/************************************************************************/ + +GDALRasterStackAlgorithm::GDALRasterStackAlgorithm() + : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL) +{ + AddProgressArg(); + AddOutputFormatArg(&m_format); + AddArg(GDAL_ARG_NAME_INPUT, 'i', + _("Input raster datasets (or specify a @ to point to a " + "file containing filenames)"), + &m_inputDatasets, GDAL_OF_RASTER) + .SetPositional() + .SetMinCount(1) + .SetAutoOpenDataset(false) + .SetMetaVar("INPUTS"); + AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER); + AddCreationOptionsArg(&m_creationOptions); + AddArg("band", 'b', _("Specify input band(s) number."), &m_bands); + AddOverwriteArg(&m_overwrite); + { + auto &arg = AddArg("resolution", 0, + _("Target resolution (in destination CRS units)"), + &m_resolution) + .SetMetaVar(",|average|highest|lowest"); + arg.AddValidationAction( + [this, &arg]() + { + const std::string val = arg.Get(); + if (val != "average" && val != "highest" && val != "lowest") + { + const auto aosTokens = + CPLStringList(CSLTokenizeString2(val.c_str(), ",", 0)); + if (aosTokens.size() != 2 || + CPLGetValueType(aosTokens[0]) == CPL_VALUE_STRING || + CPLGetValueType(aosTokens[1]) == CPL_VALUE_STRING || + CPLAtof(aosTokens[0]) <= 0 || + CPLAtof(aosTokens[1]) <= 0) + { + ReportError(CE_Failure, CPLE_AppDefined, + "resolution: two comma separated positive " + "values should be provided, or 'average', " + "'highest' or 'lowest'"); + return false; + } + } + return true; + }); + } + AddBBOXArg(&m_bbox, _("Target bounding box as xmin,ymin,xmax,ymax (in " + "destination CRS units)")); + AddArg("target-aligned-pixels", 0, + _("Round target extent to target resolution"), + &m_targetAlignedPixels) + .AddHiddenAlias("tap"); + AddArg("srcnodata", 0, _("Set nodata values for input bands."), + &m_srcNoData) + .SetMinCount(1) + .SetRepeatedArgAllowed(false); + AddArg("dstnodata", 0, + _("Set nodata values at the destination band level."), &m_dstNoData) + .SetMinCount(1) + .SetRepeatedArgAllowed(false); + AddArg("hidenodata", 0, + _("Makes the destination band not report the NoData."), + &m_hideNoData); +} + +/************************************************************************/ +/* GDALRasterStackAlgorithm::RunImpl() */ +/************************************************************************/ + +bool GDALRasterStackAlgorithm::RunImpl(GDALProgressFunc pfnProgress, + void *pProgressData) +{ + if (m_outputDataset.GetDatasetRef()) + { + ReportError(CE_Failure, CPLE_NotSupported, + "gdal raster stack does not support outputting to an " + "already opened output dataset"); + return false; + } + + std::vector ahInputDatasets; + CPLStringList aosInputDatasetNames; + bool foundByRef = false; + bool foundByName = false; + for (auto &ds : m_inputDatasets) + { + if (ds.GetDatasetRef()) + { + foundByRef = true; + ahInputDatasets.push_back( + GDALDataset::ToHandle(ds.GetDatasetRef())); + } + else if (!ds.GetName().empty()) + { + foundByName = true; + if (ds.GetName()[0] == '@') + { + auto f = VSIVirtualHandleUniquePtr( + VSIFOpenL(ds.GetName().c_str() + 1, "r")); + if (!f) + { + ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s", + ds.GetName().c_str() + 1); + return false; + } + while (const char *filename = CPLReadLineL(f.get())) + { + aosInputDatasetNames.push_back(filename); + } + } + else if (ds.GetName().find_first_of("*?[") != std::string::npos) + { + CPLStringList aosMatches(VSIGlob(ds.GetName().c_str(), nullptr, + pfnProgress, pProgressData)); + for (const char *pszStr : aosMatches) + { + aosInputDatasetNames.push_back(pszStr); + } + } + else + { + aosInputDatasetNames.push_back(ds.GetName().c_str()); + } + } + } + if (foundByName && foundByRef) + { + ReportError(CE_Failure, CPLE_NotSupported, + "Input datasets should be provided either all by reference " + "or all by name"); + return false; + } + + VSIStatBufL sStat; + if (!m_overwrite && !m_outputDataset.GetName().empty() && + (VSIStatL(m_outputDataset.GetName().c_str(), &sStat) == 0 || + std::unique_ptr( + GDALDataset::Open(m_outputDataset.GetName().c_str())))) + { + ReportError(CE_Failure, CPLE_AppDefined, + "File '%s' already exists. Specify the --overwrite " + "option to overwrite it.", + m_outputDataset.GetName().c_str()); + return false; + } + + const bool bVRTOutput = + m_outputDataset.GetName().empty() || EQUAL(m_format.c_str(), "VRT") || + EQUAL(CPLGetExtensionSafe(m_outputDataset.GetName().c_str()).c_str(), + "VRT"); + + CPLStringList aosOptions; + + aosOptions.push_back("-separate"); + + if (!m_resolution.empty()) + { + const auto aosTokens = + CPLStringList(CSLTokenizeString2(m_resolution.c_str(), ",", 0)); + if (aosTokens.size() == 2) + { + aosOptions.push_back("-tr"); + aosOptions.push_back(aosTokens[0]); + aosOptions.push_back(aosTokens[1]); + } + else + { + aosOptions.push_back("-resolution"); + aosOptions.push_back(m_resolution); + } + } + if (!m_bbox.empty()) + { + aosOptions.push_back("-te"); + aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[0])); + aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[1])); + aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[2])); + aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[3])); + } + if (m_targetAlignedPixels) + { + aosOptions.push_back("-tap"); + } + if (!m_srcNoData.empty()) + { + aosOptions.push_back("-srcnodata"); + std::string s; + for (double v : m_srcNoData) + { + if (!s.empty()) + s += " "; + s += CPLSPrintf("%.17g", v); + } + aosOptions.push_back(s); + } + if (!m_dstNoData.empty()) + { + aosOptions.push_back("-vrtnodata"); + std::string s; + for (double v : m_dstNoData) + { + if (!s.empty()) + s += " "; + s += CPLSPrintf("%.17g", v); + } + aosOptions.push_back(s); + } + if (bVRTOutput) + { + for (const auto &co : m_creationOptions) + { + aosOptions.push_back("-co"); + aosOptions.push_back(co); + } + } + for (const int b : m_bands) + { + aosOptions.push_back("-b"); + aosOptions.push_back(CPLSPrintf("%d", b)); + } + if (m_hideNoData) + { + aosOptions.push_back("-hidenodata"); + } + + GDALBuildVRTOptions *psOptions = + GDALBuildVRTOptionsNew(aosOptions.List(), nullptr); + if (bVRTOutput) + { + GDALBuildVRTOptionsSetProgress(psOptions, pfnProgress, pProgressData); + } + + auto poOutDS = std::unique_ptr(GDALDataset::FromHandle( + GDALBuildVRT(bVRTOutput ? m_outputDataset.GetName().c_str() : "", + foundByName ? aosInputDatasetNames.size() + : static_cast(m_inputDatasets.size()), + ahInputDatasets.empty() ? nullptr : ahInputDatasets.data(), + aosInputDatasetNames.List(), psOptions, nullptr))); + GDALBuildVRTOptionsFree(psOptions); + bool bOK = poOutDS != nullptr; + if (bOK) + { + if (bVRTOutput) + { + m_outputDataset.Set(std::move(poOutDS)); + } + else + { + CPLStringList aosTranslateOptions; + if (!m_format.empty()) + { + aosTranslateOptions.AddString("-of"); + aosTranslateOptions.AddString(m_format.c_str()); + } + for (const auto &co : m_creationOptions) + { + aosTranslateOptions.AddString("-co"); + aosTranslateOptions.AddString(co.c_str()); + } + + GDALTranslateOptions *psTranslateOptions = + GDALTranslateOptionsNew(aosTranslateOptions.List(), nullptr); + GDALTranslateOptionsSetProgress(psTranslateOptions, pfnProgress, + pProgressData); + + auto poFinalDS = + std::unique_ptr(GDALDataset::FromHandle( + GDALTranslate(m_outputDataset.GetName().c_str(), + GDALDataset::ToHandle(poOutDS.get()), + psTranslateOptions, nullptr))); + GDALTranslateOptionsFree(psTranslateOptions); + + bOK = poFinalDS != nullptr; + if (bOK) + { + m_outputDataset.Set(std::move(poFinalDS)); + } + } + } + + return bOK; +} + +//! @endcond diff --git a/apps/gdalalg_raster_stack.h b/apps/gdalalg_raster_stack.h new file mode 100644 index 000000000000..d60fe9b8c1ca --- /dev/null +++ b/apps/gdalalg_raster_stack.h @@ -0,0 +1,59 @@ +/****************************************************************************** + * + * Project: GDAL + * Purpose: gdal "raster stack" subcommand + * Author: Even Rouault + * + ****************************************************************************** + * Copyright (c) 2025, Even Rouault + * + * SPDX-License-Identifier: MIT + ****************************************************************************/ + +#ifndef GDALALG_RASTER_STACK_INCLUDED +#define GDALALG_RASTER_STACK_INCLUDED + +#include "gdalalgorithm.h" + +//! @cond Doxygen_Suppress + +/************************************************************************/ +/* GDALRasterStackAlgorithm */ +/************************************************************************/ + +class GDALRasterStackAlgorithm final : public GDALAlgorithm +{ + public: + static constexpr const char *NAME = "stack"; + static constexpr const char *DESCRIPTION = + "Combine together input bands into a multi-band output, either virtual " + "(VRT) or materialized."; + static constexpr const char *HELP_URL = "/programs/gdal_raster_stack.html"; + + static std::vector GetAliases() + { + return {}; + } + + explicit GDALRasterStackAlgorithm(); + + private: + bool RunImpl(GDALProgressFunc pfnProgress, void *pProgressData) override; + + std::vector m_inputDatasets{}; + std::string m_format{}; + GDALArgDatasetValue m_outputDataset{}; + std::vector m_creationOptions{}; + bool m_overwrite = false; + std::string m_resolution{}; + std::vector m_bbox{}; + bool m_targetAlignedPixels = false; + std::vector m_srcNoData{}; + std::vector m_dstNoData{}; + std::vector m_bands{}; + bool m_hideNoData = false; +}; + +//! @endcond + +#endif diff --git a/autotest/utilities/test_gdalalg_raster_mosaic.py b/autotest/utilities/test_gdalalg_raster_mosaic.py index b5d6430e40c3..bf526b9b5176 100755 --- a/autotest/utilities/test_gdalalg_raster_mosaic.py +++ b/autotest/utilities/test_gdalalg_raster_mosaic.py @@ -86,32 +86,6 @@ def test_gdalalg_raster_mosaic_overwrite(tmp_vsimem): assert ds.GetRasterBand(1).Checksum() == 4672 -def test_gdalalg_raster_mosaic_separate(tmp_vsimem): - - tmp_filename = str(tmp_vsimem / "tmp.tif") - gdal.Translate(tmp_filename, "../gcore/data/byte.tif", options="-scale 0 255 255 0") - - alg = get_mosaic_alg() - assert alg.ParseCommandLineArguments( - [ - "--separate", - "../gcore/data/byte.tif", - tmp_filename, - "", - ] - ) - assert alg.Run() - ds = alg.GetArg("output").Get().GetDataset() - assert ds.RasterCount == 2 - assert ds.RasterXSize == 20 - assert ds.RasterYSize == 20 - assert ds.GetGeoTransform() == pytest.approx( - (440720.0, 60.0, 0.0, 3751320.0, 0.0, -60.0) - ) - assert ds.GetRasterBand(1).Checksum() == 4672 - assert ds.GetRasterBand(2).Checksum() == 4563 - - def test_gdalalg_raster_mosaic_bbox(): alg = get_mosaic_alg() diff --git a/doc/source/conf.py b/doc/source/conf.py index e6b42c6dc9e3..952b1c723073 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -278,6 +278,13 @@ [author_evenr], 1, ), + ( + "programs/gdal_raster_stack", + "gdal-raster-stack", + "Combine together input bands into a multi-band output, either virtual (VRT) or materialized", + [author_evenr], + 1, + ), ( "programs/gdal_vector", "gdal-vector", diff --git a/doc/source/programs/gdal_raster.rst b/doc/source/programs/gdal_raster.rst index d3f90fa5c05c..4f6190896d98 100644 --- a/doc/source/programs/gdal_raster.rst +++ b/doc/source/programs/gdal_raster.rst @@ -19,13 +19,16 @@ Synopsis Usage: gdal raster where is one of: - - clip: Clip a raster dataset. - - convert: Convert a raster dataset. - - info: Return information on a raster dataset. - - mosaic: Build a mosaic, either virtual (VRT) or materialized. + - clip: Clip a raster dataset. + - convert: Convert a raster dataset. + - edit: Edit a raster dataset. + - info: Return information on a raster dataset. + - mosaic: Build a mosaic, either virtual (VRT) or materialized. - overview: Manage overviews of a raster dataset. - pipeline: Process a raster dataset. - reproject: Reproject a raster dataset. + - stack: Combine together input bands into a multi-band output, either virtual (VRT) or materialized. + Available sub-commands ---------------------- @@ -37,6 +40,7 @@ Available sub-commands - :ref:`gdal_raster_overview_subcommand` - :ref:`gdal_raster_pipeline_subcommand` - :ref:`gdal_raster_reproject_subcommand` +- :ref:`gdal_raster_stack_subcommand` Examples -------- diff --git a/doc/source/programs/gdal_raster_mosaic.rst b/doc/source/programs/gdal_raster_mosaic.rst index 82c1ae28b689..332bc3a45e65 100644 --- a/doc/source/programs/gdal_raster_mosaic.rst +++ b/doc/source/programs/gdal_raster_mosaic.rst @@ -37,7 +37,6 @@ Synopsis -f, --of, --format, --output-format Output format --co, --creation-option = Creation option [may be repeated] -b, --band Specify input band(s) number. [may be repeated] - --separate Place each input file into a separate band. --overwrite Whether overwriting existing output is allowed --resolution ,|average|highest|lowest> Target resolution (in destination CRS units) --bbox Target bounding box as xmin,ymin,xmax,ymax (in destination CRS units) @@ -61,10 +60,6 @@ of the command line or put in a text file (one filename per line) for very long Wildcards '*', '?' or '['] of :cpp:func:`VSIGlob` can be used, even on files located on network file systems such as /vsis3/, /vsigs/, /vsiaz/, etc. -With :option:`--separate`, each input goes into a separate band in the output dataset. Otherwise, -the files are considered as source rasters of a larger mosaic and the output file has the same number of -bands as the input files. - :program:`gdal raster mosaic` does some checks to ensure that all files that will be put in the resulting file have similar characteristics: number of bands, projection, color interpretation, etc. If not, files that do not match the common characteristics will be skipped. @@ -92,15 +87,6 @@ The following options are available: If input bands not set all bands will be added to the output. Multiple :option:`-b` switches may be used to select a set of input bands. -.. option:: --separate - - Place each input file into a separate band. See :example:`separate`. - Contrary to the default mode, it is not - required that all bands have the same datatype. - - All bands of each input file are added as separate - output bands, unless :option:`-b` is specified to select a subset of them. - .. option:: --resolution {|highest|lowest|average} In case the resolution of all input files is not the same, the :option:`--resolution` flag @@ -166,14 +152,6 @@ The following options are available: Examples -------- -.. example:: - :title: Make a RGB virtual mosaic from 3 single-band input files - :id: separate - - .. code-block:: bash - - gdal raster mosaic --separate red.tif green.tif blue.tif rgb.tif - .. example:: :title: Make a virtual mosaic with blue background colour (RGB: 0 0 255) :id: dstnodata diff --git a/doc/source/programs/gdal_raster_stack.rst b/doc/source/programs/gdal_raster_stack.rst new file mode 100644 index 000000000000..a43a702ed6a7 --- /dev/null +++ b/doc/source/programs/gdal_raster_stack.rst @@ -0,0 +1,143 @@ +.. _gdal_raster_stack_subcommand: + +================================================================================ +"gdal raster stack" sub-command +================================================================================ + +.. versionadded:: 3.11 + +.. only:: html + + Combine together input bands into a multi-band output, either virtual (VRT) or materialized. + + +.. Index:: gdal raster stack + +Synopsis +-------- + +.. code-block:: + + Usage: gdal raster stack [OPTIONS] + + Combine together input bands into a multi-band output, either virtual (VRT) or materialized. + + Positional arguments: + -i, --input Input raster datasets (or specify a @ to point to a file containing filenames) [1.. values] + -o, --output Output raster dataset (created by algorithm) [required] + + Common Options: + -h, --help Display help message and exit + --version Display GDAL version and exit + --json-usage Display usage as JSON document and exit + --drivers Display driver list as JSON document and exit + --config = Configuration option [may be repeated] + --progress Display progress bar + + Options: + -f, --of, --format, --output-format Output format + --co, --creation-option = Creation option [may be repeated] + -b, --band Specify input band(s) number. [may be repeated] + --overwrite Whether overwriting existing output is allowed + --resolution ,|average|highest|lowest> Target resolution (in destination CRS units) + --bbox Target bounding box as xmin,ymin,xmax,ymax (in destination CRS units) + --target-aligned-pixels Round target extent to target resolution + --srcnodata Set nodata values for input bands. [1.. values] + --dstnodata Set nodata values at the destination band level. [1.. values] + --hidenodata Makes the destination band not report the NoData. + + + +Description +----------- + +This program combine together input bands from GDAL datasets into a dataset, that can be +either a virtual stack in the :ref:`VRT (Virtual Dataset) ` format, +or in a more conventional raster format such as GeoTIFF. + +All bands of each input file are added as separate +output bands, unless :option:`-b` is specified to select a subset of them. + +The list of input GDAL datasets can be specified at the end +of the command line or put in a text file (one filename per line) for very long lists. +Wildcards '*', '?' or '['] of :cpp:func:`VSIGlob` can be used, even on files located +on network file systems such as /vsis3/, /vsigs/, /vsiaz/, etc. + +The following options are available: + +.. include:: gdal_options/of_raster_create_copy.rst + +.. include:: gdal_options/co.rst + +.. include:: gdal_options/overwrite.rst + +.. option:: -b + + Select an input to be processed. Bands are numbered from 1. + If input bands not set all bands will be added to the output. + Multiple :option:`-b` switches may be used to select a set of input bands. + +.. option:: --resolution {|highest|lowest|average} + + In case the resolution of all input files is not the same, the :option:`--resolution` flag + enables the user to control the way the output resolution is computed. + + `highest` will pick the smallest values of pixel dimensions within the set of source rasters. + + `lowest` will pick the largest values of pixel dimensions within the set of source rasters. + + `average` is the default and will compute an average of pixel dimensions within the set of source rasters. + + ,. The values must be expressed in georeferenced units. + Both must be positive values. + +.. option:: --bbox ,,, + + Set georeferenced extents of output file. The values must be expressed in georeferenced units. + If not specified, the extent of the output is the minimum bounding box of the set of source rasters. + Pixels within the extent of the output but not covered by a source raster will be read as valid + pixels with a value of zero unless a NODATA value is specified using :option:`--dstnodata` + or an alpha mask band is added with :option:`--addalpha`. + +.. option:: --target-aligned-pixels + + (target aligned pixels) align + the coordinates of the extent of the output file to the values of the :option:`--resolution`, + such that the aligned extent includes the minimum extent. + Alignment means that xmin / resx, ymin / resy, xmax / resx and ymax / resy are integer values. + +.. option:: --srcnodata [,]... + + Set nodata values for input bands (different values can be supplied for each band). + If the option is not specified, the intrinsic nodata settings on the source datasets + will be used (if they exist). The value set by this option is written in the NODATA element + of each ``ComplexSource`` element. + +.. option:: --dstnodata [,]... + + Set nodata values at the output band level (different values can be supplied for each band). If more + than one value is supplied, all values should be quoted to keep them together + as a single operating system argument. If the option is not specified, + intrinsic nodata settings on the first dataset will be used (if they exist). The value set by this option + is written in the ``NoDataValue`` element of each ``VRTRasterBand element``. Use a value of + `None` to ignore intrinsic nodata settings on the source datasets. + +.. option:: --hidenodata + + Even if any band contains nodata value, giving this option makes the output band + not report the NoData. Useful when you want to control the background color of + the dataset. By using along with the -addalpha option, you can prepare a + dataset which doesn't report nodata value but is transparent in areas with no + data. + + +Examples +-------- + +.. example:: + :title: Make a RGB virtual stack from 3 single-band input files + :id: separate + + .. code-block:: bash + + gdal raster stack --separate red.tif green.tif blue.tif rgb.tif diff --git a/doc/source/programs/index.rst b/doc/source/programs/index.rst index a1fe6c7e103e..10a4f65e2291 100644 --- a/doc/source/programs/index.rst +++ b/doc/source/programs/index.rst @@ -39,6 +39,7 @@ single :program:`gdal` program that accepts commands and subcommands. gdal_raster_overview_delete gdal_raster_pipeline gdal_raster_reproject + gdal_raster_stack gdal_vector gdal_vector_info gdal_vector_clip @@ -61,6 +62,7 @@ single :program:`gdal` program that accepts commands and subcommands. - :ref:`gdal_raster_overview_delete_subcommand`: Remove overviews of a raster dataset - :ref:`gdal_raster_pipeline_subcommand`: Process a raster dataset - :ref:`gdal_raster_reproject_subcommand`: Reproject a raster dataset + - :ref:`gdal_raster_stack_subcommand`: Combine together input bands into a multi-band output, either virtual (VRT) or materialized. - :ref:`gdal_vector_command`: Entry point for vector commands - :ref:`gdal_vector_info_subcommand`: Get information on a vector dataset - :ref:`gdal_vector_clip_subcommand`: Clip a vector dataset