Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster contour" subcommand
5 : * Author: Alessandro Pasotti <elpaso at itopen dot it>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Alessandro Pasotti <elpaso at itopen dot it>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include <cmath>
14 :
15 : #include "gdalalg_raster_contour.h"
16 :
17 : #include "cpl_conv.h"
18 : #include "gdal_priv.h"
19 : #include "gdal_utils.h"
20 : #include "gdal_alg.h"
21 : #include "gdal_utils_priv.h"
22 :
23 : //! @cond Doxygen_Suppress
24 :
25 : #ifndef _
26 : #define _(x) (x)
27 : #endif
28 :
29 : /************************************************************************/
30 : /* GDALRasterContourAlgorithm::GDALRasterContourAlgorithm() */
31 : /************************************************************************/
32 :
33 46 : GDALRasterContourAlgorithm::GDALRasterContourAlgorithm(bool standaloneStep)
34 : : GDALPipelineStepAlgorithm(
35 : NAME, DESCRIPTION, HELP_URL,
36 0 : ConstructorOptions()
37 46 : .SetStandaloneStep(standaloneStep)
38 46 : .SetAddAppendLayerArgument(false)
39 46 : .SetAddOverwriteLayerArgument(false)
40 46 : .SetAddUpdateArgument(false)
41 46 : .SetAddUpsertArgument(false)
42 46 : .SetAddSkipErrorsArgument(false)
43 92 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
44 : {
45 46 : m_outputLayerName = "contour";
46 :
47 46 : AddProgressArg();
48 46 : if (standaloneStep)
49 : {
50 27 : AddRasterInputArgs(false, false);
51 27 : AddVectorOutputArgs(false, false);
52 : }
53 : else
54 : {
55 19 : AddRasterHiddenInputDatasetArg();
56 19 : AddOutputLayerNameArg(/* hiddenForCLI = */ false,
57 : /* shortNameOutputLayerAllowed = */ false);
58 : }
59 :
60 : // gdal_contour specific options
61 46 : AddBandArg(&m_band).SetDefault(1);
62 :
63 : AddArg("elevation-name", 0, _("Name of the elevation field"),
64 46 : &m_elevAttributeName);
65 46 : AddArg("min-name", 0, _("Name of the minimum elevation field"), &m_amin);
66 46 : AddArg("max-name", 0, _("Name of the maximum elevation field"), &m_amax);
67 46 : AddArg("3d", 0, _("Force production of 3D vectors instead of 2D"), &m_3d);
68 :
69 : AddArg("src-nodata", 0, _("Input pixel value to treat as 'nodata'"),
70 46 : &m_sNodata);
71 92 : AddArg("interval", 0, _("Elevation interval between contours"), &m_interval)
72 92 : .SetMutualExclusionGroup("levels")
73 46 : .SetMinValueExcluded(0);
74 92 : AddArg("levels", 0, _("List of contour levels"), &m_levels)
75 46 : .SetMutualExclusionGroup("levels");
76 : AddArg("exp-base", 'e', _("Base for exponential contour level generation"),
77 92 : &m_expBase)
78 46 : .SetMutualExclusionGroup("levels");
79 92 : AddArg("offset", 0, _("Offset to apply to contour levels"), &m_offset)
80 46 : .AddAlias("off");
81 : AddArg("polygonize", 'p', _("Create polygons instead of lines"),
82 46 : &m_polygonize);
83 : AddArg("group-transactions", 0,
84 : _("Group n features per transaction (default 100 000)"),
85 92 : &m_groupTransactions)
86 46 : .SetMinValueIncluded(0);
87 46 : }
88 :
89 : /************************************************************************/
90 : /* GDALRasterContourAlgorithm::RunImpl() */
91 : /************************************************************************/
92 :
93 14 : bool GDALRasterContourAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
94 : void *pProgressData)
95 : {
96 14 : GDALPipelineStepRunContext stepCtxt;
97 14 : stepCtxt.m_pfnProgress = pfnProgress;
98 14 : stepCtxt.m_pProgressData = pProgressData;
99 14 : return RunPreStepPipelineValidations() && RunStep(stepCtxt);
100 : }
101 :
102 : /************************************************************************/
103 : /* GDALRasterContourAlgorithm::RunStep() */
104 : /************************************************************************/
105 :
106 16 : bool GDALRasterContourAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
107 : {
108 16 : CPLErrorReset();
109 :
110 16 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
111 16 : CPLAssert(poSrcDS);
112 :
113 16 : CPLAssert(!m_outputDataset.GetDatasetRef());
114 :
115 32 : CPLStringList aosOptions;
116 :
117 32 : std::string outputFilename;
118 16 : if (m_standaloneStep)
119 : {
120 14 : outputFilename = m_outputDataset.GetName();
121 14 : if (!m_format.empty())
122 : {
123 0 : aosOptions.AddString("-of");
124 0 : aosOptions.AddString(m_format.c_str());
125 : }
126 :
127 15 : for (const auto &co : m_creationOptions)
128 : {
129 1 : aosOptions.push_back("-co");
130 1 : aosOptions.push_back(co.c_str());
131 : }
132 :
133 15 : for (const auto &co : m_layerCreationOptions)
134 : {
135 1 : aosOptions.push_back("-lco");
136 1 : aosOptions.push_back(co.c_str());
137 : }
138 : }
139 : else
140 : {
141 2 : if (GetGDALDriverManager()->GetDriverByName("GPKG"))
142 : {
143 2 : aosOptions.AddString("-of");
144 2 : aosOptions.AddString("GPKG");
145 :
146 2 : outputFilename = CPLGenerateTempFilenameSafe("_contour") + ".gpkg";
147 : }
148 : else
149 : {
150 0 : aosOptions.AddString("-of");
151 0 : aosOptions.AddString("MEM");
152 : }
153 : }
154 :
155 16 : if (m_band > 0)
156 : {
157 16 : aosOptions.AddString("-b");
158 16 : aosOptions.AddString(CPLSPrintf("%d", m_band));
159 : }
160 16 : if (!m_elevAttributeName.empty())
161 : {
162 5 : aosOptions.AddString("-a");
163 5 : aosOptions.AddString(m_elevAttributeName);
164 : }
165 16 : if (!m_amin.empty())
166 : {
167 6 : aosOptions.AddString("-amin");
168 6 : aosOptions.AddString(m_amin);
169 : }
170 16 : if (!m_amax.empty())
171 : {
172 6 : aosOptions.AddString("-amax");
173 6 : aosOptions.AddString(m_amax);
174 : }
175 16 : if (m_3d)
176 : {
177 0 : aosOptions.AddString("-3d");
178 : }
179 16 : if (!std::isnan(m_sNodata))
180 : {
181 1 : aosOptions.AddString("-snodata");
182 1 : aosOptions.AddString(CPLSPrintf("%.16g", m_sNodata));
183 : }
184 16 : if (m_levels.size() > 0)
185 : {
186 12 : for (const auto &level : m_levels)
187 : {
188 8 : aosOptions.AddString("-fl");
189 8 : aosOptions.AddString(level);
190 : }
191 : }
192 16 : if (!std::isnan(m_interval))
193 : {
194 10 : aosOptions.AddString("-i");
195 10 : aosOptions.AddString(CPLSPrintf("%.16g", m_interval));
196 : }
197 16 : if (m_expBase > 0)
198 : {
199 1 : aosOptions.AddString("-e");
200 1 : aosOptions.AddString(CPLSPrintf("%d", m_expBase));
201 : }
202 16 : if (!std::isnan(m_offset))
203 : {
204 1 : aosOptions.AddString("-off");
205 1 : aosOptions.AddString(CPLSPrintf("%.16g", m_offset));
206 : }
207 16 : if (m_polygonize)
208 : {
209 5 : aosOptions.AddString("-p");
210 : }
211 16 : if (!m_outputLayerName.empty())
212 : {
213 16 : aosOptions.AddString("-nln");
214 16 : aosOptions.AddString(m_outputLayerName);
215 : }
216 :
217 : // Check that one of --interval, --levels, --exp-base is specified
218 16 : if (m_levels.size() == 0 && std::isnan(m_interval) && m_expBase == 0)
219 : {
220 1 : ReportError(
221 : CE_Failure, CPLE_AppDefined,
222 : "One of 'interval', 'levels', 'exp-base' must be specified.");
223 1 : return false;
224 : }
225 :
226 : // Check that interval is not negative
227 15 : if (!std::isnan(m_interval) && m_interval < 0)
228 : {
229 0 : ReportError(CE_Failure, CPLE_AppDefined,
230 : "Interval must be a positive number.");
231 0 : return false;
232 : }
233 :
234 15 : aosOptions.AddString(m_inputDataset[0].GetName());
235 15 : aosOptions.AddString(outputFilename);
236 :
237 15 : bool bRet = false;
238 30 : GDALContourOptionsForBinary optionsForBinary;
239 : std::unique_ptr<GDALContourOptions, decltype(&GDALContourOptionsFree)>
240 : psOptions{GDALContourOptionsNew(aosOptions.List(), &optionsForBinary),
241 15 : GDALContourOptionsFree};
242 15 : if (psOptions)
243 : {
244 15 : GDALDatasetH hSrcDS{poSrcDS};
245 15 : GDALRasterBandH hBand{nullptr};
246 15 : GDALDatasetH hDstDS{m_outputDataset.GetDatasetRef()};
247 15 : OGRLayerH hLayer{nullptr};
248 15 : char **papszStringOptions = nullptr;
249 :
250 15 : bRet = GDALContourProcessOptions(psOptions.get(), &papszStringOptions,
251 : &hSrcDS, &hBand, &hDstDS,
252 : &hLayer) == CE_None;
253 :
254 15 : if (bRet)
255 : {
256 15 : bRet = GDALContourGenerateEx(hBand, hLayer, papszStringOptions,
257 : ctxt.m_pfnProgress,
258 : ctxt.m_pProgressData) == CE_None;
259 : }
260 :
261 15 : CSLDestroy(papszStringOptions);
262 :
263 15 : auto poDstDS = GDALDataset::FromHandle(hDstDS);
264 15 : if (bRet)
265 : {
266 15 : bRet = poDstDS != nullptr;
267 : }
268 15 : if (poDstDS && !m_standaloneStep && !outputFilename.empty())
269 : {
270 2 : poDstDS->MarkSuppressOnClose();
271 2 : if (bRet)
272 2 : bRet = poDstDS->FlushCache() == CE_None;
273 : #if !defined(__APPLE__)
274 : // For some unknown reason, unlinking the file on MacOSX
275 : // leads to later "disk I/O error". See https://github.com/OSGeo/gdal/issues/13794
276 2 : VSIUnlink(outputFilename.c_str());
277 : #endif
278 : }
279 15 : m_outputDataset.Set(std::unique_ptr<GDALDataset>(poDstDS));
280 : }
281 :
282 15 : return bRet;
283 : }
284 :
285 : GDALRasterContourAlgorithmStandalone::~GDALRasterContourAlgorithmStandalone() =
286 : default;
287 :
288 : //! @endcond
|
