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 : AddArg("band", 'b', _("Specify input band(s) number."), &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("srcnodata", 0, _("Set nodata values for input bands."),
91 10 : &m_srcNoData)
92 5 : .SetMinCount(1)
93 5 : .SetRepeatedArgAllowed(false);
94 : AddArg("dstnodata", 0,
95 10 : _("Set nodata values at the destination band level."), &m_dstNoData)
96 5 : .SetMinCount(1)
97 5 : .SetRepeatedArgAllowed(false);
98 : AddArg("hidenodata", 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 : if (m_outputDataset.GetDatasetRef())
111 : {
112 0 : ReportError(CE_Failure, CPLE_NotSupported,
113 : "gdal raster stack does not support outputting to an "
114 : "already opened output dataset");
115 0 : return false;
116 : }
117 :
118 4 : std::vector<GDALDatasetH> ahInputDatasets;
119 4 : CPLStringList aosInputDatasetNames;
120 2 : bool foundByRef = false;
121 2 : bool foundByName = false;
122 5 : for (auto &ds : m_inputDatasets)
123 : {
124 3 : if (ds.GetDatasetRef())
125 : {
126 2 : foundByRef = true;
127 2 : ahInputDatasets.push_back(
128 2 : GDALDataset::ToHandle(ds.GetDatasetRef()));
129 : }
130 1 : else if (!ds.GetName().empty())
131 : {
132 1 : foundByName = true;
133 1 : if (ds.GetName()[0] == '@')
134 : {
135 : auto f = VSIVirtualHandleUniquePtr(
136 0 : VSIFOpenL(ds.GetName().c_str() + 1, "r"));
137 0 : if (!f)
138 : {
139 0 : ReportError(CE_Failure, CPLE_FileIO, "Cannot open %s",
140 0 : ds.GetName().c_str() + 1);
141 0 : return false;
142 : }
143 0 : while (const char *filename = CPLReadLineL(f.get()))
144 : {
145 0 : aosInputDatasetNames.push_back(filename);
146 0 : }
147 : }
148 1 : else if (ds.GetName().find_first_of("*?[") != std::string::npos)
149 : {
150 0 : CPLStringList aosMatches(VSIGlob(ds.GetName().c_str(), nullptr,
151 0 : pfnProgress, pProgressData));
152 0 : for (const char *pszStr : aosMatches)
153 : {
154 0 : aosInputDatasetNames.push_back(pszStr);
155 : }
156 : }
157 : else
158 : {
159 2 : std::string osDatasetName = ds.GetName();
160 1 : if (!GetReferencePathForRelativePaths().empty())
161 : {
162 0 : osDatasetName = GDALDataset::BuildFilename(
163 : osDatasetName.c_str(),
164 0 : GetReferencePathForRelativePaths().c_str(), true);
165 : }
166 1 : aosInputDatasetNames.push_back(osDatasetName.c_str());
167 : }
168 : }
169 : }
170 2 : if (foundByName && foundByRef)
171 : {
172 0 : ReportError(CE_Failure, CPLE_NotSupported,
173 : "Input datasets should be provided either all by reference "
174 : "or all by name");
175 0 : return false;
176 : }
177 :
178 : VSIStatBufL sStat;
179 3 : if (!m_overwrite && !m_outputDataset.GetName().empty() &&
180 2 : (VSIStatL(m_outputDataset.GetName().c_str(), &sStat) == 0 ||
181 3 : std::unique_ptr<GDALDataset>(
182 2 : GDALDataset::Open(m_outputDataset.GetName().c_str()))))
183 : {
184 0 : ReportError(CE_Failure, CPLE_AppDefined,
185 : "File '%s' already exists. Specify the --overwrite "
186 : "option to overwrite it.",
187 0 : m_outputDataset.GetName().c_str());
188 0 : return false;
189 : }
190 :
191 : const bool bVRTOutput =
192 3 : m_outputDataset.GetName().empty() || EQUAL(m_format.c_str(), "VRT") ||
193 5 : EQUAL(m_format.c_str(), "stream") ||
194 2 : EQUAL(CPLGetExtensionSafe(m_outputDataset.GetName().c_str()).c_str(),
195 : "VRT");
196 :
197 4 : CPLStringList aosOptions;
198 :
199 2 : aosOptions.push_back("-strict");
200 :
201 2 : aosOptions.push_back("-program_name");
202 2 : aosOptions.push_back("gdal raster stack");
203 :
204 2 : aosOptions.push_back("-separate");
205 :
206 : const auto aosTokens =
207 4 : CPLStringList(CSLTokenizeString2(m_resolution.c_str(), ",", 0));
208 2 : if (aosTokens.size() == 2)
209 : {
210 0 : aosOptions.push_back("-tr");
211 0 : aosOptions.push_back(aosTokens[0]);
212 0 : aosOptions.push_back(aosTokens[1]);
213 : }
214 : else
215 : {
216 2 : aosOptions.push_back("-resolution");
217 2 : aosOptions.push_back(m_resolution);
218 : }
219 :
220 2 : if (!m_bbox.empty())
221 : {
222 0 : aosOptions.push_back("-te");
223 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[0]));
224 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[1]));
225 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[2]));
226 0 : aosOptions.push_back(CPLSPrintf("%.17g", m_bbox[3]));
227 : }
228 2 : if (m_targetAlignedPixels)
229 : {
230 0 : aosOptions.push_back("-tap");
231 : }
232 2 : if (!m_srcNoData.empty())
233 : {
234 0 : aosOptions.push_back("-srcnodata");
235 0 : std::string s;
236 0 : for (double v : m_srcNoData)
237 : {
238 0 : if (!s.empty())
239 0 : s += " ";
240 0 : s += CPLSPrintf("%.17g", v);
241 : }
242 0 : aosOptions.push_back(s);
243 : }
244 2 : if (!m_dstNoData.empty())
245 : {
246 0 : aosOptions.push_back("-vrtnodata");
247 0 : std::string s;
248 0 : for (double v : m_dstNoData)
249 : {
250 0 : if (!s.empty())
251 0 : s += " ";
252 0 : s += CPLSPrintf("%.17g", v);
253 : }
254 0 : aosOptions.push_back(s);
255 : }
256 2 : if (bVRTOutput)
257 : {
258 2 : for (const auto &co : m_creationOptions)
259 : {
260 0 : aosOptions.push_back("-co");
261 0 : aosOptions.push_back(co);
262 : }
263 : }
264 2 : for (const int b : m_bands)
265 : {
266 0 : aosOptions.push_back("-b");
267 0 : aosOptions.push_back(CPLSPrintf("%d", b));
268 : }
269 2 : if (m_hideNoData)
270 : {
271 0 : aosOptions.push_back("-hidenodata");
272 : }
273 :
274 : GDALBuildVRTOptions *psOptions =
275 2 : GDALBuildVRTOptionsNew(aosOptions.List(), nullptr);
276 2 : if (bVRTOutput)
277 : {
278 2 : GDALBuildVRTOptionsSetProgress(psOptions, pfnProgress, pProgressData);
279 : }
280 :
281 : auto poOutDS = std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
282 2 : GDALBuildVRT(EQUAL(m_format.c_str(), "stream") ? ""
283 1 : : bVRTOutput ? m_outputDataset.GetName().c_str()
284 : : "",
285 1 : foundByName ? aosInputDatasetNames.size()
286 1 : : static_cast<int>(m_inputDatasets.size()),
287 3 : ahInputDatasets.empty() ? nullptr : ahInputDatasets.data(),
288 9 : aosInputDatasetNames.List(), psOptions, nullptr)));
289 2 : GDALBuildVRTOptionsFree(psOptions);
290 2 : bool bOK = poOutDS != nullptr;
291 2 : if (bOK)
292 : {
293 2 : if (bVRTOutput)
294 : {
295 2 : m_outputDataset.Set(std::move(poOutDS));
296 : }
297 : else
298 : {
299 0 : CPLStringList aosTranslateOptions;
300 0 : if (!m_format.empty())
301 : {
302 0 : aosTranslateOptions.AddString("-of");
303 0 : aosTranslateOptions.AddString(m_format.c_str());
304 : }
305 0 : for (const auto &co : m_creationOptions)
306 : {
307 0 : aosTranslateOptions.AddString("-co");
308 0 : aosTranslateOptions.AddString(co.c_str());
309 : }
310 :
311 : GDALTranslateOptions *psTranslateOptions =
312 0 : GDALTranslateOptionsNew(aosTranslateOptions.List(), nullptr);
313 0 : GDALTranslateOptionsSetProgress(psTranslateOptions, pfnProgress,
314 : pProgressData);
315 :
316 : auto poFinalDS =
317 : std::unique_ptr<GDALDataset>(GDALDataset::FromHandle(
318 0 : GDALTranslate(m_outputDataset.GetName().c_str(),
319 : GDALDataset::ToHandle(poOutDS.get()),
320 0 : psTranslateOptions, nullptr)));
321 0 : GDALTranslateOptionsFree(psTranslateOptions);
322 :
323 0 : bOK = poFinalDS != nullptr;
324 0 : if (bOK)
325 : {
326 0 : m_outputDataset.Set(std::move(poFinalDS));
327 : }
328 : }
329 : }
330 :
331 2 : return bOK;
332 : }
333 :
334 : //! @endcond
|
