LCOV - code coverage report
Current view: top level - apps - ogrlineref.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 578 861 67.1 %
Date: 2025-02-18 14:19:29 Functions: 12 13 92.3 %

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

Generated by: LCOV version 1.14