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