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 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "cpl_conv.h"
16 : #include "cpl_string.h"
17 : #include "gdal_version.h"
18 : #include "gdal.h"
19 : #include "gdal_alg.h"
20 : #include "gdalargumentparser.h"
21 : #include "ogr_api.h"
22 : #include "ogr_srs_api.h"
23 : #include "commonutils.h"
24 :
25 : /************************************************************************/
26 : /* GDALContourOptions */
27 : /************************************************************************/
28 :
29 : struct GDALContourOptions
30 : {
31 : int nBand = 1;
32 : double dfInterval = 0.0;
33 : double dfNoData = 0.0;
34 : double dfOffset = 0.0;
35 : double dfExpBase = 0.0;
36 : bool b3D = false;
37 : bool bPolygonize = false;
38 : bool bNoDataSet = false;
39 : bool bIgnoreNoData = false;
40 : std::string osNewLayerName = "contour";
41 : std::string osFormat;
42 : std::string osElevAttrib;
43 : std::string osElevAttribMin;
44 : std::string osElevAttribMax;
45 : std::vector<std::string> aosFixedLevels;
46 : CPLStringList aosOpenOptions;
47 : CPLStringList aosCreationOptions;
48 : bool bQuiet = false;
49 : std::string aosDestFilename;
50 : std::string aosSrcFilename;
51 : GIntBig nGroupTransactions = 100 * 1000;
52 : };
53 :
54 : /************************************************************************/
55 : /* GDALContourAppOptionsGetParser() */
56 : /************************************************************************/
57 :
58 : static std::unique_ptr<GDALArgumentParser>
59 20 : GDALContourAppOptionsGetParser(GDALContourOptions *psOptions)
60 : {
61 : auto argParser = std::make_unique<GDALArgumentParser>(
62 20 : "gdal_contour", /* bForBinary */ true);
63 :
64 20 : argParser->add_description(_("Creates contour lines from a raster file."));
65 20 : argParser->add_epilog(_(
66 : "For more details, consult the full documentation for the gdal_contour "
67 20 : "utility: http://gdal.org/gdal_contour.html"));
68 :
69 20 : argParser->add_extra_usage_hint(
70 : _("One of -i, -fl or -e must be specified."));
71 :
72 20 : argParser->add_argument("-b")
73 40 : .metavar("<name>")
74 20 : .default_value(1)
75 20 : .nargs(1)
76 20 : .scan<'i', int>()
77 20 : .store_into(psOptions->nBand)
78 20 : .help(_("Select an input band band containing the DEM data."));
79 :
80 20 : argParser->add_argument("-a")
81 40 : .metavar("<name>")
82 20 : .store_into(psOptions->osElevAttrib)
83 : .help(_("Provides a name for the attribute in which to put the "
84 20 : "elevation."));
85 :
86 20 : argParser->add_argument("-amin")
87 40 : .metavar("<name>")
88 20 : .store_into(psOptions->osElevAttribMin)
89 : .help(_("Provides a name for the attribute in which to put the minimum "
90 20 : "elevation."));
91 :
92 20 : argParser->add_argument("-amax")
93 40 : .metavar("<name>")
94 20 : .store_into(psOptions->osElevAttribMax)
95 : .help(_("Provides a name for the attribute in which to put the maximum "
96 20 : "elevation."));
97 :
98 20 : argParser->add_argument("-3d")
99 20 : .flag()
100 20 : .store_into(psOptions->b3D)
101 20 : .help(_("Force production of 3D vectors instead of 2D."));
102 :
103 20 : argParser->add_argument("-inodata")
104 20 : .flag()
105 20 : .store_into(psOptions->bIgnoreNoData)
106 : .help(_("Ignore any nodata value implied in the dataset - treat all "
107 20 : "values as valid."));
108 :
109 20 : argParser->add_argument("-snodata")
110 40 : .metavar("<value>")
111 20 : .scan<'g', double>()
112 : .action(
113 0 : [psOptions](const auto &d)
114 : {
115 0 : psOptions->bNoDataSet = true;
116 0 : psOptions->dfNoData = CPLAtofM(d.c_str());
117 20 : })
118 20 : .help(_("Input pixel value to treat as \"nodata\"."));
119 :
120 20 : argParser->add_output_format_argument(psOptions->osFormat);
121 :
122 20 : argParser->add_dataset_creation_options_argument(psOptions->aosOpenOptions);
123 :
124 : argParser->add_layer_creation_options_argument(
125 20 : psOptions->aosCreationOptions);
126 :
127 20 : auto &group = argParser->add_mutually_exclusive_group();
128 :
129 20 : group.add_argument("-i")
130 40 : .metavar("<interval>")
131 20 : .scan<'g', double>()
132 20 : .store_into(psOptions->dfInterval)
133 20 : .help(_("Elevation interval between contours."));
134 :
135 20 : group.add_argument("-e")
136 40 : .metavar("<base>")
137 20 : .scan<'g', double>()
138 20 : .store_into(psOptions->dfExpBase)
139 : .help(_("Generate levels on an exponential scale: base ^ k, for k an "
140 20 : "integer."));
141 :
142 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
143 40 : argParser->add_argument("-fl").metavar("<level>").help(
144 20 : _("Name one or more \"fixed levels\" to extract."));
145 :
146 20 : argParser->add_argument("-off")
147 40 : .metavar("<offset>")
148 20 : .scan<'g', double>()
149 20 : .store_into(psOptions->dfOffset)
150 20 : .help(_("Offset from zero relative to which to interpret intervals."));
151 :
152 20 : argParser->add_argument("-nln")
153 40 : .metavar("<name>")
154 20 : .store_into(psOptions->osNewLayerName)
155 : .help(_("Provide a name for the output vector layer. Defaults to "
156 20 : "\"contour\"."));
157 :
158 20 : argParser->add_argument("-p")
159 20 : .flag()
160 20 : .store_into(psOptions->bPolygonize)
161 20 : .help(_("Generate contour polygons instead of lines."));
162 :
163 20 : argParser->add_argument("-gt")
164 40 : .metavar("<n>|unlimited")
165 : .action(
166 6 : [psOptions](const std::string &s)
167 : {
168 3 : if (EQUAL(s.c_str(), "unlimited"))
169 1 : psOptions->nGroupTransactions = -1;
170 : else
171 2 : psOptions->nGroupTransactions = atoi(s.c_str());
172 20 : })
173 20 : .help(_("Group <n> features per transaction."));
174 :
175 20 : argParser->add_quiet_argument(&psOptions->bQuiet);
176 :
177 20 : argParser->add_argument("src_filename")
178 20 : .store_into(psOptions->aosSrcFilename)
179 20 : .help("The source raster file.");
180 :
181 20 : argParser->add_argument("dst_filename")
182 20 : .store_into(psOptions->aosDestFilename)
183 20 : .help("The destination vector file.");
184 :
185 20 : return argParser;
186 : }
187 :
188 23 : static void CreateElevAttrib(const char *pszElevAttrib, OGRLayerH hLayer)
189 : {
190 23 : OGRFieldDefnH hFld = OGR_Fld_Create(pszElevAttrib, OFTReal);
191 23 : OGRErr eErr = OGR_L_CreateField(hLayer, hFld, FALSE);
192 23 : OGR_Fld_Destroy(hFld);
193 23 : if (eErr == OGRERR_FAILURE)
194 : {
195 0 : exit(1);
196 : }
197 23 : }
198 :
199 : /************************************************************************/
200 : /* main() */
201 : /************************************************************************/
202 :
203 21 : MAIN_START(argc, argv)
204 :
205 : {
206 :
207 21 : GDALProgressFunc pfnProgress = nullptr;
208 :
209 21 : EarlySetConfigOptions(argc, argv);
210 :
211 : /* -------------------------------------------------------------------- */
212 : /* Register standard GDAL drivers, and process generic GDAL */
213 : /* command options. */
214 : /* -------------------------------------------------------------------- */
215 :
216 21 : GDALAllRegister();
217 21 : OGRRegisterAll();
218 :
219 21 : argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
220 21 : if (argc < 1)
221 1 : exit(-argc);
222 :
223 : /* -------------------------------------------------------------------- */
224 : /* Parse arguments. */
225 : /* -------------------------------------------------------------------- */
226 :
227 20 : if (argc < 2)
228 : {
229 :
230 : try
231 : {
232 0 : GDALContourOptions sOptions;
233 0 : auto argParser = GDALContourAppOptionsGetParser(&sOptions);
234 0 : fprintf(stderr, "%s\n", argParser->usage().c_str());
235 : }
236 0 : catch (const std::exception &err)
237 : {
238 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
239 0 : err.what());
240 : }
241 0 : exit(1);
242 : }
243 :
244 39 : GDALContourOptions sOptions;
245 39 : CPLStringList aosArgv;
246 :
247 : try
248 : {
249 : /* -------------------------------------------------------------------- */
250 : /* Pre-processing for custom syntax that ArgumentParser does not */
251 : /* support. */
252 : /* -------------------------------------------------------------------- */
253 170 : for (int i = 1; i < argc && argv != nullptr && argv[i] != nullptr; i++)
254 : {
255 : // argparser is confused by arguments that have at_least_one
256 : // cardinality, if they immediately precede positional arguments.
257 150 : if (EQUAL(argv[i], "-fl") && argv[i + 1])
258 : {
259 11 : if (strchr(argv[i + 1], ' '))
260 : {
261 : const CPLStringList aosTokens(
262 0 : CSLTokenizeString(argv[i + 1]));
263 0 : for (const char *pszToken : aosTokens)
264 : {
265 : // Handle min/max special values
266 0 : if (EQUAL(pszToken, "MIN"))
267 : {
268 0 : sOptions.aosFixedLevels.push_back("MIN");
269 : }
270 0 : else if (EQUAL(pszToken, "MAX"))
271 : {
272 0 : sOptions.aosFixedLevels.push_back("MAX");
273 : }
274 : else
275 : {
276 0 : sOptions.aosFixedLevels.push_back(
277 0 : std::to_string(CPLAtof(pszToken)));
278 : }
279 : }
280 0 : i += 1;
281 : }
282 : else
283 : {
284 44 : auto isNumericOrMinMax = [](const char *pszArg) -> bool
285 : {
286 44 : if (EQUAL(pszArg, "MIN") || EQUAL(pszArg, "MAX"))
287 2 : return true;
288 42 : char *pszEnd = nullptr;
289 42 : CPLStrtod(pszArg, &pszEnd);
290 42 : return pszEnd != nullptr && pszEnd[0] == '\0';
291 : };
292 :
293 44 : while (i < argc - 1 && isNumericOrMinMax(argv[i + 1]))
294 : {
295 33 : if (EQUAL(argv[i + 1], "MIN"))
296 : {
297 1 : sOptions.aosFixedLevels.push_back("MIN");
298 : }
299 32 : else if (EQUAL(argv[i + 1], "MAX"))
300 : {
301 1 : sOptions.aosFixedLevels.push_back("MAX");
302 : }
303 : else
304 : {
305 31 : sOptions.aosFixedLevels.push_back(
306 62 : std::to_string(CPLAtof(argv[i + 1])));
307 : }
308 33 : i += 1;
309 : }
310 11 : }
311 : }
312 : else
313 : {
314 139 : aosArgv.AddString(argv[i]);
315 : }
316 : }
317 :
318 39 : auto argParser = GDALContourAppOptionsGetParser(&sOptions);
319 20 : argParser->parse_args_without_binary_name(aosArgv.List());
320 :
321 21 : if (sOptions.dfInterval == 0.0 && sOptions.aosFixedLevels.empty() &&
322 1 : sOptions.dfExpBase == 0.0)
323 : {
324 1 : fprintf(stderr, "%s\n", argParser->usage().c_str());
325 1 : exit(1);
326 : }
327 : }
328 0 : catch (const std::exception &error)
329 : {
330 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", error.what());
331 0 : exit(1);
332 : }
333 :
334 38 : if (sOptions.aosSrcFilename.find("/vsistdout/") != std::string::npos ||
335 19 : sOptions.aosDestFilename.find("/vsistdout/") != std::string::npos)
336 : {
337 0 : sOptions.bQuiet = true;
338 : }
339 :
340 19 : if (!sOptions.bQuiet)
341 19 : pfnProgress = GDALTermProgress;
342 :
343 : /* -------------------------------------------------------------------- */
344 : /* Open source raster file. */
345 : /* -------------------------------------------------------------------- */
346 : GDALDatasetH hSrcDS =
347 19 : GDALOpen(sOptions.aosSrcFilename.c_str(), GA_ReadOnly);
348 19 : if (hSrcDS == nullptr)
349 0 : exit(2);
350 :
351 19 : GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, sOptions.nBand);
352 19 : if (hBand == nullptr)
353 : {
354 0 : CPLError(CE_Failure, CPLE_AppDefined,
355 : "Band %d does not exist on dataset.", sOptions.nBand);
356 0 : exit(2);
357 : }
358 :
359 19 : if (!sOptions.bNoDataSet && !sOptions.bIgnoreNoData)
360 : {
361 : int bNoDataSet;
362 19 : sOptions.dfNoData = GDALGetRasterNoDataValue(hBand, &bNoDataSet);
363 19 : sOptions.bNoDataSet = bNoDataSet;
364 : }
365 :
366 : /* -------------------------------------------------------------------- */
367 : /* Try to get a coordinate system from the raster. */
368 : /* -------------------------------------------------------------------- */
369 19 : OGRSpatialReferenceH hSRS = GDALGetSpatialRef(hSrcDS);
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Create the output file. */
373 : /* -------------------------------------------------------------------- */
374 19 : CPLString osFormat;
375 19 : if (sOptions.osFormat.empty())
376 : {
377 : const auto aoDrivers = GetOutputDriversFor(
378 19 : sOptions.aosDestFilename.c_str(), GDAL_OF_VECTOR);
379 19 : if (aoDrivers.empty())
380 : {
381 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot guess driver for %s",
382 : sOptions.aosDestFilename.c_str());
383 0 : exit(10);
384 : }
385 : else
386 : {
387 19 : if (aoDrivers.size() > 1)
388 : {
389 0 : CPLError(CE_Warning, CPLE_AppDefined,
390 : "Several drivers matching %s extension. Using %s",
391 0 : CPLGetExtensionSafe(sOptions.aosDestFilename.c_str())
392 : .c_str(),
393 0 : aoDrivers[0].c_str());
394 : }
395 19 : osFormat = aoDrivers[0];
396 : }
397 : }
398 : else
399 : {
400 0 : osFormat = sOptions.osFormat;
401 : }
402 :
403 19 : OGRSFDriverH hDriver = OGRGetDriverByName(osFormat.c_str());
404 :
405 19 : if (hDriver == nullptr)
406 : {
407 0 : fprintf(stderr, "Unable to find format driver named %s.\n",
408 : osFormat.c_str());
409 0 : exit(10);
410 : }
411 :
412 19 : OGRDataSourceH hDS = OGR_Dr_CreateDataSource(
413 : hDriver, sOptions.aosDestFilename.c_str(), sOptions.aosCreationOptions);
414 19 : if (hDS == nullptr)
415 0 : exit(1);
416 :
417 38 : OGRLayerH hLayer = OGR_DS_CreateLayer(
418 : hDS, sOptions.osNewLayerName.c_str(), hSRS,
419 19 : sOptions.bPolygonize
420 7 : ? (sOptions.b3D ? wkbMultiPolygon25D : wkbMultiPolygon)
421 12 : : (sOptions.b3D ? wkbLineString25D : wkbLineString),
422 : sOptions.aosCreationOptions);
423 19 : if (hLayer == nullptr)
424 0 : exit(1);
425 :
426 19 : if (!OGR_L_TestCapability(hLayer, OLCTransactions))
427 : {
428 16 : sOptions.nGroupTransactions = 0;
429 : }
430 :
431 19 : OGRFieldDefnH hFld = OGR_Fld_Create("ID", OFTInteger);
432 19 : OGR_Fld_SetWidth(hFld, 8);
433 19 : OGR_L_CreateField(hLayer, hFld, FALSE);
434 19 : OGR_Fld_Destroy(hFld);
435 :
436 19 : if (sOptions.bPolygonize)
437 : {
438 7 : if (!sOptions.osElevAttrib.empty())
439 : {
440 0 : sOptions.osElevAttrib.clear();
441 0 : CPLError(CE_Warning, CPLE_NotSupported,
442 : "-a is ignored in polygonal contouring mode. "
443 : "Use -amin and/or -amax instead");
444 : }
445 : }
446 : else
447 : {
448 24 : if (!sOptions.osElevAttribMin.empty() ||
449 12 : !sOptions.osElevAttribMax.empty())
450 : {
451 0 : sOptions.osElevAttribMin.clear();
452 0 : sOptions.osElevAttribMax.clear();
453 0 : CPLError(CE_Warning, CPLE_NotSupported,
454 : "-amin and/or -amax are ignored in line contouring mode. "
455 : "Use -a instead");
456 : }
457 : }
458 :
459 19 : if (!sOptions.osElevAttrib.empty())
460 : {
461 9 : CreateElevAttrib(sOptions.osElevAttrib.c_str(), hLayer);
462 : }
463 :
464 19 : if (!sOptions.osElevAttribMin.empty())
465 : {
466 7 : CreateElevAttrib(sOptions.osElevAttribMin.c_str(), hLayer);
467 : }
468 :
469 19 : if (!sOptions.osElevAttribMax.empty())
470 : {
471 7 : CreateElevAttrib(sOptions.osElevAttribMax.c_str(), hLayer);
472 : }
473 :
474 : /* -------------------------------------------------------------------- */
475 : /* Invoke. */
476 : /* -------------------------------------------------------------------- */
477 19 : int iIDField = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer), "ID");
478 19 : int iElevField = (sOptions.osElevAttrib.empty())
479 19 : ? -1
480 9 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
481 19 : sOptions.osElevAttrib.c_str());
482 :
483 : int iElevFieldMin =
484 19 : (sOptions.osElevAttribMin.empty())
485 19 : ? -1
486 7 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
487 19 : sOptions.osElevAttribMin.c_str());
488 :
489 : int iElevFieldMax =
490 19 : (sOptions.osElevAttribMax.empty())
491 19 : ? -1
492 7 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
493 19 : sOptions.osElevAttribMax.c_str());
494 :
495 19 : char **options = nullptr;
496 19 : if (!sOptions.aosFixedLevels.empty())
497 : {
498 11 : std::string values = "FIXED_LEVELS=";
499 44 : for (size_t i = 0; i < sOptions.aosFixedLevels.size(); i++)
500 : {
501 33 : if (i == sOptions.aosFixedLevels.size() - 1)
502 : {
503 11 : values = values + sOptions.aosFixedLevels[i];
504 : }
505 : else
506 : {
507 22 : values = values + sOptions.aosFixedLevels[i] + ",";
508 : }
509 : }
510 11 : options = CSLAddString(options, values.c_str());
511 : }
512 :
513 19 : if (sOptions.dfExpBase != 0.0)
514 : {
515 : options =
516 5 : CSLAppendPrintf(options, "LEVEL_EXP_BASE=%f", sOptions.dfExpBase);
517 : }
518 14 : else if (sOptions.dfInterval != 0.0)
519 : {
520 : options =
521 12 : CSLAppendPrintf(options, "LEVEL_INTERVAL=%f", sOptions.dfInterval);
522 : }
523 :
524 19 : if (sOptions.dfOffset != 0.0)
525 : {
526 1 : options = CSLAppendPrintf(options, "LEVEL_BASE=%f", sOptions.dfOffset);
527 : }
528 :
529 19 : if (sOptions.bNoDataSet)
530 : {
531 10 : options = CSLAppendPrintf(options, "NODATA=%.19g", sOptions.dfNoData);
532 : }
533 19 : if (iIDField != -1)
534 : {
535 19 : options = CSLAppendPrintf(options, "ID_FIELD=%d", iIDField);
536 : }
537 19 : if (iElevField != -1)
538 : {
539 9 : options = CSLAppendPrintf(options, "ELEV_FIELD=%d", iElevField);
540 : }
541 19 : if (iElevFieldMin != -1)
542 : {
543 7 : options = CSLAppendPrintf(options, "ELEV_FIELD_MIN=%d", iElevFieldMin);
544 : }
545 19 : if (iElevFieldMax != -1)
546 : {
547 7 : options = CSLAppendPrintf(options, "ELEV_FIELD_MAX=%d", iElevFieldMax);
548 : }
549 19 : if (sOptions.bPolygonize)
550 : {
551 7 : options = CSLAppendPrintf(options, "POLYGONIZE=YES");
552 : }
553 19 : if (sOptions.nGroupTransactions)
554 : {
555 2 : options = CSLAppendPrintf(options, "COMMIT_INTERVAL=" CPL_FRMT_GIB,
556 : sOptions.nGroupTransactions);
557 : }
558 :
559 : CPLErr eErr =
560 19 : GDALContourGenerateEx(hBand, hLayer, options, pfnProgress, nullptr);
561 :
562 19 : CSLDestroy(options);
563 19 : if (GDALClose(hDS) != CE_None)
564 0 : eErr = CE_Failure;
565 19 : GDALClose(hSrcDS);
566 :
567 19 : CSLDestroy(argv);
568 19 : GDALDestroyDriverManager();
569 19 : OGRCleanupAll();
570 :
571 19 : return (eErr == CE_None) ? 0 : 1;
572 : }
573 :
574 0 : MAIN_END
|
