Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster stack" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_stack.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "cpl_vsi_virtual.h"
17 :
18 : #include "gdal_priv.h"
19 : #include "gdal_utils.h"
20 :
21 : //! @cond Doxygen_Suppress
22 :
23 : #ifndef _
24 : #define _(x) (x)
25 : #endif
26 :
27 : /************************************************************************/
28 : /* GDALRasterStackAlgorithm::GDALRasterStackAlgorithm() */
29 : /************************************************************************/
30 :
31 5 : GDALRasterStackAlgorithm::GDALRasterStackAlgorithm()
32 5 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
33 : {
34 5 : m_supportsStreamedOutput = true;
35 :
36 5 : AddProgressArg();
37 : AddOutputFormatArg(&m_format, /* bStreamAllowed = */ true,
38 5 : /* bGDALGAllowed = */ true);
39 : AddArg(GDAL_ARG_NAME_INPUT, 'i',
40 : _("Input raster datasets (or specify a @<filename> to point to a "
41 : "file containing filenames)"),
42 10 : &m_inputDatasets, GDAL_OF_RASTER)
43 5 : .SetPositional()
44 5 : .SetMinCount(1)
45 5 : .SetAutoOpenDataset(false)
46 5 : .SetMetaVar("INPUTS");
47 5 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER);
48 5 : AddCreationOptionsArg(&m_creationOptions);
49 5 : AddBandArg(&m_bands);
50 5 : AddOverwriteArg(&m_overwrite);
51 : {
52 : auto &arg =
53 : AddArg("resolution", 0,
54 : _("Target resolution (in destination CRS units)"),
55 10 : &m_resolution)
56 5 : .SetDefault("same")
57 5 : .SetMetaVar("<xres>,<yres>|same|average|common|highest|lowest");
58 : arg.AddValidationAction(
59 1 : [this, &arg]()
60 : {
61 2 : const std::string val = arg.Get<std::string>();
62 2 : if (val != "average" && val != "highest" && val != "lowest" &&
63 2 : val != "same" && val != "common")
64 : {
65 : const auto aosTokens =
66 0 : CPLStringList(CSLTokenizeString2(val.c_str(), ",", 0));
67 0 : if (aosTokens.size() != 2 ||
68 0 : CPLGetValueType(aosTokens[0]) == CPL_VALUE_STRING ||
69 0 : CPLGetValueType(aosTokens[1]) == CPL_VALUE_STRING ||
70 0 : CPLAtof(aosTokens[0]) <= 0 ||
71 0 : CPLAtof(aosTokens[1]) <= 0)
72 : {
73 0 : ReportError(
74 : CE_Failure, CPLE_AppDefined,
75 : "resolution: two comma-separated positive "
76 : "values should be provided, or 'same', "
77 : "'average', 'common', 'highest' or 'lowest'");
78 0 : return false;
79 : }
80 : }
81 1 : return true;
82 5 : });
83 : }
84 : AddBBOXArg(&m_bbox, _("Target bounding box as xmin,ymin,xmax,ymax (in "
85 5 : "destination CRS units)"));
86 : AddArg("target-aligned-pixels", 0,
87 : _("Round target extent to target resolution"),
88 10 : &m_targetAlignedPixels)
89 5 : .AddHiddenAlias("tap");
90 : AddArg("src-nodata", 0, _("Set nodata values for input bands."),
91 10 : &m_srcNoData)
92 5 : .SetMinCount(1)
93 5 : .SetRepeatedArgAllowed(false);
94 : AddArg("dst-nodata", 0,
95 10 : _("Set nodata values at the destination band level."), &m_dstNoData)
96 5 : .SetMinCount(1)
97 5 : .SetRepeatedArgAllowed(false);
98 : AddArg("hide-nodata", 0,
99 : _("Makes the destination band not report the NoData."),
100 5 : &m_hideNoData);
101 5 : }
102 :
103 : /************************************************************************/
104 : /* GDALRasterStackAlgorithm::RunImpl() */
105 : /************************************************************************/
106 :
107 2 : bool GDALRasterStackAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
108 : void *pProgressData)
109 : {
110 2 : CPLAssert(!m_outputDataset.GetDatasetRef());
111 :
112 4 : std::vector<GDALDatasetH> ahInputDatasets;
113 4 : CPLStringList aosInputDatasetNames;
114 2 : bool foundByRef = false;
115 2 : bool foundByName = false;
116 5 : for (auto &ds : m_inputDatasets)
117 : {
118 3 : if (ds.GetDatasetRef())
119 : {
120 2 : foundByRef = true;
121 2 : ahInputDatasets.push_back(
122 2 : GDALDataset::ToHandle(ds.GetDatasetRef()));
123 : }
124 1 : else if (!ds.GetName().empty())
125 : {
126 1 : foundByName = true;
127 1 : if (ds.GetName()[0] == '@')
128 : {
129 : auto f = VSIVirtualHandleUniquePtr(
130 0 : VSIFOpenL(ds.GetName().c_str() + 1, "r"));
131 0 : if (!f)
132 : {
133 0 : ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s",
134 0 : ds.GetName().c_str() + 1);
135 0 : return false;
136 : }
137 0 : while (const char *filename = CPLReadLineL(f.get()))
138 : {
139 0 : aosInputDatasetNames.push_back(filename);
140 0 : }
141 : }
142 1 : else if (ds.GetName().find_first_of("*?[") != std::string::npos)
143 : {
144 0 : CPLStringList aosMatches(VSIGlob(ds.GetName().c_str(), nullptr,
145 0 : pfnProgress, pProgressData));
146 0 : for (const char *pszStr : aosMatches)
147 : {
148 0 : aosInputDatasetNames.push_back(pszStr);
149 : }
150 : }
151 : else
152 : {
153 2 : std::string osDatasetName = ds.GetName();
154 1 : if (!GetReferencePathForRelativePaths().empty())
155 : {
156 0 : osDatasetName = GDALDataset::BuildFilename(
157 : osDatasetName.c_str(),
158 0 : GetReferencePathForRelativePaths().c_str(), true);
159 : }
160 1 : aosInputDatasetNames.push_back(osDatasetName.c_str());
161 : }
162 : }
163 : }
164 2 : if (foundByName && foundByRef)
165 : {
166 0 : ReportError(CE_Failure, CPLE_NotSupported,
167 : "Input datasets should be provided either all by reference "
168 : "or all by name");
169 0 : return false;
170 : }
171 :
172 : const bool bVRTOutput =
173 3 : m_outputDataset.GetName().empty() || EQUAL(m_format.c_str(), "VRT") ||
174 5 : EQUAL(m_format.c_str(), "stream") ||
175 2 : EQUAL(CPLGetExtensionSafe(m_outputDataset.GetName().c_str()).c_str(),
176 : "VRT");
177 :
178 4 : CPLStringList aosOptions;
179 :
180 2 : aosOptions.push_back("-strict");
181 :
182 2 : aosOptions.push_back("-program_name");
183 2 : aosOptions.push_back("gdal raster stack");
184 :
185 2 : aosOptions.push_back("-separate");
186 :
187 : const auto aosTokens =
188 4 : CPLStringList(CSLTokenizeString2(m_resolution.c_str(), ",", 0));
189 2 : if (aosTokens.size() == 2)
190 : {
191 0 : aosOptions.push_back("-tr");
192 0 : aosOptions.push_back(aosTokens[0]);
193 0 : aosOptions.push_back(aosTokens[1]);
194 : }
195 : else
196 : {
197 2 : aosOptions.push_back("-resolution");
198 2 : aosOptions.push_back(m_resolution);
199 : }
200 :
201 2 : if (!m_bbox.empty())
202 : {
203 0 : aosOptions.push_back("-te");
204 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[0]));
205 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[1]));
206 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[2]));
207 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[3]));
208 : }
209 2 : if (m_targetAlignedPixels)
210 : {
211 0 : aosOptions.push_back("-tap");
212 : }
213 2 : if (!m_srcNoData.empty())
214 : {
215 0 : aosOptions.push_back("-srcnodata");
216 0 : std::string s;
217 0 : for (double v : m_srcNoData)
218 : {
219 0 : if (!s.empty())
220 0 : s += " ";
221 0 : s += CPLSPrintf("%.17g", v);
222 : }
223 0 : aosOptions.push_back(s);
224 : }
225 2 : if (!m_dstNoData.empty())
226 : {
227 0 : aosOptions.push_back("-vrtnodata");
228 0 : std::string s;
229 0 : for (double v : m_dstNoData)
230 : {
231 0 : if (!s.empty())
232 0 : s += " ";
233 0 : s += CPLSPrintf("%.17g", v);
234 : }
235 0 : aosOptions.push_back(s);
236 : }
237 2 : if (bVRTOutput)
238 : {
239 2 : for (const auto &co : m_creationOptions)
240 : {
241 0 : aosOptions.push_back("-co");
242 0 : aosOptions.push_back(co);
243 : }
244 : }
245 2 : for (const int b : m_bands)
246 : {
247 0 : aosOptions.push_back("-b");
248 0 : aosOptions.push_back(CPLSPrintf("%d", b));
249 : }
250 2 : if (m_hideNoData)
251 : {
252 0 : aosOptions.push_back("-hidenodata");
253 : }
254 :
255 : GDALBuildVRTOptions *psOptions =
256 2 : GDALBuildVRTOptionsNew(aosOptions.List(), nullptr);
257 2 : if (bVRTOutput)
258 : {
259 2 : GDALBuildVRTOptionsSetProgress(psOptions, pfnProgress, pProgressData);
260 : }
261 :
262 : auto poOutDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
263 2 : GDALBuildVRT(EQUAL(m_format.c_str(), "stream") ? ""
264 1 : : bVRTOutput ? m_outputDataset.GetName().c_str()
265 : : "",
266 1 : foundByName ? aosInputDatasetNames.size()
267 1 : : static_cast<int>(m_inputDatasets.size()),
268 3 : ahInputDatasets.empty() ? nullptr : ahInputDatasets.data(),
269 9 : aosInputDatasetNames.List(), psOptions, nullptr)));
270 2 : GDALBuildVRTOptionsFree(psOptions);
271 2 : bool bOK = poOutDS != nullptr;
272 2 : if (bOK)
273 : {
274 2 : if (bVRTOutput)
275 : {
276 2 : m_outputDataset.Set(std::move(poOutDS));
277 : }
278 : else
279 : {
280 0 : CPLStringList aosTranslateOptions;
281 0 : if (!m_format.empty())
282 : {
283 0 : aosTranslateOptions.AddString("-of");
284 0 : aosTranslateOptions.AddString(m_format.c_str());
285 : }
286 0 : for (const auto &co : m_creationOptions)
287 : {
288 0 : aosTranslateOptions.AddString("-co");
289 0 : aosTranslateOptions.AddString(co.c_str());
290 : }
291 :
292 : GDALTranslateOptions *psTranslateOptions =
293 0 : GDALTranslateOptionsNew(aosTranslateOptions.List(), nullptr);
294 0 : GDALTranslateOptionsSetProgress(psTranslateOptions, pfnProgress,
295 : pProgressData);
296 :
297 : auto poFinalDS =
298 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
299 0 : GDALTranslate(m_outputDataset.GetName().c_str(),
300 : GDALDataset::ToHandle(poOutDS.get()),
301 0 : psTranslateOptions, nullptr)));
302 0 : GDALTranslateOptionsFree(psTranslateOptions);
303 :
304 0 : bOK = poFinalDS != nullptr;
305 0 : if (bOK)
306 : {
307 0 : m_outputDataset.Set(std::move(poFinalDS));
308 : }
309 : }
310 : }
311 :
312 2 : return bOK;
313 : }
314 :
315 : //! @endcond
|
