From 48a7064eb2209b2ae65b1104674e0180b163ac91 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 24 Jan 2025 16:47:01 +0100 Subject: [PATCH 1/2] Rename 'gdal raster buildvrt' to 'gdal raster mosaic', and allow both virtual and materialized output --- apps/CMakeLists.txt | 2 +- apps/gdalalg_raster.cpp | 4 +- ...buildvrt.cpp => gdalalg_raster_mosaic.cpp} | 97 ++++++++---- ...ter_buildvrt.h => gdalalg_raster_mosaic.h} | 24 +-- ...ldvrt.py => test_gdalalg_raster_mosaic.py} | 138 +++++++++++------- doc/source/conf.py | 14 +- doc/source/programs/gdal_raster.rst | 4 +- ...er_buildvrt.rst => gdal_raster_mosaic.rst} | 105 +++++++------ doc/source/programs/index.rst | 4 +- 9 files changed, 241 insertions(+), 151 deletions(-) rename apps/{gdalalg_raster_buildvrt.cpp => gdalalg_raster_mosaic.cpp} (74%) rename apps/{gdalalg_raster_buildvrt.h => gdalalg_raster_mosaic.h} (67%) rename autotest/utilities/{test_gdalalg_raster_buildvrt.py => test_gdalalg_raster_mosaic.py} (76%) rename doc/source/programs/{gdal_raster_buildvrt.rst => gdal_raster_mosaic.rst} (54%) diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index a89f8c456f15..8dbea2478027 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -11,11 +11,11 @@ add_library( gdalalg_main.cpp gdalalg_pipeline.cpp gdalalg_raster.cpp - gdalalg_raster_buildvrt.cpp gdalalg_raster_info.cpp gdalalg_raster_clip.cpp gdalalg_raster_convert.cpp gdalalg_raster_edit.cpp + gdalalg_raster_mosaic.cpp gdalalg_raster_pipeline.cpp gdalalg_raster_overview_add.cpp gdalalg_raster_overview_delete.cpp diff --git a/apps/gdalalg_raster.cpp b/apps/gdalalg_raster.cpp index a77e02ddb903..84808b402f63 100644 --- a/apps/gdalalg_raster.cpp +++ b/apps/gdalalg_raster.cpp @@ -12,11 +12,11 @@ #include "gdalalgorithm.h" -#include "gdalalg_raster_buildvrt.h" #include "gdalalg_raster_info.h" #include "gdalalg_raster_clip.h" #include "gdalalg_raster_convert.h" #include "gdalalg_raster_edit.h" +#include "gdalalg_raster_mosaic.h" #include "gdalalg_raster_overview.h" #include "gdalalg_raster_pipeline.h" #include "gdalalg_raster_reproject.h" @@ -46,7 +46,7 @@ class GDALRasterAlgorithm final : public GDALAlgorithm RegisterSubAlgorithm(); RegisterSubAlgorithm(); RegisterSubAlgorithm(); - RegisterSubAlgorithm(); + RegisterSubAlgorithm(); } private: diff --git a/apps/gdalalg_raster_buildvrt.cpp b/apps/gdalalg_raster_mosaic.cpp similarity index 74% rename from apps/gdalalg_raster_buildvrt.cpp rename to apps/gdalalg_raster_mosaic.cpp index e8d32df67e5d..53e388fe15d7 100644 --- a/apps/gdalalg_raster_buildvrt.cpp +++ b/apps/gdalalg_raster_mosaic.cpp @@ -1,7 +1,7 @@ /****************************************************************************** * * Project: GDAL - * Purpose: gdal "raster buildvrt" subcommand + * Purpose: gdal "raster mosaic" subcommand * Author: Even Rouault * ****************************************************************************** @@ -10,7 +10,7 @@ * SPDX-License-Identifier: MIT ****************************************************************************/ -#include "gdalalg_raster_buildvrt.h" +#include "gdalalg_raster_mosaic.h" #include "cpl_conv.h" #include "cpl_vsi_virtual.h" @@ -25,13 +25,14 @@ #endif /************************************************************************/ -/* GDALRasterBuildVRTAlgorithm::GDALRasterBuildVRTAlgorithm() */ +/* GDALRasterMosaicAlgorithm::GDALRasterMosaicAlgorithm() */ /************************************************************************/ -GDALRasterBuildVRTAlgorithm::GDALRasterBuildVRTAlgorithm() +GDALRasterMosaicAlgorithm::GDALRasterMosaicAlgorithm() : 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)"), @@ -43,8 +44,6 @@ GDALRasterBuildVRTAlgorithm::GDALRasterBuildVRTAlgorithm() AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER); AddCreationOptionsArg(&m_creationOptions); AddArg("band", 'b', _("Specify input band(s) number."), &m_bands); - AddArg("separate", 0, _("Place each input file into a separate band."), - &m_separate); AddOverwriteArg(&m_overwrite); { auto &arg = AddArg("resolution", 0, @@ -85,29 +84,31 @@ GDALRasterBuildVRTAlgorithm::GDALRasterBuildVRTAlgorithm() &m_srcNoData) .SetMinCount(1) .SetRepeatedArgAllowed(false); - AddArg("vrtnodata", 0, _("Set nodata values at the VRT band level."), - &m_vrtNoData) + AddArg("dstnodata", 0, + _("Set nodata values at the destination band level."), &m_dstNoData) .SetMinCount(1) .SetRepeatedArgAllowed(false); - AddArg("hidenodata", 0, _("Makes the VRT band not report the NoData."), + AddArg("hidenodata", 0, + _("Makes the destination band not report the NoData."), &m_hideNoData); AddArg("addalpha", 0, - _("Adds an alpha mask band to the VRT when the source raster have " + _("Adds an alpha mask band to the destination when the source " + "raster have " "none."), &m_addAlpha); } /************************************************************************/ -/* GDALRasterBuildVRTAlgorithm::RunImpl() */ +/* GDALRasterMosaicAlgorithm::RunImpl() */ /************************************************************************/ -bool GDALRasterBuildVRTAlgorithm::RunImpl(GDALProgressFunc pfnProgress, - void *pProgressData) +bool GDALRasterMosaicAlgorithm::RunImpl(GDALProgressFunc pfnProgress, + void *pProgressData) { if (m_outputDataset.GetDatasetRef()) { ReportError(CE_Failure, CPLE_NotSupported, - "gdal raster buildvrt does not support outputting to an " + "gdal raster mosaic does not support outputting to an " "already opened output dataset"); return false; } @@ -178,6 +179,11 @@ bool GDALRasterBuildVRTAlgorithm::RunImpl(GDALProgressFunc pfnProgress, 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; if (!m_resolution.empty()) { @@ -219,11 +225,11 @@ bool GDALRasterBuildVRTAlgorithm::RunImpl(GDALProgressFunc pfnProgress, } aosOptions.push_back(s); } - if (!m_vrtNoData.empty()) + if (!m_dstNoData.empty()) { aosOptions.push_back("-vrtnodata"); std::string s; - for (double v : m_vrtNoData) + for (double v : m_dstNoData) { if (!s.empty()) s += " "; @@ -231,14 +237,13 @@ bool GDALRasterBuildVRTAlgorithm::RunImpl(GDALProgressFunc pfnProgress, } aosOptions.push_back(s); } - if (m_separate) - { - aosOptions.push_back("-separate"); - } - for (const auto &co : m_creationOptions) + if (bVRTOutput) { - aosOptions.push_back("-co"); - aosOptions.push_back(co); + for (const auto &co : m_creationOptions) + { + aosOptions.push_back("-co"); + aosOptions.push_back(co); + } } for (const int b : m_bands) { @@ -256,19 +261,57 @@ bool GDALRasterBuildVRTAlgorithm::RunImpl(GDALProgressFunc pfnProgress, GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(aosOptions.List(), nullptr); - GDALBuildVRTOptionsSetProgress(psOptions, pfnProgress, pProgressData); + if (bVRTOutput) + { + GDALBuildVRTOptionsSetProgress(psOptions, pfnProgress, pProgressData); + } auto poOutDS = std::unique_ptr(GDALDataset::FromHandle( - GDALBuildVRT(m_outputDataset.GetName().c_str(), + 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); - const bool bOK = poOutDS != nullptr; + bool bOK = poOutDS != nullptr; if (bOK) { - m_outputDataset.Set(std::move(poOutDS)); + 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; diff --git a/apps/gdalalg_raster_buildvrt.h b/apps/gdalalg_raster_mosaic.h similarity index 67% rename from apps/gdalalg_raster_buildvrt.h rename to apps/gdalalg_raster_mosaic.h index d01a562431eb..fbdbf7fed6bc 100644 --- a/apps/gdalalg_raster_buildvrt.h +++ b/apps/gdalalg_raster_mosaic.h @@ -1,7 +1,7 @@ /****************************************************************************** * * Project: GDAL - * Purpose: gdal "raster buildvrt" subcommand + * Purpose: gdal "raster mosaic" subcommand * Author: Even Rouault * ****************************************************************************** @@ -10,45 +10,45 @@ * SPDX-License-Identifier: MIT ****************************************************************************/ -#ifndef GDALALG_RASTER_BUILDVRT_INCLUDED -#define GDALALG_RASTER_BUILDVRT_INCLUDED +#ifndef GDALALG_RASTER_MOSAIC_INCLUDED +#define GDALALG_RASTER_MOSAIC_INCLUDED #include "gdalalgorithm.h" //! @cond Doxygen_Suppress /************************************************************************/ -/* GDALRasterBuildVRTAlgorithm */ +/* GDALRasterMosaicAlgorithm */ /************************************************************************/ -class GDALRasterBuildVRTAlgorithm final : public GDALAlgorithm +class GDALRasterMosaicAlgorithm final : public GDALAlgorithm { public: - static constexpr const char *NAME = "buildvrt"; - static constexpr const char *DESCRIPTION = "Build a virtual dataset (VRT)."; - static constexpr const char *HELP_URL = - "/programs/gdal_raster_buildvrt.html"; + static constexpr const char *NAME = "mosaic"; + static constexpr const char *DESCRIPTION = + "Build a mosaic, either virtual (VRT) or materialized."; + static constexpr const char *HELP_URL = "/programs/gdal_raster_mosaic.html"; static std::vector GetAliases() { return {}; } - explicit GDALRasterBuildVRTAlgorithm(); + explicit GDALRasterMosaicAlgorithm(); 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; - bool m_separate = false; std::string m_resolution{}; std::vector m_bbox{}; bool m_targetAlignedPixels = false; std::vector m_srcNoData{}; - std::vector m_vrtNoData{}; + std::vector m_dstNoData{}; std::vector m_bands{}; bool m_hideNoData = false; bool m_addAlpha = false; diff --git a/autotest/utilities/test_gdalalg_raster_buildvrt.py b/autotest/utilities/test_gdalalg_raster_mosaic.py similarity index 76% rename from autotest/utilities/test_gdalalg_raster_buildvrt.py rename to autotest/utilities/test_gdalalg_raster_mosaic.py index 5a944567ffa7..b5d6430e40c3 100755 --- a/autotest/utilities/test_gdalalg_raster_buildvrt.py +++ b/autotest/utilities/test_gdalalg_raster_mosaic.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- ############################################################################### # Project: GDAL/OGR Test Suite -# Purpose: 'gdal raster buildvrt' testing +# Purpose: 'gdal raster mosaic' testing # Author: Even Rouault # ############################################################################### @@ -16,15 +16,15 @@ from osgeo import gdal -def get_buildvrt_alg(): +def get_mosaic_alg(): reg = gdal.GetGlobalAlgorithmRegistry() raster = reg.InstantiateAlg("raster") - return raster.InstantiateSubAlgorithm("buildvrt") + return raster.InstantiateSubAlgorithm("mosaic") -def test_gdalalg_raster_buildvrt_from_dataset_handle(): +def test_gdalalg_raster_mosaic_from_dataset_handle(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set( [ gdal.Translate( @@ -46,9 +46,9 @@ def test_gdalalg_raster_buildvrt_from_dataset_handle(): assert ds.GetRasterBand(1).Checksum() == 4672 -def test_gdalalg_raster_buildvrt_from_dataset_name(): +def test_gdalalg_raster_mosaic_from_dataset_name(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["../gcore/data/byte.tif"]) alg.GetArg("output").Set("") assert alg.Run() @@ -61,23 +61,23 @@ def test_gdalalg_raster_buildvrt_from_dataset_name(): assert ds.GetRasterBand(1).Checksum() == 4672 -def test_gdalalg_raster_buildvrt_overwrite(tmp_vsimem): +def test_gdalalg_raster_mosaic_overwrite(tmp_vsimem): - out_filename = str(tmp_vsimem / "out.tif") + out_filename = str(tmp_vsimem / "out.vrt") - alg = get_buildvrt_alg() + alg = get_mosaic_alg() assert alg.ParseRunAndFinalize(["../gcore/data/utmsmall.tif", out_filename]) with gdal.Open(out_filename) as ds: assert ds.GetRasterBand(1).Checksum() == 50054 - alg = get_buildvrt_alg() + alg = get_mosaic_alg() with pytest.raises( Exception, match="already exists. Specify the --overwrite option" ): alg.ParseRunAndFinalize(["../gcore/data/byte.tif", out_filename]) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() assert alg.ParseRunAndFinalize( ["--overwrite", "../gcore/data/byte.tif", out_filename] ) @@ -86,12 +86,12 @@ def test_gdalalg_raster_buildvrt_overwrite(tmp_vsimem): assert ds.GetRasterBand(1).Checksum() == 4672 -def test_gdalalg_raster_buildvrt_separate(tmp_vsimem): +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_buildvrt_alg() + alg = get_mosaic_alg() assert alg.ParseCommandLineArguments( [ "--separate", @@ -112,9 +112,9 @@ def test_gdalalg_raster_buildvrt_separate(tmp_vsimem): assert ds.GetRasterBand(2).Checksum() == 4563 -def test_gdalalg_raster_buildvrt_bbox(): +def test_gdalalg_raster_mosaic_bbox(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() assert alg.ParseCommandLineArguments( [ "--bbox=440780.0,3750180.0,441860.0,3751260.0", @@ -132,14 +132,14 @@ def test_gdalalg_raster_buildvrt_bbox(): assert ds.GetRasterBand(1).Checksum() == 3695 -def test_gdalalg_raster_buildvrt_resolution_average(): +def test_gdalalg_raster_mosaic_resolution_average(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src2_ds = gdal.GetDriverByName("MEM").Create("", 2, 2) src2_ds.SetGeoTransform([3, 0.5, 0, 49, 0, -0.5]) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds, src2_ds]) alg.GetArg("output").Set("") assert alg.Run() @@ -149,14 +149,14 @@ def test_gdalalg_raster_buildvrt_resolution_average(): assert ds.GetGeoTransform() == pytest.approx((2.0, 0.75, 0.0, 49.0, 0.0, -0.75)) -def test_gdalalg_raster_buildvrt_resolution_highest(): +def test_gdalalg_raster_mosaic_resolution_highest(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src2_ds = gdal.GetDriverByName("MEM").Create("", 2, 2) src2_ds.SetGeoTransform([3, 0.5, 0, 49, 0, -0.5]) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds, src2_ds]) alg.GetArg("output").Set("") alg.GetArg("resolution").Set("highest") @@ -167,14 +167,14 @@ def test_gdalalg_raster_buildvrt_resolution_highest(): assert ds.GetGeoTransform() == pytest.approx((2.0, 0.5, 0.0, 49.0, 0.0, -0.5)) -def test_gdalalg_raster_buildvrt_resolution_lowest(): +def test_gdalalg_raster_mosaic_resolution_lowest(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src2_ds = gdal.GetDriverByName("MEM").Create("", 2, 2) src2_ds.SetGeoTransform([3, 0.5, 0, 49, 0, -0.5]) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds, src2_ds]) alg.GetArg("output").Set("") alg.GetArg("resolution").Set("lowest") @@ -185,14 +185,14 @@ def test_gdalalg_raster_buildvrt_resolution_lowest(): assert ds.GetGeoTransform() == pytest.approx((2.0, 1.0, 0.0, 49.0, 0.0, -1.0)) -def test_gdalalg_raster_buildvrt_resolution_custom(): +def test_gdalalg_raster_mosaic_resolution_custom(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src2_ds = gdal.GetDriverByName("MEM").Create("", 2, 2) src2_ds.SetGeoTransform([3, 0.5, 0, 49, 0, -0.5]) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds, src2_ds]) alg.GetArg("output").Set("") alg.GetArg("resolution").Set("0.5,1") @@ -203,14 +203,14 @@ def test_gdalalg_raster_buildvrt_resolution_custom(): assert ds.GetGeoTransform() == pytest.approx((2.0, 0.5, 0.0, 49.0, 0.0, -1.0)) -def test_gdalalg_raster_buildvrt_target_aligned_pixels(): +def test_gdalalg_raster_mosaic_target_aligned_pixels(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src2_ds = gdal.GetDriverByName("MEM").Create("", 2, 2) src2_ds.SetGeoTransform([3, 0.5, 0, 49, 0, -0.5]) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds, src2_ds]) alg.GetArg("output").Set("") alg.GetArg("resolution").Set("0.3,0.6") @@ -222,9 +222,9 @@ def test_gdalalg_raster_buildvrt_target_aligned_pixels(): assert ds.GetGeoTransform() == pytest.approx((1.8, 0.3, 0.0, 49.2, 0.0, -0.6)) -def test_gdalalg_raster_buildvrt_resolution_invalid(): +def test_gdalalg_raster_mosaic_resolution_invalid(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() with pytest.raises( Exception, match="resolution: two comma separated positive values should be provided, or 'average', 'highest' or 'lowest'", @@ -244,34 +244,34 @@ def test_gdalalg_raster_buildvrt_resolution_invalid(): alg.GetArg("resolution").Set("-0.5,-0.5") -def test_gdalalg_raster_buildvrt_srcnodata_vrtnodata(): +def test_gdalalg_raster_mosaic_srcnodata_dstnodata(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src1_ds.GetRasterBand(1).Fill(1) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds]) alg.GetArg("output").Set("") alg.GetArg("srcnodata").Set([1]) - alg.GetArg("vrtnodata").Set([2]) + alg.GetArg("dstnodata").Set([2]) assert alg.Run() ds = alg.GetArg("output").Get().GetDataset() assert ds.GetRasterBand(1).Checksum() == 2 assert ds.GetRasterBand(1).GetNoDataValue() == 2 -def test_gdalalg_raster_buildvrt_hidenodata(): +def test_gdalalg_raster_mosaic_hidenodata(): src1_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) src1_ds.SetGeoTransform([2, 1, 0, 49, 0, -1]) src1_ds.GetRasterBand(1).Fill(1) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([src1_ds]) alg.GetArg("output").Set("") alg.GetArg("srcnodata").Set([1]) - alg.GetArg("vrtnodata").Set([2]) + alg.GetArg("dstnodata").Set([2]) alg.GetArg("hidenodata").Set(True) assert alg.Run() ds = alg.GetArg("output").Get().GetDataset() @@ -279,9 +279,9 @@ def test_gdalalg_raster_buildvrt_hidenodata(): assert ds.GetRasterBand(1).GetNoDataValue() is None -def test_gdalalg_raster_buildvrt_addalpha(): +def test_gdalalg_raster_mosaic_addalpha(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["../gcore/data/byte.tif"]) alg.GetArg("output").Set("") alg.GetArg("addalpha").Set(True) @@ -293,9 +293,9 @@ def test_gdalalg_raster_buildvrt_addalpha(): assert ds.GetRasterBand(2).Checksum() == 4873 -def test_gdalalg_raster_buildvrt_band(): +def test_gdalalg_raster_mosaic_band(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["../gcore/data/rgbsmall.tif"]) alg.GetArg("output").Set("") alg.GetArg("band").Set([3, 2]) @@ -306,9 +306,9 @@ def test_gdalalg_raster_buildvrt_band(): assert ds.GetRasterBand(2).Checksum() == 21053 -def test_gdalalg_raster_buildvrt_glob(): +def test_gdalalg_raster_mosaic_glob(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["../gcore/data/rgbsm?ll.tif"]) alg.GetArg("output").Set("") assert alg.Run() @@ -316,12 +316,12 @@ def test_gdalalg_raster_buildvrt_glob(): assert ds.RasterCount == 3 -def test_gdalalg_raster_buildvrt_at_filename(tmp_vsimem): +def test_gdalalg_raster_mosaic_at_filename(tmp_vsimem): input_file_list = str(tmp_vsimem / "tmp.txt") gdal.FileFromMemBuffer(input_file_list, "../gcore/data/byte.tif") - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set([f"@{input_file_list}"]) alg.GetArg("output").Set("") assert alg.Run() @@ -329,35 +329,73 @@ def test_gdalalg_raster_buildvrt_at_filename(tmp_vsimem): assert ds.GetRasterBand(1).Checksum() == 4672 -def test_gdalalg_raster_buildvrt_at_filename_error(): +def test_gdalalg_raster_mosaic_at_filename_error(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["@i_do_not_exist"]) alg.GetArg("output").Set("") - with pytest.raises(Exception, match="buildvrt: Cannot open i_do_not_exist"): + with pytest.raises(Exception, match="mosaic: Cannot open i_do_not_exist"): alg.Run() -def test_gdalalg_raster_buildvrt_output_ds_alread_set(): +def test_gdalalg_raster_mosaic_output_ds_alread_set(): out_ds = gdal.GetDriverByName("MEM").Create("", 1, 1) - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["../gcore/data/byte.tif"]) alg.GetArg("output").Set(out_ds) with pytest.raises( Exception, - match="buildvrt: gdal raster buildvrt does not support outputting to an already opened output dataset", + match="mosaic: gdal raster mosaic does not support outputting to an already opened output dataset", ): alg.Run() -def test_gdalalg_raster_buildvrt_co(): +def test_gdalalg_raster_mosaic_co(): - alg = get_buildvrt_alg() + alg = get_mosaic_alg() alg.GetArg("input").Set(["../gcore/data/byte.tif"]) alg.GetArg("output").Set("") alg.GetArg("creation-option").Set(["BLOCKXSIZE=10", "BLOCKYSIZE=15"]) assert alg.Run() ds = alg.GetArg("output").Get().GetDataset() assert ds.GetRasterBand(1).GetBlockSize() == [10, 15] + + +def test_gdalalg_raster_mosaic_tif_output_implicit(tmp_vsimem): + + out_filename = str(tmp_vsimem / "out.tif") + + alg = get_mosaic_alg() + assert alg.ParseRunAndFinalize(["../gcore/data/utmsmall.tif", out_filename]) + + with gdal.Open(out_filename) as ds: + assert ds.GetRasterBand(1).Checksum() == 50054 + + +def test_gdalalg_raster_mosaic_tif_output_explicit(tmp_vsimem): + + out_filename = str(tmp_vsimem / "out.xxx") + + alg = get_mosaic_alg() + assert alg.ParseRunAndFinalize( + ["--of=GTiff", "../gcore/data/utmsmall.tif", out_filename] + ) + + with gdal.Open(out_filename) as ds: + assert ds.GetRasterBand(1).Checksum() == 50054 + + +def test_gdalalg_raster_mosaic_tif_creation_options(tmp_vsimem): + + out_filename = str(tmp_vsimem / "out.xxx") + + alg = get_mosaic_alg() + assert alg.ParseRunAndFinalize( + ["--of=GTiff", "--co=TILED=YES", "../gcore/data/utmsmall.tif", out_filename] + ) + + with gdal.Open(out_filename) as ds: + assert ds.GetRasterBand(1).Checksum() == 50054 + assert ds.GetRasterBand(1).GetBlockSize() == [256, 256] diff --git a/doc/source/conf.py b/doc/source/conf.py index 51be7d65f45c..e6b42c6dc9e3 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -215,13 +215,6 @@ [author_evenr], 1, ), - ( - "programs/gdal_raster_buildvrt", - "gdal-raster-buildvrt", - "Build a virtual dataset (VRT)", - [author_evenr], - 1, - ), ( "programs/gdal_raster_info", "gdal-raster-info", @@ -250,6 +243,13 @@ [author_evenr], 1, ), + ( + "programs/gdal_raster_mosaic", + "gdal-raster-mosaic", + "Build a mosaic", + [author_evenr], + 1, + ), ( "programs/gdal_raster_overview_add", "gdal-raster-overview-add", diff --git a/doc/source/programs/gdal_raster.rst b/doc/source/programs/gdal_raster.rst index 1ffceb39b2b9..d3f90fa5c05c 100644 --- a/doc/source/programs/gdal_raster.rst +++ b/doc/source/programs/gdal_raster.rst @@ -19,10 +19,10 @@ Synopsis Usage: gdal raster where is one of: - - buildvrt: Build a virtual dataset (VRT). - 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. - overview: Manage overviews of a raster dataset. - pipeline: Process a raster dataset. - reproject: Reproject a raster dataset. @@ -30,10 +30,10 @@ Synopsis Available sub-commands ---------------------- -- :ref:`gdal_raster_buildvrt_subcommand` - :ref:`gdal_raster_info_subcommand` - :ref:`gdal_raster_clip_subcommand` - :ref:`gdal_raster_convert_subcommand` +- :ref:`gdal_raster_mosaic_subcommand` - :ref:`gdal_raster_overview_subcommand` - :ref:`gdal_raster_pipeline_subcommand` - :ref:`gdal_raster_reproject_subcommand` diff --git a/doc/source/programs/gdal_raster_buildvrt.rst b/doc/source/programs/gdal_raster_mosaic.rst similarity index 54% rename from doc/source/programs/gdal_raster_buildvrt.rst rename to doc/source/programs/gdal_raster_mosaic.rst index ce29983e665f..82c1ae28b689 100644 --- a/doc/source/programs/gdal_raster_buildvrt.rst +++ b/doc/source/programs/gdal_raster_mosaic.rst @@ -1,69 +1,72 @@ -.. _gdal_raster_buildvrt_subcommand: +.. _gdal_raster_mosaic_subcommand: ================================================================================ -"gdal raster buildvrt" sub-command +"gdal raster mosaic" sub-command ================================================================================ .. versionadded:: 3.11 .. only:: html - Build a virtual dataset (VRT). + Build a mosaic, either virtual (VRT) or materialized -.. Index:: gdal raster buildvrt +.. Index:: gdal raster mosaic Synopsis -------- .. code-block:: - Usage: gdal raster buildvrt [OPTIONS] + Usage: gdal raster mosaic [OPTIONS] - Build a virtual dataset (VRT). + Build a mosaic, 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] + -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 - --progress Display progress bar + -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: - --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 (in destination CRS units) - --target-aligned-pixels Round target extent to target resolution - --srcnodata Set nodata values for input bands. [1.. values] - --vrtnodata Set nodata values at the VRT band level. [1.. values] - --hidenodata Makes the VRT band not report the NoData. - --addalpha Adds an alpha mask band to the VRT when the source raster have none. + -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) + --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. + --addalpha Adds an alpha mask band to the destination when the source raster have none. + Description ----------- -This program builds a :ref:`VRT (Virtual Dataset) ` that is a mosaic of a list of -input GDAL datasets. The list of input GDAL datasets can be specified at the end +This program builds a mosaic of a list of input GDAL datasets, that can be +either a virtual mosaic in the :ref:`VRT (Virtual Dataset) ` format, +or in a more conventional raster format such as GeoTIFF. + +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. -With :option:`--separate`, each input goes into a separate band in the VRT dataset. Otherwise, -the files are considered as source rasters of a larger mosaic and the VRT file has the same number of +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. -If one GDAL dataset is made of several subdatasets and has 0 raster bands, -all the subdatasets will be added to the VRT rather than the dataset itself. - -:program:`gdal raster buildvrt` does some checks to ensure that all files that will be put -in the resulting VRT have similar characteristics: number of bands, projection, color +: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. (This is only true in the default mode, and not when using the :option:`--separate` option) @@ -77,10 +80,16 @@ changed in later versions. 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 VRT. + 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 @@ -90,7 +99,7 @@ The following options are available: required that all bands have the same datatype. All bands of each input file are added as separate - VRT bands, unless :option:`-b` is specified to select a subset of them. + output bands, unless :option:`-b` is specified to select a subset of them. .. option:: --resolution {|highest|lowest|average} @@ -108,10 +117,10 @@ The following options are available: .. option:: --bbox ,,, - Set georeferenced extents of VRT file. The values must be expressed in georeferenced units. - If not specified, the extent of the VRT is the minimum bounding box of the set of source rasters. - Pixels within the extent of the VRT 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:`--vrtnodata` + 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 @@ -128,18 +137,18 @@ The following options are available: will be used (if they exist). The value set by this option is written in the NODATA element of each ``ComplexSource`` element. -.. option:: --vrtnodata [,]... +.. option:: --dstnodata [,]... - Set nodata values at the VRT band level (different values can be supplied for each band). If more + 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 (:example:`vrtnodata`). If the option is not specified, + as a single operating system argument (:example:`dstnodata`). 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:: --addalpha - Adds an alpha mask band to the VRT when the source raster have none. Mainly useful for RGB sources (or grey-level sources). + Adds an alpha mask band to the output when the source raster have none. Mainly useful for RGB sources (or grey-level sources). The alpha band is filled on-the-fly with the value 0 in areas without any source raster, and with value 255 in areas with source raster. The effect is that a RGBA viewer will render the areas without source rasters as transparent and areas with source rasters as opaque. @@ -147,7 +156,7 @@ The following options are available: .. option:: --hidenodata - Even if any band contains nodata value, giving this option makes the VRT band + 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 @@ -163,19 +172,19 @@ Examples .. code-block:: bash - gdal raster biuldvrt --separate red.tif green.tif blue.tif rgb.vrt + 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: vrtnodata + :id: dstnodata .. code-block:: bash - gdal raster buildvrt --hidenodata --vrtnodata=0,0,255 doq/*.tif doq_index.vrt + gdal raster mosaic --hidenodata --dstnodata=0,0,255 doq/*.tif doq_index.vrt .. below is an allow-list for spelling checker. .. spelling:word-list:: - buildvrt + mosaic diff --git a/doc/source/programs/index.rst b/doc/source/programs/index.rst index aebb8a2c34b6..a1fe6c7e103e 100644 --- a/doc/source/programs/index.rst +++ b/doc/source/programs/index.rst @@ -29,11 +29,11 @@ single :program:`gdal` program that accepts commands and subcommands. gdal_info gdal_convert gdal_raster - gdal_raster_buildvrt gdal_raster_info gdal_raster_clip gdal_raster_convert gdal_raster_edit + gdal_raster_mosaic gdal_raster_overview gdal_raster_overview_add gdal_raster_overview_delete @@ -51,11 +51,11 @@ single :program:`gdal` program that accepts commands and subcommands. - :ref:`gdal_info_command`: Get information on a dataset - :ref:`gdal_convert_command`: Convert a dataset - :ref:`gdal_raster_command`: Entry point for raster commands - - :ref:`gdal_raster_buildvrt_subcommand`: buildvrt: Build a virtual dataset (VRT) - :ref:`gdal_raster_info_subcommand`: Get information on a raster dataset - :ref:`gdal_raster_clip_subcommand`: Clip a raster dataset - :ref:`gdal_raster_convert_subcommand`: Convert a raster dataset - :ref:`gdal_raster_edit_subcommand`: Edit in place a raster dataset + - :ref:`gdal_raster_mosaic_subcommand`: Build a mosaic, either virtual (VRT) or materialized. - :ref:`gdal_raster_overview_subcommand`: Manage overviews of a raster dataset - :ref:`gdal_raster_overview_add_subcommand`: Add overviews to a raster dataset - :ref:`gdal_raster_overview_delete_subcommand`: Remove overviews of a raster dataset From cdb9f363b59327991ece59c1d89685586cbdbac9 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 24 Jan 2025 17:25:10 +0100 Subject: [PATCH 2/2] Add 'gdal raster stack', and remove --separate option from 'gdal raster mosaic' --- apps/CMakeLists.txt | 1 + apps/gdalalg_raster.cpp | 2 + apps/gdalalg_raster_stack.cpp | 314 ++++++++++++++++++ apps/gdalalg_raster_stack.h | 59 ++++ .../utilities/test_gdalalg_raster_mosaic.py | 26 -- doc/source/conf.py | 7 + doc/source/programs/gdal_raster.rst | 12 +- doc/source/programs/gdal_raster_mosaic.rst | 24 -- doc/source/programs/gdal_raster_stack.rst | 143 ++++++++ doc/source/programs/index.rst | 2 + 10 files changed, 536 insertions(+), 54 deletions(-) create mode 100644 apps/gdalalg_raster_stack.cpp create mode 100644 apps/gdalalg_raster_stack.h create mode 100644 doc/source/programs/gdal_raster_stack.rst 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..4cb66f440c20 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,14 +60,9 @@ 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. -(This is only true in the default mode, and not when using the :option:`--separate` option) If the inputs spatially overlap, the order of the input list is used to determine priority. Files that are listed at the end are the ones @@ -92,15 +86,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 @@ -152,7 +137,6 @@ The following options are available: The alpha band is filled on-the-fly with the value 0 in areas without any source raster, and with value 255 in areas with source raster. The effect is that a RGBA viewer will render the areas without source rasters as transparent and areas with source rasters as opaque. - This option is not compatible with :option:`--separate`. .. option:: --hidenodata @@ -166,14 +150,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