Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Contour Generator
4 : * Purpose: Contour Generator mainline.
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2003, Applied Coherent Technology (www.actgate.com).
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2018, Oslandia <infos at oslandia dot com>
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "cpl_conv.h"
32 : #include "cpl_string.h"
33 : #include "gdal_version.h"
34 : #include "gdal.h"
35 : #include "gdal_alg.h"
36 : #include "gdalargumentparser.h"
37 : #include "ogr_api.h"
38 : #include "ogr_srs_api.h"
39 : #include "commonutils.h"
40 :
41 : /************************************************************************/
42 : /* GDALContourOptions */
43 : /************************************************************************/
44 :
45 : struct GDALContourOptions
46 : {
47 : int nBand = 1;
48 : double dfInterval = 0.0;
49 : double dfNoData = 0.0;
50 : double dfOffset = 0.0;
51 : double dfExpBase = 0.0;
52 : bool b3D = false;
53 : bool bPolygonize = false;
54 : bool bNoDataSet = false;
55 : bool bIgnoreNoData = false;
56 : std::string osNewLayerName = "contour";
57 : std::string osFormat;
58 : std::string osElevAttrib;
59 : std::string osElevAttribMin;
60 : std::string osElevAttribMax;
61 : std::vector<double> adfFixedLevels;
62 : CPLStringList aosOpenOptions;
63 : CPLStringList aosCreationOptions;
64 : bool bQuiet = false;
65 : std::string aosDestFilename;
66 : std::string aosSrcFilename;
67 : };
68 :
69 : /************************************************************************/
70 : /* GDALContourAppOptionsGetParser() */
71 : /************************************************************************/
72 :
73 : static std::unique_ptr<GDALArgumentParser>
74 9 : GDALContourAppOptionsGetParser(GDALContourOptions *psOptions)
75 : {
76 : auto argParser = std::make_unique<GDALArgumentParser>(
77 9 : "gdal_contour", /* bForBinary */ true);
78 :
79 9 : argParser->add_description(_("Creates contour lines from a raster file."));
80 9 : argParser->add_epilog(_(
81 : "For more details, consult the full documentation for the gdal_contour "
82 9 : "utility: http://gdal.org/gdal_contour.html"));
83 :
84 9 : argParser->add_argument("-b")
85 18 : .metavar("<name>")
86 9 : .default_value(1)
87 9 : .nargs(1)
88 9 : .scan<'i', int>()
89 9 : .store_into(psOptions->nBand)
90 9 : .help(_("Select an input band band containing the DEM data."));
91 :
92 9 : argParser->add_argument("-a")
93 18 : .metavar("<name>")
94 9 : .store_into(psOptions->osElevAttrib)
95 : .help(_("Provides a name for the attribute in which to put the "
96 9 : "elevation."));
97 :
98 9 : argParser->add_argument("-amin")
99 18 : .metavar("<name>")
100 9 : .store_into(psOptions->osElevAttribMin)
101 : .help(_("Provides a name for the attribute in which to put the minimum "
102 9 : "elevation."));
103 :
104 9 : argParser->add_argument("-amax")
105 18 : .metavar("<name>")
106 9 : .store_into(psOptions->osElevAttribMax)
107 : .help(_("Provides a name for the attribute in which to put the maximum "
108 9 : "elevation."));
109 :
110 9 : argParser->add_argument("-3d")
111 9 : .flag()
112 9 : .store_into(psOptions->b3D)
113 9 : .help(_("Force production of 3D vectors instead of 2D."));
114 :
115 9 : argParser->add_argument("-inodata")
116 9 : .flag()
117 9 : .store_into(psOptions->bIgnoreNoData)
118 : .help(_("Ignore any nodata value implied in the dataset - treat all "
119 9 : "values as valid."));
120 :
121 9 : argParser->add_argument("-snodata")
122 18 : .metavar("<value>")
123 9 : .scan<'g', double>()
124 : .action(
125 0 : [psOptions](const auto &d)
126 : {
127 0 : psOptions->bNoDataSet = true;
128 0 : psOptions->dfNoData = CPLAtofM(d.c_str());
129 9 : })
130 9 : .help(_("Input pixel value to treat as \"nodata\"."));
131 :
132 9 : argParser->add_output_format_argument(psOptions->osFormat);
133 :
134 9 : argParser->add_argument("-dsco")
135 18 : .metavar("<NAME>=<VALUE>")
136 9 : .append()
137 0 : .action([psOptions](const std::string &s)
138 9 : { psOptions->aosOpenOptions.AddString(s.c_str()); })
139 9 : .help(_("Dataset creation option (format specific)."));
140 :
141 : argParser->add_layer_creation_options_argument(
142 9 : psOptions->aosCreationOptions);
143 :
144 9 : auto &group = argParser->add_mutually_exclusive_group();
145 :
146 9 : group.add_argument("-i")
147 18 : .metavar("<interval>")
148 9 : .scan<'g', double>()
149 9 : .store_into(psOptions->dfInterval)
150 9 : .help(_("Elevation interval between contours."));
151 :
152 9 : group.add_argument("-fl")
153 18 : .metavar("<level>")
154 9 : .nargs(argparse::nargs_pattern::at_least_one)
155 9 : .scan<'g', double>()
156 3 : .action([psOptions](const std::string &s)
157 12 : { psOptions->adfFixedLevels.push_back(CPLAtof(s.c_str())); })
158 9 : .help(_("Name one or more \"fixed levels\" to extract."));
159 :
160 9 : group.add_argument("-e")
161 18 : .metavar("<base>")
162 9 : .scan<'g', double>()
163 9 : .store_into(psOptions->dfExpBase)
164 : .help(_("Generate levels on an exponential scale: base ^ k, for k an "
165 9 : "integer."));
166 :
167 9 : argParser->add_argument("-off")
168 18 : .metavar("<offset>")
169 9 : .scan<'g', double>()
170 9 : .store_into(psOptions->dfOffset)
171 9 : .help(_("Offset from zero relative to which to interpret intervals."));
172 :
173 9 : argParser->add_argument("-nln")
174 18 : .metavar("<name>")
175 9 : .store_into(psOptions->osNewLayerName)
176 : .help(_("Provide a name for the output vector layer. Defaults to "
177 9 : "\"contour\"."));
178 :
179 9 : argParser->add_argument("-p")
180 9 : .flag()
181 9 : .store_into(psOptions->bPolygonize)
182 9 : .help(_("Generate contour polygons instead of lines."));
183 :
184 9 : argParser->add_quiet_argument(&psOptions->bQuiet);
185 :
186 9 : argParser->add_argument("src_filename")
187 9 : .store_into(psOptions->aosSrcFilename)
188 9 : .help("The source raster file.");
189 :
190 9 : argParser->add_argument("dst_filename")
191 9 : .store_into(psOptions->aosDestFilename)
192 9 : .help("The destination vector file.");
193 :
194 9 : return argParser;
195 : }
196 :
197 5 : static void CreateElevAttrib(const char *pszElevAttrib, OGRLayerH hLayer)
198 : {
199 5 : OGRFieldDefnH hFld = OGR_Fld_Create(pszElevAttrib, OFTReal);
200 5 : OGRErr eErr = OGR_L_CreateField(hLayer, hFld, FALSE);
201 5 : OGR_Fld_Destroy(hFld);
202 5 : if (eErr == OGRERR_FAILURE)
203 : {
204 0 : exit(1);
205 : }
206 5 : }
207 :
208 : /************************************************************************/
209 : /* main() */
210 : /************************************************************************/
211 :
212 9 : MAIN_START(argc, argv)
213 :
214 : {
215 :
216 9 : GDALProgressFunc pfnProgress = nullptr;
217 :
218 9 : EarlySetConfigOptions(argc, argv);
219 :
220 : /* -------------------------------------------------------------------- */
221 : /* Register standard GDAL drivers, and process generic GDAL */
222 : /* command options. */
223 : /* -------------------------------------------------------------------- */
224 :
225 9 : GDALAllRegister();
226 9 : OGRRegisterAll();
227 :
228 9 : argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
229 9 : if (argc < 1)
230 0 : exit(-argc);
231 :
232 : /* -------------------------------------------------------------------- */
233 : /* Parse arguments. */
234 : /* -------------------------------------------------------------------- */
235 :
236 9 : if (argc < 2)
237 : {
238 : try
239 : {
240 0 : GDALContourOptions sOptions;
241 0 : auto argParser = GDALContourAppOptionsGetParser(&sOptions);
242 0 : fprintf(stderr, "%s\n", argParser->usage().c_str());
243 : }
244 0 : catch (const std::exception &err)
245 : {
246 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
247 0 : err.what());
248 : }
249 0 : exit(1);
250 : }
251 :
252 17 : GDALContourOptions sOptions;
253 :
254 : try
255 : {
256 17 : auto argParser = GDALContourAppOptionsGetParser(&sOptions);
257 9 : argParser->parse_args_without_binary_name(argv + 1);
258 : }
259 0 : catch (const std::exception &error)
260 : {
261 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", error.what());
262 0 : exit(1);
263 : }
264 :
265 16 : if (sOptions.aosSrcFilename.find("/vsistdout/") != std::string::npos ||
266 8 : sOptions.aosDestFilename.find("/vsistdout/") != std::string::npos)
267 : {
268 0 : sOptions.bQuiet = true;
269 : }
270 :
271 8 : if (!sOptions.bQuiet)
272 8 : pfnProgress = GDALTermProgress;
273 :
274 : /* -------------------------------------------------------------------- */
275 : /* Open source raster file. */
276 : /* -------------------------------------------------------------------- */
277 : GDALDatasetH hSrcDS =
278 8 : GDALOpen(sOptions.aosSrcFilename.c_str(), GA_ReadOnly);
279 8 : if (hSrcDS == nullptr)
280 0 : exit(2);
281 :
282 8 : GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, sOptions.nBand);
283 8 : if (hBand == nullptr)
284 : {
285 0 : CPLError(CE_Failure, CPLE_AppDefined,
286 : "Band %d does not exist on dataset.", sOptions.nBand);
287 0 : exit(2);
288 : }
289 :
290 8 : if (!sOptions.bNoDataSet && !sOptions.bIgnoreNoData)
291 : {
292 : int bNoDataSet;
293 8 : sOptions.dfNoData = GDALGetRasterNoDataValue(hBand, &bNoDataSet);
294 8 : sOptions.bNoDataSet = bNoDataSet;
295 : }
296 :
297 : /* -------------------------------------------------------------------- */
298 : /* Try to get a coordinate system from the raster. */
299 : /* -------------------------------------------------------------------- */
300 8 : OGRSpatialReferenceH hSRS = GDALGetSpatialRef(hSrcDS);
301 :
302 : /* -------------------------------------------------------------------- */
303 : /* Create the output file. */
304 : /* -------------------------------------------------------------------- */
305 8 : CPLString osFormat;
306 8 : if (sOptions.osFormat.empty())
307 : {
308 : const auto aoDrivers = GetOutputDriversFor(
309 8 : sOptions.aosDestFilename.c_str(), GDAL_OF_VECTOR);
310 8 : if (aoDrivers.empty())
311 : {
312 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot guess driver for %s",
313 : sOptions.aosDestFilename.c_str());
314 0 : exit(10);
315 : }
316 : else
317 : {
318 8 : if (aoDrivers.size() > 1)
319 : {
320 0 : CPLError(CE_Warning, CPLE_AppDefined,
321 : "Several drivers matching %s extension. Using %s",
322 : CPLGetExtension(sOptions.aosDestFilename.c_str()),
323 0 : aoDrivers[0].c_str());
324 : }
325 8 : osFormat = aoDrivers[0];
326 : }
327 : }
328 : else
329 : {
330 0 : osFormat = sOptions.osFormat;
331 : }
332 :
333 8 : OGRSFDriverH hDriver = OGRGetDriverByName(osFormat.c_str());
334 :
335 8 : if (hDriver == nullptr)
336 : {
337 0 : fprintf(stderr, "Unable to find format driver named %s.\n",
338 : osFormat.c_str());
339 0 : exit(10);
340 : }
341 :
342 8 : OGRDataSourceH hDS = OGR_Dr_CreateDataSource(
343 : hDriver, sOptions.aosDestFilename.c_str(), sOptions.aosCreationOptions);
344 8 : if (hDS == nullptr)
345 0 : exit(1);
346 :
347 16 : OGRLayerH hLayer = OGR_DS_CreateLayer(
348 : hDS, sOptions.osNewLayerName.c_str(), hSRS,
349 8 : sOptions.bPolygonize
350 0 : ? (sOptions.b3D ? wkbMultiPolygon25D : wkbMultiPolygon)
351 8 : : (sOptions.b3D ? wkbLineString25D : wkbLineString),
352 : sOptions.aosCreationOptions);
353 8 : if (hLayer == nullptr)
354 0 : exit(1);
355 :
356 8 : OGRFieldDefnH hFld = OGR_Fld_Create("ID", OFTInteger);
357 8 : OGR_Fld_SetWidth(hFld, 8);
358 8 : OGR_L_CreateField(hLayer, hFld, FALSE);
359 8 : OGR_Fld_Destroy(hFld);
360 :
361 8 : if (sOptions.bPolygonize)
362 : {
363 0 : if (!sOptions.osElevAttrib.empty())
364 : {
365 0 : sOptions.osElevAttrib.clear();
366 0 : CPLError(CE_Warning, CPLE_NotSupported,
367 : "-a is ignored in polygonal contouring mode. "
368 : "Use -amin and/or -amax instead");
369 : }
370 : }
371 : else
372 : {
373 16 : if (!sOptions.osElevAttribMin.empty() ||
374 8 : !sOptions.osElevAttribMax.empty())
375 : {
376 0 : sOptions.osElevAttribMin.clear();
377 0 : sOptions.osElevAttribMax.clear();
378 0 : CPLError(CE_Warning, CPLE_NotSupported,
379 : "-amin and/or -amax are ignored in line contouring mode. "
380 : "Use -a instead");
381 : }
382 : }
383 :
384 8 : if (!sOptions.osElevAttrib.empty())
385 : {
386 5 : CreateElevAttrib(sOptions.osElevAttrib.c_str(), hLayer);
387 : }
388 :
389 8 : if (!sOptions.osElevAttribMin.empty())
390 : {
391 0 : CreateElevAttrib(sOptions.osElevAttribMin.c_str(), hLayer);
392 : }
393 :
394 8 : if (!sOptions.osElevAttribMax.empty())
395 : {
396 0 : CreateElevAttrib(sOptions.osElevAttribMax.c_str(), hLayer);
397 : }
398 :
399 : /* -------------------------------------------------------------------- */
400 : /* Invoke. */
401 : /* -------------------------------------------------------------------- */
402 8 : int iIDField = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer), "ID");
403 8 : int iElevField = (sOptions.osElevAttrib.empty())
404 8 : ? -1
405 5 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
406 8 : sOptions.osElevAttrib.c_str());
407 :
408 : int iElevFieldMin =
409 8 : (sOptions.osElevAttribMin.empty())
410 8 : ? -1
411 0 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
412 8 : sOptions.osElevAttribMin.c_str());
413 :
414 : int iElevFieldMax =
415 8 : (sOptions.osElevAttribMax.empty())
416 8 : ? -1
417 0 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
418 8 : sOptions.osElevAttribMax.c_str());
419 :
420 8 : char **options = nullptr;
421 8 : if (!sOptions.adfFixedLevels.empty())
422 : {
423 1 : std::string values = "FIXED_LEVELS=";
424 4 : for (size_t i = 0; i < sOptions.adfFixedLevels.size(); i++)
425 : {
426 3 : const int sz = 32;
427 3 : char *newValue = new char[sz + 1];
428 3 : if (i == sOptions.adfFixedLevels.size() - 1)
429 : {
430 1 : CPLsnprintf(newValue, sz + 1, "%f", sOptions.adfFixedLevels[i]);
431 : }
432 : else
433 : {
434 2 : CPLsnprintf(newValue, sz + 1, "%f,",
435 2 : sOptions.adfFixedLevels[i]);
436 : }
437 3 : values = values + std::string(newValue);
438 3 : delete[] newValue;
439 : }
440 1 : options = CSLAddString(options, values.c_str());
441 : }
442 7 : else if (sOptions.dfExpBase != 0.0)
443 : {
444 : options =
445 0 : CSLAppendPrintf(options, "LEVEL_EXP_BASE=%f", sOptions.dfExpBase);
446 : }
447 7 : else if (sOptions.dfInterval != 0.0)
448 : {
449 : options =
450 7 : CSLAppendPrintf(options, "LEVEL_INTERVAL=%f", sOptions.dfInterval);
451 : }
452 :
453 8 : if (sOptions.dfOffset != 0.0)
454 : {
455 0 : options = CSLAppendPrintf(options, "LEVEL_BASE=%f", sOptions.dfOffset);
456 : }
457 :
458 8 : if (sOptions.bNoDataSet)
459 : {
460 5 : options = CSLAppendPrintf(options, "NODATA=%.19g", sOptions.dfNoData);
461 : }
462 8 : if (iIDField != -1)
463 : {
464 8 : options = CSLAppendPrintf(options, "ID_FIELD=%d", iIDField);
465 : }
466 8 : if (iElevField != -1)
467 : {
468 5 : options = CSLAppendPrintf(options, "ELEV_FIELD=%d", iElevField);
469 : }
470 8 : if (iElevFieldMin != -1)
471 : {
472 0 : options = CSLAppendPrintf(options, "ELEV_FIELD_MIN=%d", iElevFieldMin);
473 : }
474 8 : if (iElevFieldMax != -1)
475 : {
476 0 : options = CSLAppendPrintf(options, "ELEV_FIELD_MAX=%d", iElevFieldMax);
477 : }
478 8 : if (sOptions.bPolygonize)
479 : {
480 0 : options = CSLAppendPrintf(options, "POLYGONIZE=YES");
481 : }
482 :
483 : CPLErr eErr =
484 8 : GDALContourGenerateEx(hBand, hLayer, options, pfnProgress, nullptr);
485 :
486 8 : CSLDestroy(options);
487 8 : if (GDALClose(hDS) != CE_None)
488 0 : eErr = CE_Failure;
489 8 : GDALClose(hSrcDS);
490 :
491 8 : CSLDestroy(argv);
492 8 : GDALDestroyDriverManager();
493 8 : OGRCleanupAll();
494 :
495 8 : return (eErr == CE_None) ? 0 : 1;
496 : }
497 :
498 0 : MAIN_END
|
