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 "ogr_api.h"
37 #include "ogr_srs_api.h"
38 #include "commonutils.h"
39
40 CPL_CVSID("$Id: gdal_contour.cpp a095bd8fcc94b5d0c06bb08dd638e8da25ea7123 2021-04-29 18:14:50 +0200 Even Rouault $")
41
42 /************************************************************************/
43 /* ArgIsNumeric() */
44 /************************************************************************/
45
ArgIsNumeric(const char * pszArg)46 static bool ArgIsNumeric( const char *pszArg )
47
48 {
49 return CPLGetValueType(pszArg) != CPL_VALUE_STRING;
50 }
51
52 /************************************************************************/
53 /* Usage() */
54 /************************************************************************/
55
Usage(const char * pszErrorMsg=nullptr)56 static void Usage(const char* pszErrorMsg = nullptr)
57
58 {
59 printf(
60 "Usage: gdal_contour [-b <band>] [-a <attribute_name>] [-amin <attribute_name>] [-amax <attribute_name>]\n"
61 " [-3d] [-inodata] [-snodata n] [-f <formatname>] [-i <interval>]\n"
62 " [[-dsco NAME=VALUE] ...] [[-lco NAME=VALUE] ...]\n"
63 " [-off <offset>] [-fl <level> <level>...] [-e <exp_base>]\n"
64 " [-nln <outlayername>] [-q] [-p]\n"
65 " <src_filename> <dst_filename>\n" );
66
67 if( pszErrorMsg != nullptr )
68 fprintf(stderr, "\nFAILURE: %s\n", pszErrorMsg);
69
70 exit( 1 );
71 }
72
CreateElevAttrib(const char * pszElevAttrib,OGRLayerH hLayer)73 static void CreateElevAttrib(const char* pszElevAttrib, OGRLayerH hLayer)
74 {
75 OGRFieldDefnH hFld = OGR_Fld_Create( pszElevAttrib, OFTReal );
76 OGRErr eErr = OGR_L_CreateField( hLayer, hFld, FALSE );
77 OGR_Fld_Destroy( hFld );
78 if( eErr == OGRERR_FAILURE )
79 {
80 exit( 1 );
81 }
82 }
83
84 /************************************************************************/
85 /* main() */
86 /************************************************************************/
87
88 #define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
89 do { if (i + nExtraArg >= argc) \
90 Usage(CPLSPrintf("%s option requires %d argument(s)", \
91 argv[i], nExtraArg)); } while( false )
92
MAIN_START(argc,argv)93 MAIN_START(argc, argv)
94
95 {
96 bool b3D = false;
97 int bNoDataSet = FALSE;
98 bool bIgnoreNoData = false;
99 int nBandIn = 1;
100 double dfInterval = 0.0;
101 double dfNoData = 0.0;
102 double dfOffset = 0.0;
103 double dfExpBase = 0.0;
104 const char *pszSrcFilename = nullptr;
105 const char *pszDstFilename = nullptr;
106 const char *pszElevAttrib = nullptr;
107 const char *pszElevAttribMin = nullptr;
108 const char *pszElevAttribMax = nullptr;
109 const char *pszFormat = nullptr;
110 char **papszDSCO = nullptr;
111 char **papszLCO = nullptr;
112 double adfFixedLevels[1000];
113 int nFixedLevelCount = 0;
114 const char *pszNewLayerName = "contour";
115 bool bQuiet = false;
116 GDALProgressFunc pfnProgress = nullptr;
117 bool bPolygonize = false;
118
119 // Check that we are running against at least GDAL 1.4.
120 // Note to developers: if we use newer API, please change the requirement.
121 if (atoi(GDALVersionInfo("VERSION_NUM")) < 1400)
122 {
123 fprintf(stderr,
124 "At least, GDAL >= 1.4.0 is required for this version of %s, "
125 "which was compiled against GDAL %s\n",
126 argv[0], GDAL_RELEASE_NAME);
127 exit(1);
128 }
129
130 GDALAllRegister();
131 OGRRegisterAll();
132
133 argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
134
135 /* -------------------------------------------------------------------- */
136 /* Parse arguments. */
137 /* -------------------------------------------------------------------- */
138 for( int i = 1; i < argc; i++ )
139 {
140 if( EQUAL(argv[i], "--utility_version") )
141 {
142 printf("%s was compiled against GDAL %s and "
143 "is running against GDAL %s\n",
144 argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
145 CSLDestroy( argv );
146 return 0;
147 }
148 else if( EQUAL(argv[i], "--help") )
149 Usage();
150 else if( EQUAL(argv[i],"-a") )
151 {
152 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
153 // coverity[tainted_data]
154 pszElevAttrib = argv[++i];
155 }
156 else if( EQUAL(argv[i],"-amin") )
157 {
158 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
159 // coverity[tainted_data]
160 pszElevAttribMin = argv[++i];
161 }
162 else if( EQUAL(argv[i],"-amax") )
163 {
164 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
165 // coverity[tainted_data]
166 pszElevAttribMax = argv[++i];
167 }
168 else if( EQUAL(argv[i],"-off") )
169 {
170 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
171 // coverity[tainted_data]
172 dfOffset = CPLAtof(argv[++i]);
173 }
174 else if( EQUAL(argv[i],"-i") )
175 {
176 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
177 // coverity[tainted_data]
178 dfInterval = CPLAtof(argv[++i]);
179 }
180 else if( EQUAL(argv[i],"-e") )
181 {
182 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
183 // coverity[tainted_data]
184 dfExpBase = CPLAtof(argv[++i]);
185 }
186 else if( EQUAL(argv[i],"-p") )
187 {
188 bPolygonize = true;
189 }
190 else if( EQUAL(argv[i],"-fl") )
191 {
192 if( i >= argc-1 )
193 Usage(CPLSPrintf("%s option requires at least 1 argument",
194 argv[i]));
195 while( i < argc-1
196 && nFixedLevelCount
197 < static_cast<int>(sizeof(adfFixedLevels)/sizeof(double))
198 && ArgIsNumeric(argv[i+1]) )
199 // coverity[tainted_data]
200 adfFixedLevels[nFixedLevelCount++] = CPLAtof(argv[++i]);
201 }
202 else if( EQUAL(argv[i],"-b") )
203 {
204 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
205 // coverity[tainted_data]
206 nBandIn = atoi(argv[++i]);
207 }
208 else if( EQUAL(argv[i],"-f") || EQUAL(argv[i],"-of") )
209 {
210 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
211 // coverity[tainted_data]
212 pszFormat = argv[++i];
213 }
214 else if( EQUAL(argv[i],"-dsco") )
215 {
216 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
217 // coverity[tainted_data]
218 papszDSCO = CSLAddString(papszDSCO, argv[++i] );
219 }
220 else if( EQUAL(argv[i],"-lco") )
221 {
222 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
223 // coverity[tainted_data]
224 papszLCO = CSLAddString(papszLCO, argv[++i] );
225 }
226 else if( EQUAL(argv[i],"-3d") )
227 {
228 b3D = true;
229 }
230 else if( EQUAL(argv[i],"-snodata") )
231 {
232 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
233 bNoDataSet = TRUE;
234 // coverity[tainted_data]
235 dfNoData = CPLAtof(argv[++i]);
236 }
237 else if( EQUAL(argv[i],"-nln") )
238 {
239 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
240 // coverity[tainted_data]
241 pszNewLayerName = argv[++i];
242 }
243 else if( EQUAL(argv[i],"-inodata") )
244 {
245 bIgnoreNoData = true;
246 }
247 else if ( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet") )
248 {
249 bQuiet = TRUE;
250 }
251 else if( pszSrcFilename == nullptr )
252 {
253 pszSrcFilename = argv[i];
254 }
255 else if( pszDstFilename == nullptr )
256 {
257 pszDstFilename = argv[i];
258 }
259 else
260 Usage("Too many command options.");
261 }
262
263 if( dfInterval == 0.0 && nFixedLevelCount == 0 && dfExpBase == 0.0 )
264 {
265 Usage("Neither -i nor -fl nor -e are specified.");
266 }
267
268 if (pszSrcFilename == nullptr)
269 {
270 Usage("Missing source filename.");
271 }
272
273 if (pszDstFilename == nullptr)
274 {
275 Usage("Missing destination filename.");
276 }
277
278 if( strcmp(pszDstFilename, "/vsistdout/") == 0 ||
279 strcmp(pszDstFilename, "/dev/stdout") == 0 )
280 {
281 bQuiet = true;
282 }
283
284 if (!bQuiet)
285 pfnProgress = GDALTermProgress;
286
287 /* -------------------------------------------------------------------- */
288 /* Open source raster file. */
289 /* -------------------------------------------------------------------- */
290 GDALDatasetH hSrcDS = GDALOpen(pszSrcFilename, GA_ReadOnly);
291 if( hSrcDS == nullptr )
292 exit( 2 );
293
294 GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, nBandIn );
295 if( hBand == nullptr )
296 {
297 CPLError( CE_Failure, CPLE_AppDefined,
298 "Band %d does not exist on dataset.",
299 nBandIn );
300 exit(2);
301 }
302
303 if( !bNoDataSet && !bIgnoreNoData )
304 dfNoData = GDALGetRasterNoDataValue( hBand, &bNoDataSet );
305
306 /* -------------------------------------------------------------------- */
307 /* Try to get a coordinate system from the raster. */
308 /* -------------------------------------------------------------------- */
309 OGRSpatialReferenceH hSRS = GDALGetSpatialRef( hSrcDS );
310
311 /* -------------------------------------------------------------------- */
312 /* Create the output file. */
313 /* -------------------------------------------------------------------- */
314 CPLString osFormat;
315 if( pszFormat == nullptr )
316 {
317 std::vector<CPLString> aoDrivers =
318 GetOutputDriversFor(pszDstFilename, GDAL_OF_VECTOR);
319 if( aoDrivers.empty() )
320 {
321 CPLError( CE_Failure, CPLE_AppDefined,
322 "Cannot guess driver for %s", pszDstFilename);
323 exit( 10 );
324 }
325 else
326 {
327 if( aoDrivers.size() > 1 )
328 {
329 CPLError( CE_Warning, CPLE_AppDefined,
330 "Several drivers matching %s extension. Using %s",
331 CPLGetExtension(pszDstFilename), aoDrivers[0].c_str() );
332 }
333 osFormat = aoDrivers[0];
334 }
335 }
336 else
337 {
338 osFormat = pszFormat;
339 }
340
341 OGRSFDriverH hDriver = OGRGetDriverByName( osFormat.c_str() );
342
343 if( hDriver == nullptr )
344 {
345 fprintf( stderr, "Unable to find format driver named %s.\n",
346 osFormat.c_str() );
347 exit( 10 );
348 }
349
350 OGRDataSourceH hDS =
351 OGR_Dr_CreateDataSource(hDriver, pszDstFilename, papszDSCO);
352 if( hDS == nullptr )
353 exit(1);
354
355 OGRLayerH hLayer =
356 OGR_DS_CreateLayer(hDS, pszNewLayerName, hSRS,
357 bPolygonize ? (b3D ? wkbMultiPolygon25D : wkbMultiPolygon )
358 : (b3D ? wkbLineString25D : wkbLineString),
359 papszLCO);
360 if( hLayer == nullptr )
361 exit( 1 );
362
363 OGRFieldDefnH hFld = OGR_Fld_Create("ID", OFTInteger);
364 OGR_Fld_SetWidth( hFld, 8 );
365 OGR_L_CreateField( hLayer, hFld, FALSE );
366 OGR_Fld_Destroy( hFld );
367
368 if( bPolygonize )
369 {
370 if( pszElevAttrib )
371 {
372 pszElevAttrib = nullptr;
373 CPLError(CE_Warning, CPLE_NotSupported,
374 "-a is ignored in polygonal contouring mode. "
375 "Use -amin and/or -amax instead");
376 }
377 }
378 else
379 {
380 if( pszElevAttribMin != nullptr || pszElevAttribMax != nullptr )
381 {
382 pszElevAttribMin = nullptr;
383 pszElevAttribMax = nullptr;
384 CPLError(CE_Warning, CPLE_NotSupported,
385 "-amin and/or -amax are ignored in line contouring mode. "
386 "Use -a instead");
387 }
388 }
389
390 if( pszElevAttrib )
391 {
392 CreateElevAttrib( pszElevAttrib, hLayer );
393 }
394
395 if( pszElevAttribMin )
396 {
397 CreateElevAttrib( pszElevAttribMin, hLayer );
398 }
399
400 if( pszElevAttribMax )
401 {
402 CreateElevAttrib( pszElevAttribMax, hLayer );
403 }
404
405 /* -------------------------------------------------------------------- */
406 /* Invoke. */
407 /* -------------------------------------------------------------------- */
408 int iIDField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hLayer ), "ID" );
409 int iElevField = (pszElevAttrib == nullptr) ? -1 :
410 OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hLayer ),
411 pszElevAttrib );
412
413 int iElevFieldMin = (pszElevAttribMin == nullptr) ? -1 :
414 OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hLayer ),
415 pszElevAttribMin );
416
417 int iElevFieldMax = (pszElevAttribMax == nullptr) ? -1 :
418 OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( hLayer ),
419 pszElevAttribMax );
420
421 char** options = nullptr;
422 if ( nFixedLevelCount > 0 ) {
423 std::string values = "FIXED_LEVELS=";
424 for ( int i = 0; i < nFixedLevelCount; i++ ) {
425 const int sz = 32;
426 char* newValue = new char[sz+1];
427 if ( i == nFixedLevelCount - 1 ) {
428 CPLsnprintf( newValue, sz+1, "%f", adfFixedLevels[i] );
429 }
430 else {
431 CPLsnprintf( newValue, sz+1, "%f,", adfFixedLevels[i] );
432 }
433 values = values + std::string( newValue );
434 delete[] newValue;
435 }
436 options = CSLAddString( options, values.c_str() );
437 }
438 else if ( dfExpBase != 0.0 ) {
439 options = CSLAppendPrintf( options, "LEVEL_EXP_BASE=%f", dfExpBase );
440 }
441 else if ( dfInterval != 0.0 ) {
442 options = CSLAppendPrintf( options, "LEVEL_INTERVAL=%f", dfInterval );
443 }
444
445 if ( dfOffset != 0.0 ) {
446 options = CSLAppendPrintf( options, "LEVEL_BASE=%f", dfOffset );
447 }
448
449 if ( bNoDataSet ) {
450 options = CSLAppendPrintf( options, "NODATA=%.19g", dfNoData );
451 }
452 if ( iIDField != -1 ) {
453 options = CSLAppendPrintf( options, "ID_FIELD=%d", iIDField );
454 }
455 if ( iElevField != -1 ) {
456 options = CSLAppendPrintf( options, "ELEV_FIELD=%d", iElevField );
457 }
458 if ( iElevFieldMin != -1 ) {
459 options = CSLAppendPrintf( options, "ELEV_FIELD_MIN=%d", iElevFieldMin );
460 }
461 if ( iElevFieldMax != -1 ) {
462 options = CSLAppendPrintf( options, "ELEV_FIELD_MAX=%d", iElevFieldMax );
463 }
464 if ( bPolygonize ) {
465 options = CSLAppendPrintf( options, "POLYGONIZE=YES" );
466 }
467
468 CPLErr eErr = GDALContourGenerateEx( hBand, hLayer, options, pfnProgress, nullptr );
469
470 CSLDestroy( options );
471 OGR_DS_Destroy( hDS );
472 GDALClose( hSrcDS );
473
474 CSLDestroy( argv );
475 CSLDestroy( papszDSCO );
476 CSLDestroy( papszLCO );
477 GDALDestroyDriverManager();
478 OGRCleanupAll();
479
480 return (eErr == CE_None) ? 0 : 1;
481 }
482 MAIN_END
483