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