1 /******************************************************************************
2 * Project: ogr linear referencing utility
3 * Purpose: main source file
4 * Author: Dmitry Baryshnikov (aka Bishop), polimax@mail.ru
5 *
6 ******************************************************************************
7 * Copyright (C) 2014 NextGIS
8 *
9 * Permission is hereby granted, free of charge, to any person obtaining a
10 * copy of this software and associated documentation files (the "Software"),
11 * to deal in the Software without restriction, including without limitation
12 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13 * and/or sell copies of the Software, and to permit persons to whom the
14 * Software is furnished to do so, subject to the following conditions:
15 *
16 * The above copyright notice and this permission notice shall be included
17 * in all copies or substantial portions of the Software.
18 *
19 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25 * DEALINGS IN THE SOFTWARE.
26 ****************************************************************************/
27
28 #include "cpl_port.h"
29
30 #include "commonutils.h"
31 #include "cpl_conv.h"
32 #include "cpl_error.h"
33 #include "cpl_string.h"
34 #include "gdal_version.h"
35 #include "gdal.h"
36 #include "gdal_alg.h"
37 #include "ogr_api.h"
38 #include "ogr_geos.h"
39 #include "ogr_p.h"
40 #include "ogrsf_frmts.h"
41
42 #include <limits>
43 #include <map>
44 #include <set>
45 #include <vector>
46
47 CPL_CVSID("$Id: ogrlineref.cpp ba6b40f6e6a9d787c27271bc6cacb3fe5733bc27 2019-08-12 22:23:27 +0200 Even Rouault $")
48
49 #if defined(HAVE_GEOS)
50 #if GEOS_VERSION_MAJOR > 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 2)
51 #define HAVE_GEOS_PROJECT
52 #endif
53 #endif
54
55 #define FIELD_START "beg"
56 #define FIELD_FINISH "end"
57 #define FIELD_SCALE_FACTOR "scale"
58 constexpr double DELTA = 0.00000001; // - delta
59 #ifdef HAVE_GEOS_PROJECT
60 constexpr double TOLERANCE_DEGREE = 0.00008983153;
61 constexpr double TOLERANCE_METER = 10.0;
62 #endif
63
64 enum operation
65 {
66 op_unknown = 0,
67 op_create,
68 op_get_pos,
69 op_get_coord,
70 op_get_subline
71 };
72
73 typedef struct _curve_data
74 {
75 OGRLineString* pPart;
76 double dfBeg, dfEnd, dfFactor;
IsInside_curve_data77 bool IsInside( const double& dfDist ) const {
78 return (dfDist + DELTA >= dfBeg) && (dfDist - DELTA <= dfEnd);
79 }
80 } CURVE_DATA;
81
82 /************************************************************************/
83 /* Usage() */
84 /************************************************************************/
85 static void Usage( const char* pszAdditionalMsg,
86 bool bShort = true ) CPL_NO_RETURN;
87
Usage(const char * pszAdditionalMsg,bool bShort)88 static void Usage( const char* pszAdditionalMsg, bool bShort )
89 {
90 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
91
92 printf("Usage: ogrlineref [--help-general] [-progress] [-quiet]\n"
93 " [-f format_name] [[-dsco NAME=VALUE] ...] [[-lco NAME=VALUE]...]\n"
94 " [-create]\n"
95 " [-l src_line_datasource_name] [-ln layer_name] [-lf field_name]\n"
96 " [-p src_repers_datasource_name] [-pn layer_name] [-pm pos_field_name] [-pf field_name]\n"
97 " [-r src_parts_datasource_name] [-rn layer_name]\n"
98 " [-o dst_datasource_name] [-on layer_name] [-of field_name] [-s step]\n"
99 " [-get_pos] [-x long] [-y lat]\n"
100 " [-get_coord] [-m position] \n"
101 " [-get_subline] [-mb position] [-me position]\n");
102
103 if( bShort )
104 {
105 printf("\nNote: ogrlineref --long-usage for full help.\n");
106 if( pszAdditionalMsg )
107 fprintf(stderr, "\nFAILURE: %s\n", pszAdditionalMsg);
108 exit(1);
109 }
110
111 printf(
112 "\n -f format_name: output file format name, possible values are:\n");
113
114 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
115 {
116 GDALDriver *poDriver = poR->GetDriver(iDriver);
117
118 if( CPLTestBool(CSLFetchNameValueDef(poDriver->GetMetadata(),
119 GDAL_DCAP_CREATE, "FALSE")) )
120 printf(" -f \"%s\"\n", poDriver->GetDescription());
121 }
122
123 printf(
124 " -progress: Display progress on terminal. Only works if input layers have the \n"
125 " \"fast feature count\" capability\n"
126 " -dsco NAME=VALUE: Dataset creation option (format specific)\n"
127 " -lco NAME=VALUE: Layer creation option (format specific)\n"
128 " -l src_line_datasource_name: Datasource of line path name\n"
129 " -ln layer_name: Layer name in datasource (optional)\n"
130 " -lf field_name: Field name for unique paths in layer (optional)\n"
131 " -p src_repers_datasource_name: Datasource of repers name\n"
132 " -pn layer_name: Layer name in datasource (optional)\n"
133 " -pm pos_field_name: Line position field name\n"
134 " -pf field_name: Field name for correspondence repers of separate paths in layer (optional)\n"
135 " -r src_parts_datasource_name: Parts datasource name\n"
136 " -rn layer_name: Layer name in datasource (optional)\n"
137 " -o dst_datasource_name: Parts datasource name\n"
138 " -on layer_name: Layer name in datasource (optional)\n"
139 " -of field_name: Field name for correspondence parts of separate paths in layer (optional)\n"
140 " -s step: part size in m\n"
141 );
142
143 if( pszAdditionalMsg )
144 fprintf(stderr, "\nFAILURE: %s\n", pszAdditionalMsg);
145
146 exit(1);
147 }
148
Usage(bool bShort=true)149 static void Usage( bool bShort = true )
150 {
151 Usage(nullptr, bShort);
152 }
153
154 /************************************************************************/
155 /* SetupTargetLayer() */
156 /************************************************************************/
157
SetupTargetLayer(OGRLayer * poSrcLayer,GDALDataset * poDstDS,char ** papszLCO,const char * pszNewLayerName,const char * pszOutputSepFieldName=nullptr)158 static OGRLayer* SetupTargetLayer( OGRLayer * poSrcLayer, GDALDataset *poDstDS,
159 char **papszLCO, const char *pszNewLayerName,
160 const char* pszOutputSepFieldName = nullptr )
161 {
162 const CPLString szLayerName =
163 pszNewLayerName == nullptr
164 ? CPLGetBasename(poDstDS->GetDescription())
165 : pszNewLayerName;
166
167 /* -------------------------------------------------------------------- */
168 /* Get other info. */
169 /* -------------------------------------------------------------------- */
170 OGRFeatureDefn *poSrcFDefn = poSrcLayer->GetLayerDefn();
171
172 /* -------------------------------------------------------------------- */
173 /* Find requested geometry fields. */
174 /* -------------------------------------------------------------------- */
175
176 OGRSpatialReference *poOutputSRS = poSrcLayer->GetSpatialRef();
177
178 /* -------------------------------------------------------------------- */
179 /* Find the layer. */
180 /* -------------------------------------------------------------------- */
181
182 // GetLayerByName() can instantiate layers that would have been
183 // 'hidden' otherwise, for example, non-spatial tables in a
184 // PostGIS-enabled database, so this apparently useless command is
185 // not useless... (#4012)
186 CPLPushErrorHandler(CPLQuietErrorHandler);
187 OGRLayer *poDstLayer = poDstDS->GetLayerByName(szLayerName);
188 CPLPopErrorHandler();
189 CPLErrorReset();
190
191 if( poDstLayer != nullptr )
192 {
193 const int nLayerCount = poDstDS->GetLayerCount();
194 int iLayer = -1; // Used after for.
195 for( iLayer = 0; iLayer < nLayerCount; iLayer++ )
196 {
197 OGRLayer *poLayer = poDstDS->GetLayer(iLayer);
198 if( poLayer == poDstLayer )
199 break;
200 }
201
202 if( iLayer == nLayerCount )
203 // Should not happen with an ideal driver.
204 poDstLayer = nullptr;
205 }
206
207 /* -------------------------------------------------------------------- */
208 /* If the layer does not exist, then create it. */
209 /* -------------------------------------------------------------------- */
210 if( poDstLayer == nullptr )
211 {
212 if( !poDstDS->TestCapability(ODsCCreateLayer) )
213 {
214 fprintf(stderr,
215 "Layer %s not found, and "
216 "CreateLayer not supported by driver.\n",
217 szLayerName.c_str());
218 return nullptr;
219 }
220
221 OGRwkbGeometryType eGType = wkbLineString;
222
223 CPLErrorReset();
224
225 if( poDstDS->TestCapability(ODsCCreateGeomFieldAfterCreateLayer) )
226 {
227 eGType = wkbNone;
228 }
229
230 poDstLayer =
231 poDstDS->CreateLayer(szLayerName, poOutputSRS,
232 static_cast<OGRwkbGeometryType>(eGType),
233 papszLCO);
234
235 if( poDstLayer == nullptr )
236 return nullptr;
237
238 if( poDstDS->TestCapability(ODsCCreateGeomFieldAfterCreateLayer) )
239 {
240 OGRGeomFieldDefn oGFldDefn(poSrcFDefn->GetGeomFieldDefn(0));
241 if( poOutputSRS != nullptr )
242 oGFldDefn.SetSpatialRef(poOutputSRS);
243 oGFldDefn.SetType(wkbLineString);
244 poDstLayer->CreateGeomField(&oGFldDefn);
245 }
246 }
247
248 /* -------------------------------------------------------------------- */
249 /* Otherwise we will append to it, if append was requested. */
250 /* -------------------------------------------------------------------- */
251 else
252 {
253 fprintf(stderr, "FAILED: Layer %s already exists.\n",
254 szLayerName.c_str());
255 return nullptr;
256 }
257
258 // Create beg, end, scale factor fields.
259 OGRFieldDefn oFieldDefn_Beg(FIELD_START, OFTReal);
260 if( poDstLayer->CreateField(&oFieldDefn_Beg) != OGRERR_NONE )
261 {
262 CPLError(CE_Failure, CPLE_AppDefined, "Create %s field failed!",
263 oFieldDefn_Beg.GetNameRef());
264 return nullptr;
265 }
266
267 OGRFieldDefn oFieldDefn_End(FIELD_FINISH, OFTReal);
268 if( poDstLayer->CreateField(&oFieldDefn_End) != OGRERR_NONE )
269 {
270 CPLError(CE_Failure, CPLE_AppDefined, "Create %s field failed!",
271 oFieldDefn_End.GetNameRef());
272 return nullptr;
273 }
274
275 OGRFieldDefn oFieldDefn_SF(FIELD_SCALE_FACTOR, OFTReal);
276 if( poDstLayer->CreateField(&oFieldDefn_SF) != OGRERR_NONE )
277 {
278 CPLError(CE_Failure, CPLE_AppDefined, "Create %s field failed!",
279 oFieldDefn_SF.GetNameRef());
280 return nullptr;
281 }
282
283 if( pszOutputSepFieldName != nullptr )
284 {
285 OGRFieldDefn oSepField(pszOutputSepFieldName, OFTString);
286 oSepField.SetWidth(254);
287 if( poDstLayer->CreateField(&oSepField) != OGRERR_NONE )
288 {
289 CPLError(CE_Failure, CPLE_AppDefined, "Create %s field failed!",
290 oSepField.GetNameRef());
291 return nullptr;
292 }
293 }
294
295 // Now that we've created a field, GetLayerDefn() won't return NULL.
296 OGRFeatureDefn *poDstFDefn = poDstLayer->GetLayerDefn();
297
298 // Sanity check: if it fails, the driver is buggy.
299 if( poDstFDefn != nullptr && poDstFDefn->GetFieldCount() != 3 )
300 {
301 CPLError(CE_Warning, CPLE_AppDefined,
302 "The output driver has claimed to have added the %s field, "
303 "but it did not!",
304 oFieldDefn_Beg.GetNameRef());
305 }
306
307 return poDstLayer;
308 }
309
310 /* -------------------------------------------------------------------- */
311 /* CheckDestDataSourceNameConsistency() */
312 /* -------------------------------------------------------------------- */
313
314 static
CheckDestDataSourceNameConsistency(const char * pszDestFilename,const char * pszDriverName)315 void CheckDestDataSourceNameConsistency(const char* pszDestFilename,
316 const char* pszDriverName)
317 {
318 char* pszDestExtension = CPLStrdup(CPLGetExtension(pszDestFilename));
319
320 // TODO: Would be good to have driver metadata like for GDAL drivers.
321 static const char* apszExtensions[][2] = { { "shp" , "ESRI Shapefile" },
322 { "dbf" , "ESRI Shapefile" },
323 { "sqlite" , "SQLite" },
324 { "db" , "SQLite" },
325 { "mif" , "MapInfo File" },
326 { "tab" , "MapInfo File" },
327 { "s57" , "S57" },
328 { "bna" , "BNA" },
329 { "csv" , "CSV" },
330 { "gml" , "GML" },
331 { "kml" , "KML" },
332 { "kmz" , "LIBKML" },
333 { "json" , "GeoJSON" },
334 { "geojson", "GeoJSON" },
335 { "dxf" , "DXF" },
336 { "gdb" , "FileGDB" },
337 { "pix" , "PCIDSK" },
338 { "sql" , "PGDump" },
339 { "gtm" , "GPSTrackMaker" },
340 { "gmt" , "GMT" },
341 { "pdf" , "PDF" },
342 { nullptr, nullptr }
343 };
344 static const char* apszBeginName[][2] = { { "PG:" , "PG" },
345 { "MySQL:" , "MySQL" },
346 { "CouchDB:" , "CouchDB" },
347 { "GFT:" , "GFT" },
348 { "MSSQL:" , "MSSQLSpatial" },
349 { "ODBC:" , "ODBC" },
350 { "OCI:" , "OCI" },
351 { "SDE:" , "SDE" },
352 { "WFS:" , "WFS" },
353 { nullptr, nullptr }
354 };
355
356 for( int i = 0; apszExtensions[i][0] != nullptr; i++ )
357 {
358 if( EQUAL(pszDestExtension, apszExtensions[i][0]) &&
359 !EQUAL(pszDriverName, apszExtensions[i][1]) )
360 {
361 fprintf(stderr,
362 "Warning: The target file has a '%s' extension, "
363 "which is normally used by the %s driver,\n"
364 "but the requested output driver is %s. "
365 "Is it really what you want ?\n",
366 pszDestExtension,
367 apszExtensions[i][1],
368 pszDriverName);
369 break;
370 }
371 }
372
373 for( int i = 0; apszBeginName[i][0] != nullptr; i++ )
374 {
375 if( EQUALN(pszDestFilename, apszBeginName[i][0],
376 strlen(apszBeginName[i][0])) &&
377 !EQUAL(pszDriverName, apszBeginName[i][1]) )
378 {
379 fprintf(stderr,
380 "Warning: The target file has a name which is normally "
381 "recognized by the %s driver,\n"
382 "but the requested output driver is %s. "
383 "Is it really what you want ?\n",
384 apszBeginName[i][1],
385 pszDriverName);
386 break;
387 }
388 }
389
390 CPLFree(pszDestExtension);
391 }
392
393 //------------------------------------------------------------------------
394 // AddFeature
395 //------------------------------------------------------------------------
396
AddFeature(OGRLayer * const poOutLayer,OGRLineString * pPart,double dfFrom,double dfTo,double dfScaleFactor,bool bQuiet,const char * pszOutputSepFieldName=nullptr,const char * pszOutputSepFieldValue=nullptr)397 static OGRErr AddFeature( OGRLayer* const poOutLayer, OGRLineString* pPart,
398 double dfFrom, double dfTo, double dfScaleFactor,
399 bool bQuiet,
400 const char* pszOutputSepFieldName = nullptr,
401 const char* pszOutputSepFieldValue = nullptr )
402 {
403 OGRFeature *poFeature =
404 OGRFeature::CreateFeature(poOutLayer->GetLayerDefn());
405
406 poFeature->SetField(FIELD_START, dfFrom);
407 poFeature->SetField(FIELD_FINISH, dfTo);
408 poFeature->SetField(FIELD_SCALE_FACTOR, dfScaleFactor);
409
410 if( pszOutputSepFieldName != nullptr )
411 {
412 poFeature->SetField(pszOutputSepFieldName, pszOutputSepFieldValue);
413 }
414
415 poFeature->SetGeometryDirectly(pPart);
416
417 if( poOutLayer->CreateFeature(poFeature) != OGRERR_NONE )
418 {
419 if( !bQuiet )
420 printf("Failed to create feature in shapefile.\n");
421 return OGRERR_FAILURE;
422 }
423
424 OGRFeature::DestroyFeature(poFeature);
425
426 return OGRERR_NONE;
427 }
428
429 //------------------------------------------------------------------------
430 // CreateSubline
431 //------------------------------------------------------------------------
CreateSubline(OGRLayer * const poPkLayer,double dfPosBeg,double dfPosEnd,OGRLayer * const poOutLayer,CPL_UNUSED int bDisplayProgress,bool bQuiet)432 static OGRErr CreateSubline( OGRLayer* const poPkLayer,
433 double dfPosBeg,
434 double dfPosEnd,
435 OGRLayer* const poOutLayer,
436 CPL_UNUSED int bDisplayProgress,
437 bool bQuiet )
438 {
439 // Get step
440 poPkLayer->ResetReading();
441 OGRFeature* pFeature = poPkLayer->GetNextFeature();
442 if( nullptr != pFeature )
443 {
444 // FIXME: Clang Static Analyzer rightly found that the following
445 // code is dead
446 // dfBeg = pFeature->GetFieldAsDouble(FIELD_START);
447 // dfEnd = pFeature->GetFieldAsDouble(FIELD_FINISH);
448 OGRFeature::DestroyFeature(pFeature);
449 }
450 else
451 {
452 fprintf(stderr, "Get step for positions %f - %f failed\n",
453 dfPosBeg, dfPosEnd);
454 return OGRERR_FAILURE;
455 }
456 // Get second part.
457 double dfBeg = 0.0;
458 double dfEnd = 0.0;
459 pFeature = poPkLayer->GetNextFeature();
460 if( nullptr != pFeature )
461 {
462 dfBeg = pFeature->GetFieldAsDouble(FIELD_START);
463 dfEnd = pFeature->GetFieldAsDouble(FIELD_FINISH);
464 OGRFeature::DestroyFeature(pFeature);
465 }
466 else
467 {
468 fprintf(stderr, "Get step for positions %f - %f failed\n",
469 dfPosBeg, dfPosEnd);
470 return OGRERR_FAILURE;
471 }
472 const double dfStep = dfEnd - dfBeg;
473
474 // Round input to step
475 const double dfPosBegLow = floor(dfPosBeg / dfStep) * dfStep;
476 const double dfPosEndHigh = ceil(dfPosEnd / dfStep) * dfStep;
477
478 CPLString szAttributeFilter;
479 szAttributeFilter.Printf("%s >= %f AND %s <= %f",
480 FIELD_START, dfPosBegLow,
481 FIELD_FINISH, dfPosEndHigh);
482 // TODO: ExecuteSQL should be faster.
483 poPkLayer->SetAttributeFilter(szAttributeFilter);
484 poPkLayer->ResetReading();
485
486 std::map<double, OGRFeature *> moParts;
487
488 while( (pFeature = poPkLayer->GetNextFeature()) != nullptr )
489 {
490 double dfStart = pFeature->GetFieldAsDouble(FIELD_START);
491 moParts[dfStart] = pFeature;
492 }
493
494
495 if( moParts.empty() )
496 {
497 fprintf(stderr, "Get parts for positions %f - %f failed\n",
498 dfPosBeg, dfPosEnd);
499 return OGRERR_FAILURE;
500 }
501
502 OGRLineString SubLine;
503 if( moParts.size() == 1 )
504 {
505 std::map<double, OGRFeature *>::iterator IT = moParts.begin();
506 const double dfStart = IT->first;
507 double dfPosBegCorr = dfPosBeg - dfStart;
508 const double dfSF = IT->second->GetFieldAsDouble(FIELD_SCALE_FACTOR);
509 dfPosBegCorr *= dfSF;
510
511 const double dfPosEndCorr = (dfPosEnd - dfStart) * dfSF;
512
513 OGRLineString *pLine = IT->second->GetGeometryRef()->toLineString();
514
515 OGRLineString *pSubLine =
516 pLine->getSubLine(dfPosBegCorr, dfPosEndCorr, FALSE);
517
518 OGRFeature::DestroyFeature(IT->second);
519 // Store.
520 return AddFeature(poOutLayer, pSubLine, dfPosBeg, dfPosEnd,
521 1.0, bQuiet);
522 }
523 else
524 {
525 int nCounter = static_cast<int>(moParts.size());
526 std::map<double, OGRFeature *>::iterator IT = moParts.begin();
527 OGRLineString *pOutLine = new OGRLineString();
528 // Get first part.
529 const double dfStart = IT->first;
530 double dfPosBegCorr = dfPosBeg - dfStart;
531 double dfSF = IT->second->GetFieldAsDouble(FIELD_SCALE_FACTOR);
532 dfPosBegCorr *= dfSF;
533
534 OGRLineString *pLine = IT->second->GetGeometryRef()->toLineString();
535
536 OGRLineString *pSubLine =
537 pLine->getSubLine(dfPosBegCorr, pLine->get_Length(), FALSE);
538
539 pOutLine->addSubLineString(pSubLine);
540 delete pSubLine;
541 OGRFeature::DestroyFeature(IT->second);
542
543 ++IT;
544 nCounter--;
545
546 while( nCounter > 1 )
547 {
548 pLine = IT->second->GetGeometryRef()->toLineString();
549 pOutLine->addSubLineString(pLine);
550 OGRFeature::DestroyFeature(IT->second);
551 ++IT;
552 nCounter--;
553 }
554
555 // Get last part
556 double dfPosEndCorr = dfPosEnd - IT->first;
557 dfSF = IT->second->GetFieldAsDouble(FIELD_SCALE_FACTOR);
558 dfPosEndCorr *= dfSF;
559
560 pLine = IT->second->GetGeometryRef()->toLineString();
561
562 pSubLine = pLine->getSubLine(0, dfPosEndCorr, FALSE);
563
564 pOutLine->addSubLineString(pSubLine);
565 delete pSubLine;
566
567 OGRFeature::DestroyFeature(IT->second);
568 // Store
569 return AddFeature(poOutLayer, pOutLine, dfPosBeg, dfPosEnd,
570 1.0, bQuiet);
571 }
572 }
573
574 //------------------------------------------------------------------------
575 // Project
576 //------------------------------------------------------------------------
577 #ifdef HAVE_GEOS_PROJECT
Project(OGRLineString * pLine,OGRPoint * pPoint)578 static double Project( OGRLineString* pLine, OGRPoint* pPoint )
579 {
580 if( nullptr == pLine || nullptr == pPoint )
581 return -1;
582 OGRPoint TestPoint;
583 pLine->StartPoint(&TestPoint);
584 if( TestPoint.Equals(pPoint) )
585 return 0.0;
586 pLine->EndPoint(&TestPoint);
587 if( TestPoint.Equals(pPoint) )
588 return pLine->get_Length();
589
590 return pLine->Project(pPoint);
591 }
592 #endif
593
594 //------------------------------------------------------------------------
595 // CreatePartsFromLineString
596 //------------------------------------------------------------------------
597 #ifdef HAVE_GEOS_PROJECT
CreatePartsFromLineString(OGRLineString * pPathGeom,OGRLayer * const poPkLayer,int nMValField,double dfStep,OGRLayer * const poOutLayer,int bDisplayProgress,bool bQuiet,const char * pszOutputSepFieldName=nullptr,const char * pszOutputSepFieldValue=nullptr)598 static OGRErr CreatePartsFromLineString(
599 OGRLineString* pPathGeom, OGRLayer* const poPkLayer, int nMValField,
600 double dfStep, OGRLayer* const poOutLayer, int bDisplayProgress,
601 bool bQuiet, const char* pszOutputSepFieldName = nullptr,
602 const char* pszOutputSepFieldValue = nullptr )
603 {
604 // Check repers/milestones/reference points type
605 OGRwkbGeometryType eGeomType = poPkLayer->GetGeomType();
606 if( wkbFlatten(eGeomType) != wkbPoint )
607 {
608 fprintf(stderr, "Unsupported geometry type %s for path\n",
609 OGRGeometryTypeToName(eGeomType));
610 return OGRERR_FAILURE;
611 }
612
613 double dfTolerance = 1.0;
614 OGRSpatialReference* pSpaRef = pPathGeom->getSpatialReference();
615 if( pSpaRef->IsGeographic() )
616 {
617 dfTolerance = TOLERANCE_DEGREE;
618 }
619 else
620 {
621 dfTolerance = TOLERANCE_METER;
622 }
623
624 // Create sorted list of repers.
625 std::map<double, OGRPoint*> moRepers;
626 poPkLayer->ResetReading();
627 OGRFeature* pReperFeature = nullptr;
628 double dfTestDistance = 0.0;
629 while( (pReperFeature = poPkLayer->GetNextFeature()) != nullptr )
630 {
631 const double dfReperPos = pReperFeature->GetFieldAsDouble(nMValField);
632 OGRGeometry* pGeom = pReperFeature->GetGeometryRef();
633 if( nullptr != pGeom )
634 {
635 OGRPoint* pPt = pGeom->clone()->toPoint();
636 if( !bQuiet )
637 {
638 if( moRepers.find(dfReperPos) != moRepers.end() )
639 {
640 CPLError(
641 CE_Warning, CPLE_AppDefined,
642 "The distance %f is already present in repers file!",
643 dfReperPos);
644 }
645 }
646 // Check if reper is inside the path
647 dfTestDistance = Project(pPathGeom, pPt);
648 if( dfTestDistance < 0 )
649 {
650 if( !bQuiet )
651 {
652 CPLError(CE_Warning, CPLE_AppDefined,
653 "The distance %f is out of path!", dfReperPos);
654 }
655 delete pPt;
656 }
657 else
658 {
659 const double dfDist = pPathGeom->Distance(pPt);
660 if( dfDist < dfTolerance )
661 moRepers[dfReperPos] = pPt;
662 else
663 delete pPt;
664 }
665 }
666 OGRFeature::DestroyFeature(pReperFeature);
667 }
668
669 if( moRepers.size() < 2 )
670 {
671 fprintf(stderr, "Not enough repers to proceed\n");
672 return OGRERR_FAILURE;
673 }
674
675 // Check direction.
676 if( !bQuiet )
677 {
678 fprintf(stdout, "Check path direction\n");
679 }
680
681 // Get distance along path from pt1 and pt2.
682 // If pt1 distance > pt2 distance, reverse path
683 std::map<double, OGRPoint*>::const_iterator IT;
684 IT = moRepers.begin();
685 double dfPosition = IT->first;
686 double dfBeginPosition = IT->first;
687 OGRPoint *pt1 = IT->second;
688 ++IT;
689 OGRPoint *pt2 = IT->second;
690
691 double dfDistance1 = Project(pPathGeom, pt1);
692 double dfDistance2 = Project(pPathGeom, pt2);
693
694 if( dfDistance1 > dfDistance2 )
695 {
696 if( !bQuiet )
697 {
698 fprintf(stderr,
699 "Warning: The path is opposite the repers direction. "
700 "Let's reverse path\n");
701 }
702 pPathGeom->reversePoints();
703
704 dfDistance1 = Project(pPathGeom, pt1);
705 dfDistance2 = Project(pPathGeom, pt2);
706 }
707
708 OGRLineString* pPart = nullptr;
709
710 std::vector<CURVE_DATA> astSubLines;
711
712 if( !bQuiet )
713 {
714 fprintf(stdout, "Create parts\n");
715 }
716
717 // Get first part
718 // If first point is not at the beginning of the path
719 // The first part should be from the beginning of the path
720 // to the first point. length == part.getLength
721 OGRPoint *pPtBeg = nullptr;
722 OGRPoint *pPtEnd = nullptr;
723 double dfPtBegPosition = 0.0;
724 double dfPtEndPosition = 0.0;
725
726 if( dfDistance1 > DELTA )
727 {
728 pPart = pPathGeom->getSubLine(0, dfDistance1, FALSE);
729 if( nullptr != pPart )
730 {
731 double dfLen = pPart->get_Length();
732 if( pSpaRef->IsGeographic() )
733 {
734 //convert to UTM/WGS84
735 OGRPoint pt;
736 pPart->Value(dfLen / 2, &pt);
737 const int nZoneEnv =
738 static_cast<int>(30 + (pt.getX() + 3.0) / 6.0 + 0.5);
739 int nEPSG;
740 if( pt.getY() > 0 )
741 {
742 nEPSG = 32600 + nZoneEnv;
743 }
744 else
745 {
746 nEPSG = 32700 + nZoneEnv;
747 }
748 OGRSpatialReference SpatRef;
749 SpatRef.importFromEPSG(nEPSG);
750 OGRGeometry *pTransformPart = pPart->clone();
751 if( pTransformPart->transformTo(&SpatRef) == OGRERR_NONE )
752 {
753 OGRLineString* pTransformPartLS =
754 pTransformPart->toLineString();
755 dfLen = pTransformPartLS->get_Length();
756 }
757
758 CURVE_DATA data = { pPart, dfPosition - dfLen, dfPosition,
759 pPart->get_Length() / dfLen };
760 astSubLines.push_back(data);
761
762 pPtBeg = new OGRPoint();
763 pPart->getPoint(0, pPtBeg);
764 dfPtBegPosition = dfPosition - dfLen;
765
766 // AddFeature(poOutLayer, pPart, dfPosition - dfLen, dfPosition,
767 // pPart->get_Length() / dfLen, bQuiet);
768 delete pTransformPart;
769 }
770 else
771 {
772 CURVE_DATA data = {
773 pPart, dfPosition - dfLen, dfPosition, 1.0
774 };
775 astSubLines.push_back(data);
776 // AddFeature(poOutLayer, pPart, dfPosition - dfLen, dfPosition,
777 // 1.0, bQuiet);
778 pPtBeg = new OGRPoint();
779 pPart->getPoint(0, pPtBeg);
780 dfPtBegPosition = dfPosition - dfLen;
781 }
782 }
783 }
784
785 if( dfDistance2 - dfDistance1 > DELTA )
786 {
787 pPart = pPathGeom->getSubLine(dfDistance1, dfDistance2, FALSE);
788 if( nullptr != pPart )
789 {
790 CURVE_DATA data = {
791 pPart, dfPosition, IT->first,
792 pPart->get_Length() / (IT->first - dfPosition)
793 };
794 astSubLines.push_back(data);
795 // AddFeature(poOutLayer, pPart, dfPosition, IT->first,
796 // pPart->get_Length() / (IT->first - dfPosition),
797 // bQuiet);
798 }
799 }
800
801 GDALProgressFunc pfnProgress = nullptr;
802 void *pProgressArg = nullptr;
803
804 double dfFactor = 1.0 / moRepers.size();
805 if( bDisplayProgress )
806 {
807 pfnProgress = GDALScaledProgress;
808 pProgressArg =
809 GDALCreateScaledProgress(0.0, 1.0, GDALTermProgress, nullptr);
810 }
811
812 int nCount = 2;
813 dfDistance1 = dfDistance2;
814 dfPosition = IT->first;
815 ++IT; // Get third point
816
817 double dfEndPosition = 0.0;
818 while( IT != moRepers.end() )
819 {
820 if( bDisplayProgress )
821 {
822 pfnProgress(nCount * dfFactor, "", pProgressArg);
823 nCount++;
824 }
825
826 dfEndPosition = IT->first;
827
828 dfDistance2 = Project(pPathGeom, IT->second);
829
830 if(dfDistance2 - dfDistance1 > DELTA)
831 {
832 pPart = pPathGeom->getSubLine(dfDistance1, dfDistance2, FALSE);
833 if( nullptr != pPart )
834 {
835 CURVE_DATA data = {
836 pPart, dfPosition, IT->first,
837 pPart->get_Length() / (IT->first - dfPosition)
838 };
839 astSubLines.push_back(data);
840 // AddFeature(poOutLayer, pPart, dfPosition, IT->first,
841 // pPart->get_Length() / (IT->first - dfPosition),
842 // bQuiet);
843 dfDistance1 = dfDistance2;
844 dfPosition = IT->first;
845 }
846 }
847
848 ++IT;
849 }
850
851 // Get last part
852 if( pPathGeom->get_Length() - dfDistance1 > DELTA )
853 {
854 pPart =
855 pPathGeom->getSubLine(dfDistance1, pPathGeom->get_Length(), FALSE);
856 if( nullptr != pPart )
857 {
858 double dfLen = pPart->get_Length();
859 if( pSpaRef->IsGeographic() )
860 {
861 //convert to UTM/WGS84
862 OGRPoint pt;
863 pPart->Value(dfLen / 2, &pt);
864 const int nZoneEnv =
865 static_cast<int>(30 + (pt.getX() + 3.0) / 6.0 + 0.5);
866 int nEPSG;
867 if( pt.getY() > 0 )
868 {
869 nEPSG = 32600 + nZoneEnv;
870 }
871 else
872 {
873 nEPSG = 32700 + nZoneEnv;
874 }
875 OGRSpatialReference SpatRef;
876 SpatRef.importFromEPSG(nEPSG);
877 OGRGeometry *pTransformPart = pPart->clone();
878 if( pTransformPart->transformTo(&SpatRef) == OGRERR_NONE )
879 {
880 OGRLineString* pTransformPartLS =
881 pTransformPart->toLineString();
882 dfLen = pTransformPartLS->get_Length();
883 }
884 CURVE_DATA data = {
885 pPart, dfPosition, dfPosition + dfLen,
886 pPart->get_Length() / dfLen
887 };
888 astSubLines.push_back(data);
889 // AddFeature(poOutLayer, pPart, dfPosition, dfPosition + dfLen,
890 // pPart->get_Length() / dfLen, bQuiet);
891
892 pPtEnd = new OGRPoint();
893 pPart->getPoint(pPart->getNumPoints() - 1, pPtEnd);
894 dfPtEndPosition = dfPosition + dfLen;
895
896 delete pTransformPart;
897 }
898 else
899 {
900 CURVE_DATA data = {
901 pPart, dfPosition, dfPosition + dfLen, 1.0
902 };
903 astSubLines.push_back(data);
904 // AddFeature(poOutLayer, pPart, dfPosition - dfLen, dfPosition,
905 // 1.0, bQuiet);
906 pPtEnd = new OGRPoint();
907 pPart->getPoint(pPart->getNumPoints() - 1, pPtEnd);
908 dfPtEndPosition = dfPosition + dfLen;
909 }
910 }
911 }
912
913 // Create pickets
914 if( !bQuiet )
915 {
916 fprintf(stdout, "\nCreate pickets\n");
917 }
918
919 const double dfRoundBeg =
920 pPtBeg != nullptr
921 ? ceil(dfPtBegPosition / dfStep) * dfStep
922 : ceil(dfBeginPosition / dfStep) * dfStep;
923
924 if( pPtEnd != nullptr )
925 dfEndPosition = dfPtEndPosition;
926
927 dfFactor = dfStep / (dfEndPosition - dfRoundBeg);
928 nCount = 0;
929 for( std::map<double, OGRPoint*>::iterator oIter = moRepers.begin();
930 oIter != moRepers.end();
931 ++oIter )
932 {
933 delete oIter->second;
934 }
935
936 moRepers.clear();
937
938 if( pPtBeg != nullptr )
939 moRepers[dfPtBegPosition] = pPtBeg;
940 if( pPtEnd != nullptr )
941 moRepers[dfPtEndPosition] = pPtEnd;
942
943 for( double dfDist = dfRoundBeg; dfDist <= dfEndPosition; dfDist += dfStep )
944 {
945 if( bDisplayProgress )
946 {
947 pfnProgress(nCount * dfFactor, "", pProgressArg);
948 nCount++;
949 }
950
951 for( int j = 0; j < static_cast<int>(astSubLines.size()); j++ )
952 {
953 if( astSubLines[j].IsInside(dfDist) )
954 {
955 const double dfRealDist =
956 (dfDist - astSubLines[j].dfBeg) * astSubLines[j].dfFactor;
957 OGRPoint *pReperPoint = new OGRPoint();
958 astSubLines[j].pPart->Value(dfRealDist, pReperPoint);
959
960 moRepers[dfDist] = pReperPoint;
961 break;
962 }
963 }
964 }
965
966 for( int i = 0; i < static_cast<int>(astSubLines.size()); i++ )
967 {
968 delete astSubLines[i].pPart;
969 }
970 astSubLines.clear();
971
972 if( !bQuiet )
973 {
974 fprintf(stdout, "\nCreate sublines\n");
975 }
976
977 IT = moRepers.begin();
978 dfFactor = 1.0 / moRepers.size();
979 nCount = 0;
980 dfDistance1 = 0;
981 dfPosition = IT->first;
982
983 while( IT != moRepers.end() )
984 {
985 if( bDisplayProgress )
986 {
987 pfnProgress(nCount * dfFactor, "", pProgressArg);
988 nCount++;
989 }
990
991 dfDistance2 = Project(pPathGeom, IT->second);
992
993 if( dfDistance2 - dfDistance1 > DELTA )
994 {
995 pPart = pPathGeom->getSubLine(dfDistance1, dfDistance2, FALSE);
996 if( nullptr != pPart )
997 {
998 AddFeature(poOutLayer, pPart, dfPosition, IT->first,
999 pPart->get_Length() / (IT->first - dfPosition),
1000 bQuiet, pszOutputSepFieldName,
1001 pszOutputSepFieldValue);
1002 dfDistance1 = dfDistance2;
1003 dfPosition = IT->first;
1004 }
1005 }
1006
1007 ++IT;
1008 }
1009
1010 for( std::map<double, OGRPoint*>::iterator oIter = moRepers.begin();
1011 oIter != moRepers.end();
1012 ++oIter)
1013 {
1014 delete oIter->second;
1015 }
1016
1017 if( !bQuiet )
1018 {
1019 fprintf(stdout, "\nSuccess!\n\n");
1020 }
1021
1022 if( nullptr != pProgressArg )
1023 {
1024 GDALDestroyScaledProgress(pProgressArg);
1025 }
1026
1027 return OGRERR_NONE;
1028 }
1029 #endif
1030
1031 //------------------------------------------------------------------------
1032 // CreateParts
1033 //------------------------------------------------------------------------
1034 #ifdef HAVE_GEOS_PROJECT
CreateParts(OGRLayer * const poLnLayer,OGRLayer * const poPkLayer,int nMValField,double dfStep,OGRLayer * const poOutLayer,int bDisplayProgress,bool bQuiet,const char * pszOutputSepFieldName=nullptr,const char * pszOutputSepFieldValue=nullptr)1035 static OGRErr CreateParts(
1036 OGRLayer* const poLnLayer, OGRLayer* const poPkLayer, int nMValField,
1037 double dfStep, OGRLayer* const poOutLayer, int bDisplayProgress,
1038 bool bQuiet,
1039 const char* pszOutputSepFieldName = nullptr,
1040 const char* pszOutputSepFieldValue = nullptr )
1041 {
1042 OGRErr eRetCode = OGRERR_FAILURE;
1043
1044 //check path and get first line
1045 OGRwkbGeometryType eGeomType = poLnLayer->GetGeomType();
1046 if( wkbFlatten(eGeomType) != wkbLineString &&
1047 wkbFlatten(eGeomType) != wkbMultiLineString )
1048 {
1049 fprintf(stderr, "Unsupported geometry type %s for path\n",
1050 OGRGeometryTypeToName(eGeomType));
1051 return eRetCode;
1052 }
1053
1054 poLnLayer->ResetReading();
1055 // Get first geometry
1056 // TODO: Attribute filter for path geometry.
1057 OGRFeature* pPathFeature = poLnLayer->GetNextFeature();
1058 if( nullptr != pPathFeature )
1059 {
1060 OGRGeometry* pGeom = pPathFeature->GetGeometryRef();
1061
1062 if( pGeom != nullptr &&
1063 wkbFlatten(pGeom->getGeometryType()) == wkbMultiLineString )
1064 {
1065 if( !bQuiet )
1066 {
1067 fprintf(stdout, "\nThe geometry " CPL_FRMT_GIB
1068 " is wkbMultiLineString type\n",
1069 pPathFeature->GetFID());
1070 }
1071
1072 OGRGeometryCollection* pGeomColl = pGeom->toGeometryCollection();
1073 for( int i = 0; i < pGeomColl->getNumGeometries(); ++i )
1074 {
1075 OGRLineString* pPath =
1076 pGeomColl->getGeometryRef(i)->clone()->toLineString();
1077 pPath->assignSpatialReference(pGeomColl->getSpatialReference());
1078 eRetCode = CreatePartsFromLineString(
1079 pPath, poPkLayer, nMValField, dfStep, poOutLayer,
1080 bDisplayProgress, bQuiet, pszOutputSepFieldName,
1081 pszOutputSepFieldValue);
1082
1083 if( eRetCode != OGRERR_NONE )
1084 {
1085 OGRFeature::DestroyFeature(pPathFeature);
1086 return eRetCode;
1087 }
1088 }
1089 }
1090 else if( pGeom != nullptr &&
1091 wkbFlatten(pGeom->getGeometryType()) == wkbLineString )
1092 {
1093 OGRLineString* pGeomClone = pGeom->clone()->toLineString();
1094 eRetCode = CreatePartsFromLineString(
1095 pGeomClone, poPkLayer, nMValField, dfStep, poOutLayer,
1096 bDisplayProgress, bQuiet, pszOutputSepFieldName,
1097 pszOutputSepFieldValue);
1098 delete pGeomClone;
1099 }
1100
1101 OGRFeature::DestroyFeature(pPathFeature);
1102 }
1103
1104 // Should never reach
1105
1106 return eRetCode;
1107 }
1108 #endif
1109
1110 //------------------------------------------------------------------------
1111 // CreatePartsMultiple
1112 //------------------------------------------------------------------------
1113 #ifdef HAVE_GEOS_PROJECT
CreatePartsMultiple(OGRLayer * const poLnLayer,const char * pszLineSepFieldName,OGRLayer * const poPkLayer,const char * pszPicketsSepFieldName,int nMValField,double dfStep,OGRLayer * const poOutLayer,const char * pszOutputSepFieldName,int bDisplayProgress,bool bQuiet)1114 static OGRErr CreatePartsMultiple(
1115 OGRLayer* const poLnLayer, const char* pszLineSepFieldName,
1116 OGRLayer* const poPkLayer, const char* pszPicketsSepFieldName,
1117 int nMValField, double dfStep, OGRLayer* const poOutLayer,
1118 const char* pszOutputSepFieldName, int bDisplayProgress, bool bQuiet )
1119 {
1120 // Read all separate field values into array
1121 OGRFeatureDefn *pDefn = poLnLayer->GetLayerDefn();
1122 const int nLineSepFieldInd = pDefn->GetFieldIndex(pszLineSepFieldName);
1123 if( nLineSepFieldInd == -1 )
1124 {
1125 fprintf(stderr, "The field %s not found\n", pszLineSepFieldName);
1126 return OGRERR_FAILURE;
1127 }
1128
1129 poLnLayer->ResetReading();
1130
1131 std::set<CPLString> asIDs;
1132 OGRFeature* pFeature = nullptr;
1133 while( (pFeature = poLnLayer->GetNextFeature()) != nullptr )
1134 {
1135 CPLString sID = pFeature->GetFieldAsString(nLineSepFieldInd);
1136 asIDs.insert(sID);
1137
1138 OGRFeature::DestroyFeature(pFeature);
1139 }
1140
1141 for( std::set<CPLString>::const_iterator it = asIDs.begin();
1142 it != asIDs.end(); ++it )
1143 {
1144 // Create select clause
1145 CPLString sLineWhere;
1146 sLineWhere.Printf("%s = '%s'", pszLineSepFieldName, it->c_str());
1147 poLnLayer->SetAttributeFilter(sLineWhere);
1148
1149 CPLString sPkWhere;
1150 sPkWhere.Printf("%s = '%s'", pszPicketsSepFieldName, it->c_str());
1151 poPkLayer->SetAttributeFilter(sPkWhere);
1152
1153 if( !bQuiet )
1154 {
1155 fprintf(stdout, "The %s %s\n", pszPicketsSepFieldName, it->c_str());
1156 }
1157
1158 // Don't check success as we want to try all paths
1159 CreateParts(poLnLayer, poPkLayer, nMValField, dfStep, poOutLayer,
1160 bDisplayProgress, bQuiet, pszOutputSepFieldName, *it);
1161 }
1162
1163 return OGRERR_NONE;
1164 }
1165 #endif
1166
1167 //------------------------------------------------------------------------
1168 // GetPosition
1169 //------------------------------------------------------------------------
1170 #ifdef HAVE_GEOS_PROJECT
GetPosition(OGRLayer * const poPkLayer,double dfX,double dfY,int,int bQuiet)1171 static OGRErr GetPosition( OGRLayer* const poPkLayer,
1172 double dfX,
1173 double dfY,
1174 int /* bDisplayProgress */,
1175 int bQuiet)
1176 {
1177 // Create point
1178 OGRPoint pt;
1179 pt.setX(dfX);
1180 pt.setY(dfY);
1181 pt.assignSpatialReference(poPkLayer->GetSpatialRef());
1182
1183 poPkLayer->ResetReading();
1184 OGRLineString *pCloserPart = nullptr;
1185 double dfBeg = 0.0;
1186 double dfScale = 0.0;
1187 double dfMinDistance = std::numeric_limits<double>::max();
1188 OGRFeature* pFeature = nullptr;
1189 while( (pFeature = poPkLayer->GetNextFeature()) != nullptr )
1190 {
1191 OGRGeometry* pCurrentGeom = pFeature->GetGeometryRef();
1192 if( pCurrentGeom != nullptr )
1193 {
1194 double dfCurrentDistance = pCurrentGeom->Distance(&pt);
1195 if( dfCurrentDistance < dfMinDistance )
1196 {
1197 dfMinDistance = dfCurrentDistance;
1198 if( pCloserPart != nullptr )
1199 delete pCloserPart;
1200 pCloserPart = pFeature->StealGeometry()->toLineString();
1201 dfBeg = pFeature->GetFieldAsDouble(FIELD_START);
1202 dfScale = pFeature->GetFieldAsDouble(FIELD_SCALE_FACTOR);
1203 }
1204 }
1205 OGRFeature::DestroyFeature(pFeature);
1206 }
1207
1208 if( nullptr == pCloserPart )
1209 {
1210 fprintf(stderr, "Filed to find closest part\n");
1211 return OGRERR_FAILURE;
1212 }
1213 // Now we have closest part
1214 // Get real distance
1215 const double dfRealDist = Project(pCloserPart, &pt);
1216 delete pCloserPart;
1217 // Compute reference distance
1218 const double dfRefDist = dfBeg + dfRealDist / dfScale;
1219 if( bQuiet )
1220 {
1221 fprintf(stdout, "%s", CPLSPrintf("%f\n", dfRefDist));
1222 }
1223 else
1224 {
1225 fprintf(stdout, "%s", CPLSPrintf(
1226 "The position for coordinates lat:%f, long:%f is %f\n",
1227 dfY, dfX, dfRefDist));
1228 }
1229
1230 return OGRERR_NONE;
1231 }
1232 #endif
1233
1234 //------------------------------------------------------------------------
1235 // GetCoordinates
1236 //------------------------------------------------------------------------
GetCoordinates(OGRLayer * const poPkLayer,double dfPos,int,bool bQuiet)1237 static OGRErr GetCoordinates( OGRLayer* const poPkLayer,
1238 double dfPos,
1239 /* CPL_UNUSED */ int /* bDisplayProgress */,
1240 bool bQuiet )
1241 {
1242 CPLString szAttributeFilter;
1243 szAttributeFilter.Printf(
1244 "%s < %f AND %s > %f", FIELD_START, dfPos, FIELD_FINISH, dfPos);
1245 // TODO: ExecuteSQL should be faster.
1246 poPkLayer->SetAttributeFilter(szAttributeFilter);
1247 poPkLayer->ResetReading();
1248
1249 bool bHaveCoords = false;
1250 OGRFeature* pFeature = nullptr;
1251 while( (pFeature = poPkLayer->GetNextFeature()) != nullptr )
1252 {
1253 bHaveCoords = true;
1254 const double dfStart = pFeature->GetFieldAsDouble(FIELD_START);
1255 double dfPosCorr = dfPos - dfStart;
1256 const double dfSF = pFeature->GetFieldAsDouble(FIELD_SCALE_FACTOR);
1257 dfPosCorr *= dfSF;
1258 OGRLineString *pLine = pFeature->GetGeometryRef()->toLineString();
1259
1260 OGRPoint pt;
1261 pLine->Value(dfPosCorr, &pt);
1262
1263 if( bQuiet )
1264 {
1265 fprintf(stdout, "%s",
1266 CPLSPrintf("%f,%f,%f\n", pt.getX(), pt.getY(), pt.getZ()));
1267 }
1268 else
1269 {
1270 fprintf(stdout, "%s", CPLSPrintf(
1271 "The position for distance %f is lat:%f, long:%f, height:%f\n",
1272 dfPos, pt.getY(), pt.getX(), pt.getZ()));
1273 }
1274 OGRFeature::DestroyFeature(pFeature);
1275 }
1276
1277 if( bHaveCoords )
1278 {
1279 return OGRERR_NONE;
1280 }
1281 else
1282 {
1283 fprintf(stderr, "Get coordinates for position %f failed\n", dfPos);
1284 return OGRERR_FAILURE;
1285 }
1286 }
1287
1288 /************************************************************************/
1289 /* main() */
1290 /************************************************************************/
1291
1292 #define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
1293 do { if( iArg + nExtraArg >= nArgc ) \
1294 Usage(CPLSPrintf("%s option requires %d argument(s)", \
1295 papszArgv[iArg], nExtraArg)); } while( false )
1296
MAIN_START(nArgc,papszArgv)1297 MAIN_START(nArgc, papszArgv)
1298
1299 {
1300 OGRErr eErr = OGRERR_NONE;
1301 bool bQuiet = false;
1302 const char *pszFormat = "ESRI Shapefile";
1303
1304 const char *pszOutputDataSource = nullptr;
1305 const char *pszLineDataSource = nullptr;
1306 const char *pszPicketsDataSource = nullptr;
1307 const char *pszPartsDataSource = nullptr;
1308 char *pszOutputLayerName = nullptr;
1309 const char *pszLineLayerName = nullptr;
1310 #ifdef HAVE_GEOS_PROJECT
1311 const char *pszPicketsLayerName = nullptr;
1312 const char *pszPicketsMField = nullptr;
1313 #endif
1314 const char *pszPartsLayerName = nullptr;
1315
1316 #ifdef HAVE_GEOS_PROJECT
1317 const char *pszLineSepFieldName = nullptr;
1318 const char *pszPicketsSepFieldName = nullptr;
1319 const char *pszOutputSepFieldName = "uniq_uid";
1320 #endif
1321
1322 char **papszDSCO = nullptr;
1323 char **papszLCO = nullptr;
1324
1325 operation stOper = op_unknown;
1326 #ifdef HAVE_GEOS_PROJECT
1327 double dfX = -100000000.0;
1328 double dfY = -100000000.0;
1329 #endif
1330 double dfPos = -100000000.0;
1331
1332 int bDisplayProgress = FALSE;
1333
1334 double dfPosBeg = -100000000.0;
1335 double dfPosEnd = -100000000.0;
1336 #ifdef HAVE_GEOS_PROJECT
1337 double dfStep = -100000000.0;
1338 #endif
1339
1340 // Check strict compilation and runtime library version as we use C++ API.
1341 if( ! GDAL_CHECK_VERSION(papszArgv[0]) )
1342 exit(1);
1343
1344 EarlySetConfigOptions(nArgc, papszArgv);
1345
1346 OGRRegisterAll();
1347
1348 /* -------------------------------------------------------------------- */
1349 /* Processing command line arguments. */
1350 /* -------------------------------------------------------------------- */
1351 nArgc = OGRGeneralCmdLineProcessor( nArgc, &papszArgv, 0 );
1352
1353 if( nArgc < 1 )
1354 exit( -nArgc );
1355
1356 for( int iArg = 1; iArg < nArgc; iArg++ )
1357 {
1358 if( EQUAL(papszArgv[iArg], "--utility_version") )
1359 {
1360 printf("%s was compiled against GDAL %s and "
1361 "is running against GDAL %s\n",
1362 papszArgv[0], GDAL_RELEASE_NAME,
1363 GDALVersionInfo("RELEASE_NAME"));
1364 CSLDestroy( papszArgv );
1365 return 0;
1366 }
1367 else if( EQUAL(papszArgv[iArg],"--help") )
1368 {
1369 Usage();
1370 }
1371 else if( EQUAL(papszArgv[iArg], "--long-usage") )
1372 {
1373 Usage(false);
1374 }
1375
1376 else if( EQUAL(papszArgv[iArg], "-q") ||
1377 EQUAL(papszArgv[iArg], "-quiet") )
1378 {
1379 bQuiet = true;
1380 }
1381 else if( (EQUAL(papszArgv[iArg], "-f") ||
1382 EQUAL(papszArgv[iArg], "-of")) )
1383 {
1384 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1385 // bFormatExplicitlySet = TRUE;
1386 // coverity[tainted_data]
1387 pszFormat = papszArgv[++iArg];
1388 }
1389 else if( EQUAL(papszArgv[iArg], "-dsco") )
1390 {
1391 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1392 // coverity[tainted_data]
1393 papszDSCO = CSLAddString(papszDSCO, papszArgv[++iArg]);
1394 }
1395 else if( EQUAL(papszArgv[iArg],"-lco") )
1396 {
1397 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1398 // coverity[tainted_data]
1399 papszLCO = CSLAddString(papszLCO, papszArgv[++iArg]);
1400 }
1401 else if( EQUAL(papszArgv[iArg], "-create") )
1402 {
1403 stOper = op_create;
1404 }
1405 else if( EQUAL(papszArgv[iArg], "-get_pos") )
1406 {
1407 stOper = op_get_pos;
1408 }
1409 else if( EQUAL(papszArgv[iArg], "-get_coord") )
1410 {
1411 stOper = op_get_coord;
1412 }
1413 else if( EQUAL(papszArgv[iArg], "-get_subline") )
1414 {
1415 stOper = op_get_subline;
1416 }
1417 else if( EQUAL(papszArgv[iArg], "-l") )
1418 {
1419 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1420 // coverity[tainted_data]
1421 pszLineDataSource = papszArgv[++iArg];
1422 }
1423 else if( EQUAL(papszArgv[iArg], "-ln") )
1424 {
1425 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1426 // coverity[tainted_data]
1427 pszLineLayerName = papszArgv[++iArg];
1428 }
1429 else if( EQUAL(papszArgv[iArg], "-lf") )
1430 {
1431 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1432 #ifdef HAVE_GEOS_PROJECT
1433 // coverity[tainted_data]
1434 pszLineSepFieldName = papszArgv[++iArg];
1435 #else
1436 fprintf(stderr,
1437 "GEOS support not enabled or incompatible version.\n");
1438 exit(1);
1439 #endif
1440 }
1441 else if( EQUAL(papszArgv[iArg], "-p") )
1442 {
1443 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1444 // coverity[tainted_data]
1445 pszPicketsDataSource = papszArgv[++iArg];
1446 }
1447 else if( EQUAL(papszArgv[iArg], "-pn") )
1448 {
1449 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1450 #ifdef HAVE_GEOS_PROJECT
1451 // coverity[tainted_data]
1452 pszPicketsLayerName = papszArgv[++iArg];
1453 #else
1454 fprintf(stderr,
1455 "GEOS support not enabled or incompatible version.\n");
1456 exit(1);
1457 #endif
1458 }
1459 else if( EQUAL(papszArgv[iArg], "-pm") )
1460 {
1461 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1462 #ifdef HAVE_GEOS_PROJECT
1463 // coverity[tainted_data]
1464 pszPicketsMField = papszArgv[++iArg];
1465 #else
1466 fprintf(stderr,
1467 "GEOS support not enabled or incompatible version.\n");
1468 exit(1);
1469 #endif
1470 }
1471 else if( EQUAL(papszArgv[iArg], "-pf") )
1472 {
1473 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1474 #ifdef HAVE_GEOS_PROJECT
1475 // coverity[tainted_data]
1476 pszPicketsSepFieldName = papszArgv[++iArg];
1477 #else
1478 fprintf(stderr,
1479 "GEOS support not enabled or incompatible version.\n");
1480 exit(1);
1481 #endif
1482 }
1483 else if( EQUAL(papszArgv[iArg], "-r") )
1484 {
1485 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1486 // coverity[tainted_data]
1487 pszPartsDataSource = papszArgv[++iArg];
1488 }
1489 else if( EQUAL(papszArgv[iArg], "-rn") )
1490 {
1491 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1492 // coverity[tainted_data]
1493 pszPartsLayerName = papszArgv[++iArg];
1494 }
1495 else if( EQUAL(papszArgv[iArg], "-o") )
1496 {
1497 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1498 // coverity[tainted_data]
1499 pszOutputDataSource = papszArgv[++iArg];
1500 }
1501 else if( EQUAL(papszArgv[iArg], "-on") )
1502 {
1503 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1504 // coverity[tainted_data]
1505 pszOutputLayerName = CPLStrdup(papszArgv[++iArg]);
1506 }
1507 else if( EQUAL(papszArgv[iArg], "-of") )
1508 {
1509 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1510 #ifdef HAVE_GEOS_PROJECT
1511 // coverity[tainted_data]
1512 pszOutputSepFieldName = papszArgv[++iArg];
1513 #else
1514 fprintf(stderr,
1515 "GEOS support not enabled or incompatible version.\n");
1516 exit(1);
1517 #endif
1518 }
1519 else if( EQUAL(papszArgv[iArg], "-x") )
1520 {
1521 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1522 #ifdef HAVE_GEOS_PROJECT
1523 // coverity[tainted_data]
1524 dfX = CPLAtofM(papszArgv[++iArg]);
1525 #else
1526 fprintf(stderr,
1527 "GEOS support not enabled or incompatible version.\n");
1528 exit(1);
1529 #endif
1530 }
1531 else if( EQUAL(papszArgv[iArg], "-y") )
1532 {
1533 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1534 #ifdef HAVE_GEOS_PROJECT
1535 // coverity[tainted_data]
1536 dfY = CPLAtofM(papszArgv[++iArg]);
1537 #else
1538 fprintf(stderr,
1539 "GEOS support not enabled or incompatible version.\n");
1540 exit(1);
1541 #endif
1542 }
1543 else if( EQUAL(papszArgv[iArg], "-m") )
1544 {
1545 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1546 // coverity[tainted_data]
1547 dfPos = CPLAtofM(papszArgv[++iArg]);
1548 }
1549 else if( EQUAL(papszArgv[iArg], "-mb") )
1550 {
1551 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1552 // coverity[tainted_data]
1553 dfPosBeg = CPLAtofM(papszArgv[++iArg]);
1554 }
1555 else if( EQUAL(papszArgv[iArg], "-me") )
1556 {
1557 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1558 // coverity[tainted_data]
1559 dfPosEnd = CPLAtofM(papszArgv[++iArg]);
1560 }
1561 else if( EQUAL(papszArgv[iArg], "-s") )
1562 {
1563 CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
1564 #ifdef HAVE_GEOS_PROJECT
1565 // coverity[tainted_data]
1566 dfStep = CPLAtofM(papszArgv[++iArg]);
1567 #else
1568 fprintf(stderr,
1569 "GEOS support not enabled or incompatible version.\n");
1570 exit(1);
1571 #endif
1572 }
1573 else if( EQUAL(papszArgv[iArg], "-progress") )
1574 {
1575 bDisplayProgress = TRUE;
1576 }
1577 else if( papszArgv[iArg][0] == '-' )
1578 {
1579 Usage(CPLSPrintf("Unknown option name '%s'", papszArgv[iArg]));
1580 }
1581 }
1582
1583 if( stOper == op_create )
1584 {
1585 #ifdef HAVE_GEOS_PROJECT
1586 if( pszOutputDataSource == nullptr )
1587 Usage("no output datasource provided");
1588 else if( pszLineDataSource == nullptr )
1589 Usage("no path datasource provided");
1590 else if( pszPicketsDataSource == nullptr )
1591 Usage("no repers datasource provided");
1592 else if( pszPicketsMField == nullptr )
1593 Usage("no position field provided");
1594 else if( dfStep == -100000000.0 )
1595 Usage("no step provided");
1596
1597 /* -------------------------------------------------------------------- */
1598 /* Open data source. */
1599 /* -------------------------------------------------------------------- */
1600
1601 GDALDataset *poLnDS = reinterpret_cast<GDALDataset *>(
1602 OGROpen(pszLineDataSource, FALSE, nullptr));
1603
1604 /* -------------------------------------------------------------------- */
1605 /* Report failure */
1606 /* -------------------------------------------------------------------- */
1607 if( poLnDS == nullptr )
1608 {
1609 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1610
1611 fprintf(stderr, "FAILURE:\n"
1612 "Unable to open path datasource `%s' with "
1613 "the following drivers.\n",
1614 pszLineDataSource);
1615
1616 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1617 {
1618 fprintf(stderr, " -> %s\n",
1619 poR->GetDriver(iDriver)->GetDescription());
1620 }
1621
1622 exit(1);
1623 }
1624
1625 GDALDataset *poPkDS = reinterpret_cast<GDALDataset *>(
1626 OGROpen(pszPicketsDataSource, FALSE, nullptr));
1627 /* -------------------------------------------------------------------- */
1628 /* Report failure */
1629 /* -------------------------------------------------------------------- */
1630 if( poPkDS == nullptr )
1631 {
1632 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1633
1634 fprintf(stderr, "FAILURE:\n"
1635 "Unable to open repers datasource `%s' "
1636 "with the following drivers.\n",
1637 pszPicketsDataSource);
1638
1639 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1640 {
1641 fprintf(stderr, " -> %s\n",
1642 poR->GetDriver(iDriver)->GetDescription());
1643 }
1644
1645 exit(1);
1646 }
1647
1648 /* -------------------------------------------------------------------- */
1649 /* Find the output driver. */
1650 /* -------------------------------------------------------------------- */
1651
1652 if( !bQuiet )
1653 CheckDestDataSourceNameConsistency(pszOutputDataSource, pszFormat);
1654
1655 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1656
1657 GDALDriver *poDriver = poR->GetDriverByName(pszFormat);
1658 if( poDriver == nullptr )
1659 {
1660 fprintf(stderr, "Unable to find driver `%s'.\n", pszFormat);
1661 fprintf(stderr, "The following drivers are available:\n");
1662
1663 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1664 {
1665 fprintf(stderr, " -> `%s'\n",
1666 poR->GetDriver(iDriver)->GetDescription());
1667 }
1668 exit(1);
1669 }
1670
1671 if( !CPLTestBool(CSLFetchNameValueDef(poDriver->GetMetadata(),
1672 GDAL_DCAP_CREATE, "FALSE")) )
1673 {
1674 fprintf(stderr,
1675 "%s driver does not support data source creation.\n",
1676 pszFormat);
1677 exit(1);
1678 }
1679
1680 /* -------------------------------------------------------------------- */
1681 /* Create the output data source. */
1682 /* -------------------------------------------------------------------- */
1683 GDALDataset *poODS = poDriver->Create(
1684 pszOutputDataSource, 0, 0, 0, GDT_Unknown, papszDSCO);
1685 if( poODS == nullptr )
1686 {
1687 fprintf( stderr, "%s driver failed to create %s\n",
1688 pszFormat, pszOutputDataSource );
1689 exit( 1 );
1690 }
1691
1692 OGRLayer *poLnLayer =
1693 pszLineLayerName == nullptr
1694 ? poLnDS->GetLayer(0)
1695 : poLnDS->GetLayerByName(pszLineLayerName);
1696
1697 if( poLnLayer == nullptr )
1698 {
1699 fprintf(stderr, "Get path layer failed.\n");
1700 exit(1);
1701 }
1702
1703 OGRLayer *poPkLayer =
1704 pszPicketsLayerName == nullptr
1705 ? poPkDS->GetLayer(0)
1706 : poPkDS->GetLayerByName(pszPicketsLayerName);
1707
1708 if(poPkLayer == nullptr)
1709 {
1710 fprintf(stderr, "Get repers layer failed.\n");
1711 exit(1);
1712 }
1713
1714 OGRFeatureDefn *poPkFDefn = poPkLayer->GetLayerDefn();
1715 int nMValField = poPkFDefn->GetFieldIndex( pszPicketsMField );
1716
1717 OGRLayer *poOutLayer = nullptr;
1718 if( pszLineSepFieldName != nullptr &&
1719 pszPicketsSepFieldName != nullptr )
1720 {
1721 poOutLayer = SetupTargetLayer(poLnLayer, poODS, papszLCO,
1722 pszOutputLayerName,
1723 pszOutputSepFieldName);
1724 if( poOutLayer == nullptr )
1725 {
1726 fprintf(stderr, "Create output layer failed.\n");
1727 exit(1);
1728 }
1729
1730 // Do the work
1731 eErr = CreatePartsMultiple(
1732 poLnLayer, pszLineSepFieldName, poPkLayer,
1733 pszPicketsSepFieldName, nMValField, dfStep, poOutLayer,
1734 pszOutputSepFieldName, bDisplayProgress, bQuiet);
1735 }
1736 else
1737 {
1738 poOutLayer = SetupTargetLayer(poLnLayer, poODS, papszLCO,
1739 pszOutputLayerName);
1740 if( poOutLayer == nullptr )
1741 {
1742 fprintf(stderr, "Create output layer failed.\n");
1743 exit(1);
1744 }
1745
1746 // Do the work
1747 eErr = CreateParts(poLnLayer, poPkLayer, nMValField, dfStep,
1748 poOutLayer, bDisplayProgress, bQuiet);
1749 }
1750
1751 GDALClose(poLnDS);
1752 GDALClose(poPkDS);
1753 GDALClose(poODS);
1754
1755 if( nullptr != pszOutputLayerName )
1756 CPLFree(pszOutputLayerName);
1757 #else // HAVE_GEOS_PROJECT
1758 fprintf(stderr, "GEOS support not enabled or incompatible version.\n");
1759 exit(1);
1760 #endif // HAVE_GEOS_PROJECT
1761 }
1762 else if(stOper == op_get_pos)
1763 {
1764 #ifdef HAVE_GEOS_PROJECT
1765 if( pszPartsDataSource == nullptr )
1766 Usage("no parts datasource provided");
1767 else if( dfX == -100000000.0 || dfY == -100000000.0 )
1768 Usage("no coordinates provided");
1769
1770 GDALDataset *poPartsDS = reinterpret_cast<GDALDataset *>(
1771 OGROpen(pszPartsDataSource, FALSE, nullptr));
1772 /* -------------------------------------------------------------------- */
1773 /* Report failure */
1774 /* -------------------------------------------------------------------- */
1775 if( poPartsDS == nullptr )
1776 {
1777 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1778
1779 fprintf(stderr, "FAILURE:\n"
1780 "Unable to open parts datasource `%s' with "
1781 "the following drivers.\n",
1782 pszPicketsDataSource);
1783
1784 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1785 {
1786 fprintf(stderr, " -> %s\n",
1787 poR->GetDriver(iDriver)->GetDescription() );
1788 }
1789
1790 exit(1);
1791 }
1792
1793 OGRLayer *poPartsLayer =
1794 pszPartsLayerName == nullptr
1795 ? poPartsDS->GetLayer(0)
1796 : poPartsDS->GetLayerByName(pszPartsLayerName);
1797
1798 if( poPartsLayer == nullptr )
1799 {
1800 fprintf(stderr, "Get parts layer failed.\n");
1801 exit(1);
1802 }
1803
1804 // Do the work
1805 eErr = GetPosition(poPartsLayer, dfX, dfY, bDisplayProgress, bQuiet);
1806
1807 GDALClose(poPartsDS);
1808
1809 #else // HAVE_GEOS_PROJECT
1810 fprintf(stderr, "GEOS support not enabled or incompatible version.\n");
1811 exit(1);
1812 #endif // HAVE_GEOS_PROJECT
1813 }
1814 else if( stOper == op_get_coord )
1815 {
1816 if( pszPartsDataSource == nullptr )
1817 Usage("no parts datasource provided");
1818 else if( dfPos == -100000000.0 )
1819 Usage("no position provided");
1820
1821 GDALDataset *poPartsDS = reinterpret_cast<GDALDataset *>(
1822 OGROpen(pszPartsDataSource, FALSE, nullptr));
1823 /* -------------------------------------------------------------------- */
1824 /* Report failure */
1825 /* -------------------------------------------------------------------- */
1826 if( poPartsDS == nullptr )
1827 {
1828 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1829
1830 fprintf(stderr, "FAILURE:\n"
1831 "Unable to open parts datasource `%s' with "
1832 "the following drivers.\n",
1833 pszPicketsDataSource);
1834
1835 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1836 {
1837 fprintf(stderr, " -> %s\n",
1838 poR->GetDriver(iDriver)->GetDescription());
1839 }
1840
1841 exit(1);
1842 }
1843
1844 OGRLayer *poPartsLayer =
1845 pszPartsLayerName == nullptr
1846 ? poPartsDS->GetLayer(0)
1847 : poPartsDS->GetLayerByName(pszPartsLayerName);
1848
1849 if( poPartsLayer == nullptr )
1850 {
1851 fprintf(stderr, "Get parts layer failed.\n");
1852 exit(1);
1853 }
1854 // Do the work
1855 eErr = GetCoordinates(poPartsLayer, dfPos, bDisplayProgress, bQuiet);
1856
1857 GDALClose(poPartsDS);
1858 }
1859 else if( stOper == op_get_subline )
1860 {
1861 if( pszOutputDataSource == nullptr )
1862 Usage("no output datasource provided");
1863 else if( pszPartsDataSource == nullptr )
1864 Usage("no parts datasource provided");
1865 else if( dfPosBeg == -100000000.0 )
1866 Usage("no begin position provided");
1867 else if( dfPosEnd == -100000000.0 )
1868 Usage("no end position provided");
1869
1870 // Open data source.
1871 GDALDataset *poPartsDS = reinterpret_cast<GDALDataset *>(
1872 OGROpen(pszPartsDataSource, FALSE, nullptr));
1873
1874 // Report failure.
1875 if( poPartsDS == nullptr )
1876 {
1877 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1878
1879 fprintf(stderr, "FAILURE:\n"
1880 "Unable to open parts datasource `%s' with "
1881 "the following drivers.\n",
1882 pszLineDataSource);
1883
1884 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1885 {
1886 fprintf(stderr, " -> %s\n",
1887 poR->GetDriver(iDriver)->GetDescription());
1888 }
1889
1890 exit(1);
1891 }
1892
1893 // Find the output driver.
1894
1895 if( !bQuiet )
1896 CheckDestDataSourceNameConsistency(pszOutputDataSource, pszFormat);
1897
1898 OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
1899
1900 GDALDriver *poDriver = poR->GetDriverByName(pszFormat);
1901 if( poDriver == nullptr )
1902 {
1903 fprintf(stderr, "Unable to find driver `%s'.\n", pszFormat);
1904 fprintf(stderr, "The following drivers are available:\n");
1905
1906 for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
1907 {
1908 fprintf(stderr, " -> `%s'\n",
1909 poR->GetDriver(iDriver)->GetDescription());
1910 }
1911 exit(1);
1912 }
1913
1914 if( !CPLTestBool(CSLFetchNameValueDef(poDriver->GetMetadata(),
1915 GDAL_DCAP_CREATE, "FALSE")) )
1916 {
1917 fprintf(stderr,
1918 "%s driver does not support data source creation.\n",
1919 pszFormat);
1920 exit(1);
1921 }
1922
1923 // Create the output data source.
1924 GDALDataset *poODS =
1925 poDriver->Create(pszOutputDataSource, 0, 0, 0,
1926 GDT_Unknown, papszDSCO);
1927 if( poODS == nullptr )
1928 {
1929 fprintf(stderr, "%s driver failed to create %s\n",
1930 pszFormat, pszOutputDataSource);
1931 exit(1);
1932 }
1933
1934 OGRLayer *poPartsLayer =
1935 pszLineLayerName == nullptr
1936 ? poPartsDS->GetLayer(0)
1937 : poPartsDS->GetLayerByName(pszLineLayerName);
1938
1939 if( poPartsLayer == nullptr )
1940 {
1941 fprintf(stderr, "Get parts layer failed.\n");
1942 exit(1);
1943 }
1944
1945 OGRLayer *poOutLayer =
1946 SetupTargetLayer(poPartsLayer, poODS, papszLCO, pszOutputLayerName);
1947 if( poOutLayer == nullptr )
1948 {
1949 fprintf(stderr, "Create output layer failed.\n");
1950 exit(1);
1951 }
1952
1953 // Do the work.
1954 eErr = CreateSubline(poPartsLayer, dfPosBeg, dfPosEnd, poOutLayer,
1955 bDisplayProgress, bQuiet);
1956
1957 GDALClose(poPartsDS);
1958 GDALClose(poODS);
1959
1960 if( nullptr != pszOutputLayerName )
1961 CPLFree(pszOutputLayerName);
1962 }
1963 else
1964 {
1965 Usage("no operation provided");
1966 }
1967
1968 CSLDestroy( papszArgv );
1969 CSLDestroy( papszDSCO );
1970 CSLDestroy( papszLCO );
1971
1972 OGRCleanupAll();
1973
1974 #ifdef DBMALLOC
1975 malloc_dump(1);
1976 #endif
1977
1978 return eErr == OGRERR_NONE ? 0 : 1;
1979 }
1980 MAIN_END
1981