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