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