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<double> adfFixedLevels;
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").scan<'g', double>().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 21 : static void CreateElevAttrib(const char *pszElevAttrib, OGRLayerH hLayer)
189 : {
190 21 : OGRFieldDefnH hFld = OGR_Fld_Create(pszElevAttrib, OFTReal);
191 21 : OGRErr eErr = OGR_L_CreateField(hLayer, hFld, FALSE);
192 21 : OGR_Fld_Destroy(hFld);
193 21 : if (eErr == OGRERR_FAILURE)
194 : {
195 0 : exit(1);
196 : }
197 21 : }
198 :
199 : /************************************************************************/
200 : /* main() */
201 : /************************************************************************/
202 :
203 20 : MAIN_START(argc, argv)
204 :
205 : {
206 :
207 20 : GDALProgressFunc pfnProgress = nullptr;
208 :
209 20 : EarlySetConfigOptions(argc, argv);
210 :
211 : /* -------------------------------------------------------------------- */
212 : /* Register standard GDAL drivers, and process generic GDAL */
213 : /* command options. */
214 : /* -------------------------------------------------------------------- */
215 :
216 20 : GDALAllRegister();
217 20 : OGRRegisterAll();
218 :
219 20 : argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
220 20 : if (argc < 1)
221 0 : 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 38 : GDALContourOptions sOptions;
245 38 : CPLStringList aosArgv;
246 :
247 : try
248 : {
249 : /* -------------------------------------------------------------------- */
250 : /* Pre-processing for custom syntax that ArgumentParser does not */
251 : /* support. */
252 : /* -------------------------------------------------------------------- */
253 163 : 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 143 : if (EQUAL(argv[i], "-fl") && argv[i + 1])
258 : {
259 10 : 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 0 : sOptions.adfFixedLevels.push_back(CPLAtof(pszToken));
266 : }
267 0 : i += 1;
268 : }
269 : else
270 : {
271 38 : auto isNumeric = [](const char *pszArg) -> bool
272 : {
273 38 : char *pszEnd = nullptr;
274 38 : CPLStrtod(pszArg, &pszEnd);
275 38 : return pszEnd != nullptr && pszEnd[0] == '\0';
276 : };
277 :
278 38 : while (i < argc - 1 && isNumeric(argv[i + 1]))
279 : {
280 28 : sOptions.adfFixedLevels.push_back(CPLAtof(argv[i + 1]));
281 28 : i += 1;
282 : }
283 10 : }
284 : }
285 : else
286 : {
287 133 : aosArgv.AddString(argv[i]);
288 : }
289 : }
290 :
291 38 : auto argParser = GDALContourAppOptionsGetParser(&sOptions);
292 20 : argParser->parse_args_without_binary_name(aosArgv.List());
293 :
294 20 : if (sOptions.dfInterval == 0.0 && sOptions.adfFixedLevels.empty() &&
295 1 : sOptions.dfExpBase == 0.0)
296 : {
297 1 : fprintf(stderr, "%s\n", argParser->usage().c_str());
298 1 : exit(1);
299 : }
300 : }
301 0 : catch (const std::exception &error)
302 : {
303 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", error.what());
304 0 : exit(1);
305 : }
306 :
307 36 : if (sOptions.aosSrcFilename.find("/vsistdout/") != std::string::npos ||
308 18 : sOptions.aosDestFilename.find("/vsistdout/") != std::string::npos)
309 : {
310 0 : sOptions.bQuiet = true;
311 : }
312 :
313 18 : if (!sOptions.bQuiet)
314 18 : pfnProgress = GDALTermProgress;
315 :
316 : /* -------------------------------------------------------------------- */
317 : /* Open source raster file. */
318 : /* -------------------------------------------------------------------- */
319 : GDALDatasetH hSrcDS =
320 18 : GDALOpen(sOptions.aosSrcFilename.c_str(), GA_ReadOnly);
321 18 : if (hSrcDS == nullptr)
322 0 : exit(2);
323 :
324 18 : GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, sOptions.nBand);
325 18 : if (hBand == nullptr)
326 : {
327 0 : CPLError(CE_Failure, CPLE_AppDefined,
328 : "Band %d does not exist on dataset.", sOptions.nBand);
329 0 : exit(2);
330 : }
331 :
332 18 : if (!sOptions.bNoDataSet && !sOptions.bIgnoreNoData)
333 : {
334 : int bNoDataSet;
335 18 : sOptions.dfNoData = GDALGetRasterNoDataValue(hBand, &bNoDataSet);
336 18 : sOptions.bNoDataSet = bNoDataSet;
337 : }
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Try to get a coordinate system from the raster. */
341 : /* -------------------------------------------------------------------- */
342 18 : OGRSpatialReferenceH hSRS = GDALGetSpatialRef(hSrcDS);
343 :
344 : /* -------------------------------------------------------------------- */
345 : /* Create the output file. */
346 : /* -------------------------------------------------------------------- */
347 18 : CPLString osFormat;
348 18 : if (sOptions.osFormat.empty())
349 : {
350 : const auto aoDrivers = GetOutputDriversFor(
351 18 : sOptions.aosDestFilename.c_str(), GDAL_OF_VECTOR);
352 18 : if (aoDrivers.empty())
353 : {
354 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot guess driver for %s",
355 : sOptions.aosDestFilename.c_str());
356 0 : exit(10);
357 : }
358 : else
359 : {
360 18 : if (aoDrivers.size() > 1)
361 : {
362 0 : CPLError(CE_Warning, CPLE_AppDefined,
363 : "Several drivers matching %s extension. Using %s",
364 : CPLGetExtension(sOptions.aosDestFilename.c_str()),
365 0 : aoDrivers[0].c_str());
366 : }
367 18 : osFormat = aoDrivers[0];
368 : }
369 : }
370 : else
371 : {
372 0 : osFormat = sOptions.osFormat;
373 : }
374 :
375 18 : OGRSFDriverH hDriver = OGRGetDriverByName(osFormat.c_str());
376 :
377 18 : if (hDriver == nullptr)
378 : {
379 0 : fprintf(stderr, "Unable to find format driver named %s.\n",
380 : osFormat.c_str());
381 0 : exit(10);
382 : }
383 :
384 18 : OGRDataSourceH hDS = OGR_Dr_CreateDataSource(
385 : hDriver, sOptions.aosDestFilename.c_str(), sOptions.aosCreationOptions);
386 18 : if (hDS == nullptr)
387 0 : exit(1);
388 :
389 36 : OGRLayerH hLayer = OGR_DS_CreateLayer(
390 : hDS, sOptions.osNewLayerName.c_str(), hSRS,
391 18 : sOptions.bPolygonize
392 6 : ? (sOptions.b3D ? wkbMultiPolygon25D : wkbMultiPolygon)
393 12 : : (sOptions.b3D ? wkbLineString25D : wkbLineString),
394 : sOptions.aosCreationOptions);
395 18 : if (hLayer == nullptr)
396 0 : exit(1);
397 :
398 18 : if (!OGR_L_TestCapability(hLayer, OLCTransactions))
399 : {
400 15 : sOptions.nGroupTransactions = 0;
401 : }
402 :
403 18 : OGRFieldDefnH hFld = OGR_Fld_Create("ID", OFTInteger);
404 18 : OGR_Fld_SetWidth(hFld, 8);
405 18 : OGR_L_CreateField(hLayer, hFld, FALSE);
406 18 : OGR_Fld_Destroy(hFld);
407 :
408 18 : if (sOptions.bPolygonize)
409 : {
410 6 : if (!sOptions.osElevAttrib.empty())
411 : {
412 0 : sOptions.osElevAttrib.clear();
413 0 : CPLError(CE_Warning, CPLE_NotSupported,
414 : "-a is ignored in polygonal contouring mode. "
415 : "Use -amin and/or -amax instead");
416 : }
417 : }
418 : else
419 : {
420 24 : if (!sOptions.osElevAttribMin.empty() ||
421 12 : !sOptions.osElevAttribMax.empty())
422 : {
423 0 : sOptions.osElevAttribMin.clear();
424 0 : sOptions.osElevAttribMax.clear();
425 0 : CPLError(CE_Warning, CPLE_NotSupported,
426 : "-amin and/or -amax are ignored in line contouring mode. "
427 : "Use -a instead");
428 : }
429 : }
430 :
431 18 : if (!sOptions.osElevAttrib.empty())
432 : {
433 9 : CreateElevAttrib(sOptions.osElevAttrib.c_str(), hLayer);
434 : }
435 :
436 18 : if (!sOptions.osElevAttribMin.empty())
437 : {
438 6 : CreateElevAttrib(sOptions.osElevAttribMin.c_str(), hLayer);
439 : }
440 :
441 18 : if (!sOptions.osElevAttribMax.empty())
442 : {
443 6 : CreateElevAttrib(sOptions.osElevAttribMax.c_str(), hLayer);
444 : }
445 :
446 : /* -------------------------------------------------------------------- */
447 : /* Invoke. */
448 : /* -------------------------------------------------------------------- */
449 18 : int iIDField = OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer), "ID");
450 18 : int iElevField = (sOptions.osElevAttrib.empty())
451 18 : ? -1
452 9 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
453 18 : sOptions.osElevAttrib.c_str());
454 :
455 : int iElevFieldMin =
456 18 : (sOptions.osElevAttribMin.empty())
457 18 : ? -1
458 6 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
459 18 : sOptions.osElevAttribMin.c_str());
460 :
461 : int iElevFieldMax =
462 18 : (sOptions.osElevAttribMax.empty())
463 18 : ? -1
464 6 : : OGR_FD_GetFieldIndex(OGR_L_GetLayerDefn(hLayer),
465 18 : sOptions.osElevAttribMax.c_str());
466 :
467 18 : char **options = nullptr;
468 18 : if (!sOptions.adfFixedLevels.empty())
469 : {
470 10 : std::string values = "FIXED_LEVELS=";
471 38 : for (size_t i = 0; i < sOptions.adfFixedLevels.size(); i++)
472 : {
473 28 : const int sz = 32;
474 28 : char *newValue = new char[sz + 1];
475 28 : if (i == sOptions.adfFixedLevels.size() - 1)
476 : {
477 10 : CPLsnprintf(newValue, sz + 1, "%f", sOptions.adfFixedLevels[i]);
478 : }
479 : else
480 : {
481 18 : CPLsnprintf(newValue, sz + 1, "%f,",
482 18 : sOptions.adfFixedLevels[i]);
483 : }
484 28 : values = values + std::string(newValue);
485 28 : delete[] newValue;
486 : }
487 10 : options = CSLAddString(options, values.c_str());
488 : }
489 :
490 18 : if (sOptions.dfExpBase != 0.0)
491 : {
492 : options =
493 5 : CSLAppendPrintf(options, "LEVEL_EXP_BASE=%f", sOptions.dfExpBase);
494 : }
495 13 : else if (sOptions.dfInterval != 0.0)
496 : {
497 : options =
498 12 : CSLAppendPrintf(options, "LEVEL_INTERVAL=%f", sOptions.dfInterval);
499 : }
500 :
501 18 : if (sOptions.dfOffset != 0.0)
502 : {
503 1 : options = CSLAppendPrintf(options, "LEVEL_BASE=%f", sOptions.dfOffset);
504 : }
505 :
506 18 : if (sOptions.bNoDataSet)
507 : {
508 10 : options = CSLAppendPrintf(options, "NODATA=%.19g", sOptions.dfNoData);
509 : }
510 18 : if (iIDField != -1)
511 : {
512 18 : options = CSLAppendPrintf(options, "ID_FIELD=%d", iIDField);
513 : }
514 18 : if (iElevField != -1)
515 : {
516 9 : options = CSLAppendPrintf(options, "ELEV_FIELD=%d", iElevField);
517 : }
518 18 : if (iElevFieldMin != -1)
519 : {
520 6 : options = CSLAppendPrintf(options, "ELEV_FIELD_MIN=%d", iElevFieldMin);
521 : }
522 18 : if (iElevFieldMax != -1)
523 : {
524 6 : options = CSLAppendPrintf(options, "ELEV_FIELD_MAX=%d", iElevFieldMax);
525 : }
526 18 : if (sOptions.bPolygonize)
527 : {
528 6 : options = CSLAppendPrintf(options, "POLYGONIZE=YES");
529 : }
530 18 : if (sOptions.nGroupTransactions)
531 : {
532 2 : options = CSLAppendPrintf(options, "COMMIT_INTERVAL=" CPL_FRMT_GIB,
533 : sOptions.nGroupTransactions);
534 : }
535 :
536 : CPLErr eErr =
537 18 : GDALContourGenerateEx(hBand, hLayer, options, pfnProgress, nullptr);
538 :
539 18 : CSLDestroy(options);
540 18 : if (GDALClose(hDS) != CE_None)
541 0 : eErr = CE_Failure;
542 18 : GDALClose(hSrcDS);
543 :
544 18 : CSLDestroy(argv);
545 18 : GDALDestroyDriverManager();
546 18 : OGRCleanupAll();
547 :
548 18 : return (eErr == CE_None) ? 0 : 1;
549 : }
550 :
551 0 : MAIN_END
|
