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