Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Interlis 1 Translator
       4             :  * Purpose:  Implements OGRILI1Layer class.
       5             :  * Author:   Pirmin Kalberer, Sourcepole AG
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2004, Pirmin Kalberer, Sourcepole AG
       9             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_conv.h"
      15             : #include "cpl_string.h"
      16             : #include "ogr_geos.h"
      17             : #include "ogr_ili1.h"
      18             : 
      19             : #include <map>
      20             : #include <vector>
      21             : 
      22             : /************************************************************************/
      23             : /*                           OGRILI1Layer()                              */
      24             : /************************************************************************/
      25             : 
      26         152 : OGRILI1Layer::OGRILI1Layer(OGRFeatureDefn *poFeatureDefnIn,
      27             :                            const GeomFieldInfos &oGeomFieldInfosIn,
      28         152 :                            OGRILI1DataSource *poDSIn)
      29             :     : poFeatureDefn(poFeatureDefnIn), oGeomFieldInfos(oGeomFieldInfosIn),
      30             :       nFeatures(0), papoFeatures(nullptr), nFeatureIdx(0), bGeomsJoined(FALSE),
      31         152 :       poDS(poDSIn)
      32             : {
      33         152 :     SetDescription(poFeatureDefn->GetName());
      34         152 :     poFeatureDefn->Reference();
      35         152 : }
      36             : 
      37             : /************************************************************************/
      38             : /*                           ~OGRILI1Layer()                           */
      39             : /************************************************************************/
      40             : 
      41         304 : OGRILI1Layer::~OGRILI1Layer()
      42             : {
      43         418 :     for (int i = 0; i < nFeatures; i++)
      44             :     {
      45         266 :         delete papoFeatures[i];
      46             :     }
      47         152 :     CPLFree(papoFeatures);
      48             : 
      49         152 :     if (poFeatureDefn)
      50         152 :         poFeatureDefn->Release();
      51         304 : }
      52             : 
      53         266 : OGRErr OGRILI1Layer::AddFeature(OGRFeature *poFeature)
      54             : {
      55         266 :     nFeatures++;
      56             : 
      57         266 :     papoFeatures = static_cast<OGRFeature **>(
      58         266 :         CPLRealloc(papoFeatures, sizeof(void *) * nFeatures));
      59             : 
      60         266 :     papoFeatures[nFeatures - 1] = poFeature;
      61             : 
      62         266 :     return OGRERR_NONE;
      63             : }
      64             : 
      65             : /************************************************************************/
      66             : /*                            ResetReading()                            */
      67             : /************************************************************************/
      68             : 
      69          49 : void OGRILI1Layer::ResetReading()
      70             : {
      71          49 :     nFeatureIdx = 0;
      72          49 : }
      73             : 
      74             : /************************************************************************/
      75             : /*                           GetNextFeature()                           */
      76             : /************************************************************************/
      77             : 
      78          51 : OGRFeature *OGRILI1Layer::GetNextFeature()
      79             : {
      80          51 :     if (!bGeomsJoined)
      81          36 :         JoinGeomLayers();
      82             : 
      83          51 :     while (nFeatureIdx < nFeatures)
      84             :     {
      85          34 :         OGRFeature *poFeature = GetNextFeatureRef();
      86          34 :         if (poFeature)
      87          34 :             return poFeature->Clone();
      88             :     }
      89          17 :     return nullptr;
      90             : }
      91             : 
      92         109 : OGRFeature *OGRILI1Layer::GetNextFeatureRef()
      93             : {
      94         109 :     OGRFeature *poFeature = nullptr;
      95         109 :     if (nFeatureIdx < nFeatures)
      96             :     {
      97         102 :         poFeature = papoFeatures[nFeatureIdx++];
      98             :         // apply filters
      99         204 :         if ((m_poFilterGeom == nullptr ||
     100         204 :              FilterGeometry(poFeature->GetGeometryRef())) &&
     101         102 :             (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
     102         102 :             return poFeature;
     103             :     }
     104           7 :     return nullptr;
     105             : }
     106             : 
     107             : /************************************************************************/
     108             : /*                             GetFeatureRef()                          */
     109             : /************************************************************************/
     110             : 
     111           0 : OGRFeature *OGRILI1Layer::GetFeatureRef(GIntBig nFID)
     112             : 
     113             : {
     114           0 :     ResetReading();
     115             : 
     116           0 :     OGRFeature *poFeature = nullptr;
     117           0 :     while ((poFeature = GetNextFeatureRef()) != nullptr)
     118             :     {
     119           0 :         if (poFeature->GetFID() == nFID)
     120           0 :             return poFeature;
     121             :     }
     122             : 
     123           0 :     return nullptr;
     124             : }
     125             : 
     126          19 : OGRFeature *OGRILI1Layer::GetFeatureRef(const char *fid)
     127             : 
     128             : {
     129          19 :     ResetReading();
     130             : 
     131          19 :     OGRFeature *poFeature = nullptr;
     132          41 :     while ((poFeature = GetNextFeatureRef()) != nullptr)
     133             :     {
     134          41 :         if (!strcmp(poFeature->GetFieldAsString(0), fid))
     135          19 :             return poFeature;
     136             :     }
     137             : 
     138           0 :     return nullptr;
     139             : }
     140             : 
     141             : /************************************************************************/
     142             : /*                          GetFeatureCount()                           */
     143             : /************************************************************************/
     144             : 
     145         189 : GIntBig OGRILI1Layer::GetFeatureCount(int bForce)
     146             : {
     147         189 :     if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr
     148             :         /* && poAreaLineLayer == NULL*/)
     149             :     {
     150         189 :         return nFeatures;
     151             :     }
     152             :     else
     153             :     {
     154           0 :         return OGRLayer::GetFeatureCount(bForce);
     155             :     }
     156             : }
     157             : 
     158             : #ifndef d2str_defined
     159             : #define d2str_defined
     160             : 
     161         346 : static const char *d2str(double val)
     162             : {
     163         346 :     if (val == (int)val)
     164         342 :         return CPLSPrintf("%d", (int)val);
     165           4 :     if (fabs(val) < 370)
     166           4 :         return CPLSPrintf("%.16g", val);
     167           0 :     if (fabs(val) > 100000000.0)
     168           0 :         return CPLSPrintf("%.16g", val);
     169             : 
     170           0 :     return CPLSPrintf("%.3f", val);
     171             : }
     172             : #endif
     173             : 
     174          37 : static void AppendCoordinateList(OGRLineString *poLine, OGRILI1DataSource *poDS)
     175             : {
     176          37 :     const bool b3D = CPL_TO_BOOL(wkbHasZ(poLine->getGeometryType()));
     177             : 
     178         165 :     for (int iPoint = 0; iPoint < poLine->getNumPoints(); iPoint++)
     179             :     {
     180         128 :         if (iPoint == 0)
     181          37 :             VSIFPrintfL(poDS->GetTransferFile(), "STPT");
     182             :         else
     183          91 :             VSIFPrintfL(poDS->GetTransferFile(), "LIPT");
     184         128 :         VSIFPrintfL(poDS->GetTransferFile(), " %s",
     185             :                     d2str(poLine->getX(iPoint)));
     186         128 :         VSIFPrintfL(poDS->GetTransferFile(), " %s",
     187             :                     d2str(poLine->getY(iPoint)));
     188         128 :         if (b3D)
     189          63 :             VSIFPrintfL(poDS->GetTransferFile(), " %s",
     190             :                         d2str(poLine->getZ(iPoint)));
     191         128 :         VSIFPrintfL(poDS->GetTransferFile(), "\n");
     192             :     }
     193          37 :     VSIFPrintfL(poDS->GetTransferFile(), "ELIN\n");
     194          37 : }
     195             : 
     196           1 : static void AppendCompoundCurve(OGRCompoundCurve *poCC, OGRILI1DataSource *poDS)
     197             : {
     198           3 :     for (int iMember = 0; iMember < poCC->getNumCurves(); iMember++)
     199             :     {
     200           2 :         OGRCurve *poGeometry = poCC->getCurve(iMember);
     201           2 :         int b3D = wkbHasZ(poGeometry->getGeometryType());
     202           3 :         int bIsArc = (poGeometry->getGeometryType() == wkbCircularString ||
     203           1 :                       poGeometry->getGeometryType() == wkbCircularStringZ);
     204           2 :         OGRSimpleCurve *poLine = (OGRSimpleCurve *)poGeometry;
     205           7 :         for (int iPoint = 0; iPoint < poLine->getNumPoints(); iPoint++)
     206             :         {
     207             :             // Skip last point in curve member
     208           7 :             if (iPoint == poLine->getNumPoints() - 1 &&
     209           2 :                 iMember < poCC->getNumCurves() - 1)
     210           1 :                 continue;
     211           4 :             if (iMember == 0 && iPoint == 0)
     212           1 :                 VSIFPrintfL(poDS->GetTransferFile(), "STPT");
     213           3 :             else if (bIsArc && iPoint == 1)
     214           1 :                 VSIFPrintfL(poDS->GetTransferFile(), "ARCP");
     215             :             else
     216           2 :                 VSIFPrintfL(poDS->GetTransferFile(), "LIPT");
     217           4 :             VSIFPrintfL(poDS->GetTransferFile(), " %s",
     218             :                         d2str(poLine->getX(iPoint)));
     219           4 :             VSIFPrintfL(poDS->GetTransferFile(), " %s",
     220             :                         d2str(poLine->getY(iPoint)));
     221           4 :             if (b3D)
     222           0 :                 VSIFPrintfL(poDS->GetTransferFile(), " %s",
     223             :                             d2str(poLine->getZ(iPoint)));
     224           4 :             VSIFPrintfL(poDS->GetTransferFile(), "\n");
     225             :         }
     226             :     }
     227           1 :     VSIFPrintfL(poDS->GetTransferFile(), "ELIN\n");
     228           1 : }
     229             : 
     230         101 : int OGRILI1Layer::GeometryAppend(OGRGeometry *poGeometry)
     231             : {
     232             : #ifdef DEBUG_VERBOSE
     233             :     CPLDebug("OGR_ILI", "OGRILI1Layer::GeometryAppend OGRGeometryType: %s",
     234             :              OGRGeometryTypeToName(poGeometry->getGeometryType()));
     235             : #endif
     236             :     /* -------------------------------------------------------------------- */
     237             :     /*      2D Point                                                        */
     238             :     /* -------------------------------------------------------------------- */
     239         101 :     if (poGeometry->getGeometryType() == wkbPoint)
     240             :     {
     241             :         /* embedded in from non-geometry fields */
     242             :     }
     243             :     /* -------------------------------------------------------------------- */
     244             :     /*      3D Point                                                        */
     245             :     /* -------------------------------------------------------------------- */
     246          90 :     else if (poGeometry->getGeometryType() == wkbPoint25D)
     247             :     {
     248             :         /* embedded in from non-geometry fields */
     249             :     }
     250             : 
     251             :     /* -------------------------------------------------------------------- */
     252             :     /*      LineString and LinearRing                                       */
     253             :     /* -------------------------------------------------------------------- */
     254         143 :     else if (poGeometry->getGeometryType() == wkbLineString ||
     255          62 :              poGeometry->getGeometryType() == wkbLineString25D)
     256             :     {
     257          37 :         AppendCoordinateList(poGeometry->toLineString(), poDS);
     258             :     }
     259             : 
     260             :     /* -------------------------------------------------------------------- */
     261             :     /*      Polygon                                                         */
     262             :     /* -------------------------------------------------------------------- */
     263          79 :     else if (poGeometry->getGeometryType() == wkbPolygon ||
     264          35 :              poGeometry->getGeometryType() == wkbPolygon25D)
     265             :     {
     266          18 :         OGRPolygon *poPolygon = poGeometry->toPolygon();
     267          36 :         for (auto &&poRing : poPolygon)
     268             :         {
     269          18 :             if (!GeometryAppend(poRing))
     270           0 :                 return FALSE;
     271             :         }
     272             :     }
     273             : 
     274             :     /* -------------------------------------------------------------------- */
     275             :     /*      MultiPolygon                                                    */
     276             :     /* -------------------------------------------------------------------- */
     277          26 :     else if (wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon ||
     278          20 :              wkbFlatten(poGeometry->getGeometryType()) == wkbMultiLineString ||
     279          14 :              wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPoint ||
     280           8 :              wkbFlatten(poGeometry->getGeometryType()) ==
     281           2 :                  wkbGeometryCollection ||
     282          47 :              wkbFlatten(poGeometry->getGeometryType()) == wkbMultiCurve ||
     283           1 :              wkbFlatten(poGeometry->getGeometryType()) == wkbMultiCurveZ)
     284             :     {
     285          25 :         OGRGeometryCollection *poGC = poGeometry->toGeometryCollection();
     286             : 
     287             : #if 0
     288             :         // TODO: Why this large NOP block?
     289             :         if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon )
     290             :         {
     291             :         }
     292             :         else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiLineString )
     293             :         {
     294             :         }
     295             :         else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPoint )
     296             :         {
     297             :         }
     298             :         else
     299             :         {
     300             :         }
     301             : #endif
     302             : 
     303          62 :         for (auto &&poMember : poGC)
     304             :         {
     305          37 :             if (!GeometryAppend(poMember))
     306           0 :                 return FALSE;
     307             :         }
     308             :     }
     309           1 :     else if (poGeometry->getGeometryType() == wkbCompoundCurve ||
     310           0 :              poGeometry->getGeometryType() == wkbCompoundCurveZ)
     311             :     {
     312           1 :         AppendCompoundCurve(poGeometry->toCompoundCurve(), poDS);
     313             :     }
     314             :     else
     315             :     {
     316           0 :         CPLError(CE_Warning, CPLE_AppDefined,
     317             :                  "Skipping unknown geometry type '%s'",
     318           0 :                  OGRGeometryTypeToName(poGeometry->getGeometryType()));
     319           0 :         return FALSE;
     320             :     }
     321             : 
     322         101 :     return TRUE;
     323             : }
     324             : 
     325             : /************************************************************************/
     326             : /*                           ICreateFeature()                            */
     327             : /************************************************************************/
     328             : 
     329          84 : OGRErr OGRILI1Layer::ICreateFeature(OGRFeature *poFeature)
     330             : {
     331             :     // FIXME: tid not thread safe
     332             :     static long tid = -1;  // system generated TID (must be unique within table)
     333          84 :     VSILFILE *fp = poDS->GetTransferFile();
     334          84 :     VSIFPrintfL(fp, "OBJE");
     335             : 
     336         168 :     if (poFeatureDefn->GetFieldCount() &&
     337          84 :         !EQUAL(poFeatureDefn->GetFieldDefn(0)->GetNameRef(), "TID"))
     338             :     {
     339             :         // Input is not generated from an Interlis 1 source
     340          84 :         if (poFeature->GetFID() != OGRNullFID)
     341           0 :             tid = (int)poFeature->GetFID();
     342             :         else
     343          84 :             ++tid;
     344          84 :         VSIFPrintfL(fp, " %ld", tid);
     345             :         // Embedded geometry
     346          84 :         if (poFeature->GetGeometryRef() != nullptr)
     347             :         {
     348          46 :             OGRGeometry *poGeometry = poFeature->GetGeometryRef();
     349             :             // 2D Point
     350          46 :             if (poGeometry->getGeometryType() == wkbPoint)
     351             :             {
     352           5 :                 OGRPoint *poPoint = poGeometry->toPoint();
     353             : 
     354           5 :                 VSIFPrintfL(fp, " %s", d2str(poPoint->getX()));
     355           5 :                 VSIFPrintfL(fp, " %s", d2str(poPoint->getY()));
     356             :             }
     357             :             // 3D Point
     358          41 :             else if (poGeometry->getGeometryType() == wkbPoint25D)
     359             :             {
     360           3 :                 OGRPoint *poPoint = poGeometry->toPoint();
     361             : 
     362           3 :                 VSIFPrintfL(fp, " %s", d2str(poPoint->getX()));
     363           3 :                 VSIFPrintfL(fp, " %s", d2str(poPoint->getY()));
     364           3 :                 VSIFPrintfL(fp, " %s", d2str(poPoint->getZ()));
     365             :             }
     366             :         }
     367             :     }
     368             : 
     369             :     // Write all fields.
     370         436 :     for (int iField = 0; iField < poFeatureDefn->GetFieldCount(); iField++)
     371             :     {
     372         352 :         if (poFeature->IsFieldSetAndNotNull(iField))
     373             :         {
     374         253 :             const char *pszRaw = poFeature->GetFieldAsString(iField);
     375         253 :             if (poFeatureDefn->GetFieldDefn(iField)->GetType() == OFTString)
     376             :             {
     377             :                 // Interlis 1 encoding is ISO 8859-1 (Latin1) -> Recode from
     378             :                 // UTF-8
     379             :                 char *pszString =
     380          58 :                     CPLRecode(pszRaw, CPL_ENC_UTF8, CPL_ENC_ISO8859_1);
     381             :                 // Replace spaces
     382         234 :                 for (size_t i = 0; pszString[i] != '\0'; i++)
     383             :                 {
     384         176 :                     if (pszString[i] == ' ')
     385           2 :                         pszString[i] = '_';
     386             :                 }
     387          58 :                 VSIFPrintfL(fp, " %s", pszString);
     388          58 :                 CPLFree(pszString);
     389             :             }
     390             :             else
     391             :             {
     392         195 :                 VSIFPrintfL(fp, " %s", pszRaw);
     393             :             }
     394             :         }
     395             :         else
     396             :         {
     397          99 :             VSIFPrintfL(fp, " @");
     398             :         }
     399             :     }
     400          84 :     VSIFPrintfL(fp, "\n");
     401             : 
     402             :     // Write out Geometry
     403          84 :     if (poFeature->GetGeometryRef() != nullptr)
     404             :     {
     405          46 :         GeometryAppend(poFeature->GetGeometryRef());
     406             :     }
     407             : 
     408          84 :     return OGRERR_NONE;
     409             : }
     410             : 
     411             : /************************************************************************/
     412             : /*                           TestCapability()                           */
     413             : /************************************************************************/
     414             : 
     415         107 : int OGRILI1Layer::TestCapability(CPL_UNUSED const char *pszCap)
     416             : {
     417         107 :     if (EQUAL(pszCap, OLCCurveGeometries))
     418          37 :         return TRUE;
     419          70 :     if (EQUAL(pszCap, OLCZGeometries))
     420           0 :         return TRUE;
     421             : 
     422          70 :     if (EQUAL(pszCap, OLCCreateField) || EQUAL(pszCap, OLCSequentialWrite))
     423          32 :         return poDS->GetTransferFile() != nullptr;
     424             : 
     425          38 :     return FALSE;
     426             : }
     427             : 
     428             : /************************************************************************/
     429             : /*                            CreateField()                             */
     430             : /************************************************************************/
     431             : 
     432         112 : OGRErr OGRILI1Layer::CreateField(const OGRFieldDefn *poField,
     433             :                                  int /* bApproxOK */)
     434             : {
     435         112 :     poFeatureDefn->AddFieldDefn(poField);
     436             : 
     437         112 :     return OGRERR_NONE;
     438             : }
     439             : 
     440             : /************************************************************************/
     441             : /*                         Internal routines                            */
     442             : /************************************************************************/
     443             : 
     444          36 : void OGRILI1Layer::JoinGeomLayers()
     445             : {
     446          36 :     bGeomsJoined = true;
     447             : 
     448             :     CPLConfigOptionSetter oSetter("OGR_ARC_STEPSIZE", "0.96",
     449          72 :                                   /* bSetOnlyIfUndefined = */ true);
     450             : 
     451          49 :     for (GeomFieldInfos::const_iterator it = oGeomFieldInfos.begin();
     452          62 :          it != oGeomFieldInfos.end(); ++it)
     453             :     {
     454          13 :         OGRFeatureDefn *geomFeatureDefn = it->second.GetGeomTableDefnRef();
     455          13 :         if (geomFeatureDefn)
     456             :         {
     457          14 :             CPLDebug("OGR_ILI", "Join geometry table %s of field '%s'",
     458          14 :                      geomFeatureDefn->GetName(), it->first.c_str());
     459             :             OGRILI1Layer *poGeomLayer =
     460           7 :                 poDS->GetLayerByName(geomFeatureDefn->GetName());
     461             :             const int nGeomFieldIndex =
     462           7 :                 GetLayerDefn()->GetGeomFieldIndex(it->first.c_str());
     463           7 :             if (it->second.iliGeomType == "Surface")
     464             :             {
     465           5 :                 JoinSurfaceLayer(poGeomLayer, nGeomFieldIndex);
     466             :             }
     467           2 :             else if (it->second.iliGeomType == "Area")
     468             :             {
     469           4 :                 CPLString pointField = it->first + "__Point";
     470             :                 const int nPointFieldIndex =
     471           2 :                     GetLayerDefn()->GetGeomFieldIndex(pointField.c_str());
     472           2 :                 PolygonizeAreaLayer(poGeomLayer, nGeomFieldIndex,
     473             :                                     nPointFieldIndex);
     474             :             }
     475             :         }
     476             :     }
     477          36 : }
     478             : 
     479           5 : void OGRILI1Layer::JoinSurfaceLayer(OGRILI1Layer *poSurfaceLineLayer,
     480             :                                     int nSurfaceFieldIndex)
     481             : {
     482           5 :     CPLDebug("OGR_ILI", "Joining surface layer %s with geometries",
     483           5 :              GetLayerDefn()->GetName());
     484             :     OGRwkbGeometryType geomType =
     485           5 :         GetLayerDefn()->GetGeomFieldDefn(nSurfaceFieldIndex)->GetType();
     486             : 
     487          10 :     std::map<OGRFeature *, std::vector<OGRCurve *>> oMapFeatureToGeomSet;
     488             : 
     489           5 :     poSurfaceLineLayer->ResetReading();
     490             : 
     491             :     // First map: for each target curvepolygon, find all belonging curves
     492          24 :     while (OGRFeature *linefeature = poSurfaceLineLayer->GetNextFeatureRef())
     493             :     {
     494             :         // OBJE entries with same _RefTID are polygon rings of same feature
     495          19 :         OGRFeature *feature = nullptr;
     496          19 :         if (poFeatureDefn->GetFieldDefn(0)->GetType() == OFTString)
     497             :         {
     498          19 :             feature = GetFeatureRef(linefeature->GetFieldAsString(1));
     499             :         }
     500             :         else
     501             :         {
     502           0 :             GIntBig reftid = linefeature->GetFieldAsInteger64(1);
     503           0 :             feature = GetFeatureRef(reftid);
     504             :         }
     505          19 :         if (feature)
     506             :         {
     507          19 :             OGRGeometry *poGeom = linefeature->GetGeomFieldRef(0);
     508          19 :             OGRMultiCurve *curves = dynamic_cast<OGRMultiCurve *>(poGeom);
     509          19 :             if (curves)
     510             :             {
     511          38 :                 for (auto &&curve : curves)
     512             :                 {
     513          19 :                     if (!curve->IsEmpty())
     514          19 :                         oMapFeatureToGeomSet[feature].push_back(curve);
     515             :                 }
     516             :             }
     517             :         }
     518             :         else
     519             :         {
     520           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     521             :                      "Couldn't join feature FID " CPL_FRMT_GIB,
     522             :                      linefeature->GetFieldAsInteger64(1));
     523             :         }
     524          19 :     }
     525             : 
     526             :     // Now for each target polygon, assemble the curves together.
     527             :     std::map<OGRFeature *, std::vector<OGRCurve *>>::const_iterator oIter =
     528           5 :         oMapFeatureToGeomSet.begin();
     529          18 :     for (; oIter != oMapFeatureToGeomSet.end(); ++oIter)
     530             :     {
     531          13 :         OGRFeature *feature = oIter->first;
     532          26 :         std::vector<OGRCurve *> oCurves = oIter->second;
     533             : 
     534          26 :         std::vector<OGRCurve *> oSetDestCurves;
     535          13 :         double dfLargestArea = 0.0;
     536          13 :         OGRCurve *poLargestCurve = nullptr;
     537             :         while (true)
     538             :         {
     539          28 :             std::vector<OGRCurve *>::iterator oIterCurves = oCurves.begin();
     540          28 :             if (oIterCurves == oCurves.end())
     541          13 :                 break;
     542             : 
     543          30 :             OGRPoint endPointCC;
     544          15 :             OGRCompoundCurve *poCC = new OGRCompoundCurve();
     545             : 
     546          15 :             bool bFirst = true;
     547             :             while (true)
     548             :             {
     549          19 :                 bool bNewCurveAdded = false;
     550          19 :                 const double dfEps = 1e-14;
     551          20 :                 for (oIterCurves = oCurves.begin();
     552          21 :                      oIterCurves != oCurves.end(); ++oIterCurves)
     553             :                 {
     554          20 :                     OGRCurve *curve = *oIterCurves;
     555          20 :                     OGRPoint startPoint;
     556          20 :                     OGRPoint endPoint;
     557          20 :                     curve->StartPoint(&startPoint);
     558          20 :                     curve->EndPoint(&endPoint);
     559          25 :                     if (bFirst ||
     560           5 :                         (fabs(startPoint.getX() - endPointCC.getX()) < dfEps &&
     561           4 :                          fabs(startPoint.getY() - endPointCC.getY()) < dfEps))
     562             :                     {
     563          19 :                         bFirst = false;
     564             : 
     565          19 :                         curve->EndPoint(&endPointCC);
     566             : 
     567             :                         const OGRwkbGeometryType eCurveType =
     568          19 :                             wkbFlatten(curve->getGeometryType());
     569          19 :                         if (eCurveType == wkbCompoundCurve)
     570             :                         {
     571             :                             OGRCompoundCurve *poCCSub =
     572          11 :                                 curve->toCompoundCurve();
     573          24 :                             for (auto &&subCurve : poCCSub)
     574             :                             {
     575          13 :                                 poCC->addCurve(subCurve);
     576             :                             }
     577             :                         }
     578             :                         else
     579             :                         {
     580           8 :                             poCC->addCurve(curve);
     581             :                         }
     582          19 :                         oCurves.erase(oIterCurves);
     583          19 :                         bNewCurveAdded = true;
     584          19 :                         break;
     585             :                     }
     586           1 :                     else if (fabs(endPoint.getX() - endPointCC.getX()) <
     587           1 :                                  dfEps &&
     588           0 :                              fabs(endPoint.getY() - endPointCC.getY()) < dfEps)
     589             :                     {
     590           0 :                         curve->StartPoint(&endPointCC);
     591             : 
     592             :                         const OGRwkbGeometryType eCurveType =
     593           0 :                             wkbFlatten(curve->getGeometryType());
     594           0 :                         if (eCurveType == wkbLineString ||
     595             :                             eCurveType == wkbCircularString)
     596             :                         {
     597             :                             OGRSimpleCurve *poSC =
     598           0 :                                 curve->toSimpleCurve()->clone();
     599           0 :                             poSC->reversePoints();
     600           0 :                             poCC->addCurveDirectly(poSC);
     601             :                         }
     602           0 :                         else if (eCurveType == wkbCompoundCurve)
     603             :                         {
     604             :                             // Reverse the order of the elements of the
     605             :                             // compound curve
     606             :                             OGRCompoundCurve *poCCSub =
     607           0 :                                 curve->toCompoundCurve();
     608           0 :                             for (int i = poCCSub->getNumCurves() - 1; i >= 0;
     609             :                                  --i)
     610             :                             {
     611             :                                 OGRSimpleCurve *poSC = poCCSub->getCurve(i)
     612           0 :                                                            ->toSimpleCurve()
     613           0 :                                                            ->clone();
     614           0 :                                 poSC->reversePoints();
     615           0 :                                 poCC->addCurveDirectly(poSC);
     616             :                             }
     617             :                         }
     618             : 
     619           0 :                         oCurves.erase(oIterCurves);
     620           0 :                         bNewCurveAdded = true;
     621           0 :                         break;
     622             :                     }
     623             :                 }
     624          19 :                 if (!bNewCurveAdded || oCurves.empty() || poCC->get_IsClosed())
     625          15 :                     break;
     626           4 :             }
     627             : 
     628          15 :             if (!poCC->get_IsClosed())
     629             :             {
     630           0 :                 char *pszJSon = poCC->exportToJson();
     631           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     632             :                          "A ring %s for feature " CPL_FRMT_GIB " in layer %s "
     633             :                          "was not closed. Dropping it",
     634           0 :                          pszJSon, feature->GetFID(), GetName());
     635           0 :                 delete poCC;
     636           0 :                 CPLFree(pszJSon);
     637             :             }
     638             :             else
     639             :             {
     640          15 :                 double dfArea = poCC->get_Area();
     641          15 :                 if (dfArea >= dfLargestArea)
     642             :                 {
     643          15 :                     dfLargestArea = dfArea;
     644          15 :                     poLargestCurve = poCC;
     645             :                 }
     646          15 :                 oSetDestCurves.push_back(poCC);
     647             :             }
     648          15 :         }
     649             : 
     650             :         // Now build the final polygon by first inserting the largest ring.
     651             :         OGRCurvePolygon *poPoly =
     652          13 :             (geomType == wkbPolygon) ? new OGRPolygon() : new OGRCurvePolygon();
     653          13 :         if (poLargestCurve)
     654             :         {
     655             :             std::vector<OGRCurve *>::iterator oIterCurves =
     656          13 :                 oSetDestCurves.begin();
     657          15 :             for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
     658             :             {
     659          15 :                 OGRCurve *poCurve = *oIterCurves;
     660          15 :                 if (poCurve == poLargestCurve)
     661             :                 {
     662          13 :                     oSetDestCurves.erase(oIterCurves);
     663          13 :                     break;
     664             :                 }
     665             :             }
     666             : 
     667          13 :             if (geomType == wkbPolygon)
     668             :             {
     669           6 :                 poLargestCurve = OGRCurve::CastToLinearRing(poLargestCurve);
     670             :             }
     671          13 :             OGRErr error = poPoly->addRingDirectly(poLargestCurve);
     672          13 :             if (error != OGRERR_NONE)
     673             :             {
     674           0 :                 char *pszJSon = poLargestCurve->exportToJson();
     675           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     676             :                          "Cannot add ring %s to feature " CPL_FRMT_GIB
     677             :                          " in layer %s",
     678           0 :                          pszJSon, feature->GetFID(), GetName());
     679           0 :                 CPLFree(pszJSon);
     680             :             }
     681             : 
     682          13 :             oIterCurves = oSetDestCurves.begin();
     683          15 :             for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
     684             :             {
     685           2 :                 OGRCurve *poCurve = *oIterCurves;
     686           2 :                 if (geomType == wkbPolygon)
     687             :                 {
     688           1 :                     poCurve = OGRCurve::CastToLinearRing(poCurve);
     689             :                 }
     690           2 :                 error = poPoly->addRingDirectly(poCurve);
     691           2 :                 if (error != OGRERR_NONE)
     692             :                 {
     693           0 :                     char *pszJSon = poCurve->exportToJson();
     694           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
     695             :                              "Cannot add ring %s to feature " CPL_FRMT_GIB
     696             :                              " in layer %s",
     697           0 :                              pszJSon, feature->GetFID(), GetName());
     698           0 :                     CPLFree(pszJSon);
     699             :                 }
     700             :             }
     701             :         }
     702             : 
     703          13 :         feature->SetGeomFieldDirectly(nSurfaceFieldIndex, poPoly);
     704             :     }
     705             : 
     706           5 :     ResetReading();
     707           5 : }
     708             : 
     709           2 : OGRMultiPolygon *OGRILI1Layer::Polygonize(OGRGeometryCollection *poLines,
     710             :                                           bool
     711             : #if defined(HAVE_GEOS)
     712             :                                               fix_crossing_lines
     713             : #endif
     714             : )
     715             : {
     716           2 :     if (poLines->getNumGeometries() == 0)
     717             :     {
     718           0 :         return new OGRMultiPolygon();
     719             :     }
     720             : 
     721             : #if defined(HAVE_GEOS)
     722           2 :     GEOSGeom *ahInGeoms = nullptr;
     723           2 :     OGRGeometryCollection *poNoncrossingLines = poLines;
     724           2 :     GEOSGeom hResultGeom = nullptr;
     725           2 :     OGRGeometry *poMP = nullptr;
     726             : 
     727           2 :     if (fix_crossing_lines && poLines->getNumGeometries() > 0)
     728             :     {
     729           0 :         CPLDebug("OGR_ILI", "Fixing crossing lines");
     730             :         // A union of the geometry collection with one line fixes
     731             :         // invalid geometries.
     732           0 :         OGRGeometry *poUnion = poLines->Union(poLines->getGeometryRef(0));
     733           0 :         if (poUnion != nullptr)
     734             :         {
     735           0 :             if (wkbFlatten(poUnion->getGeometryType()) ==
     736           0 :                     wkbGeometryCollection ||
     737           0 :                 wkbFlatten(poUnion->getGeometryType()) == wkbMultiLineString)
     738             :             {
     739           0 :                 poNoncrossingLines =
     740           0 :                     dynamic_cast<OGRGeometryCollection *>(poUnion);
     741           0 :                 CPLDebug("OGR_ILI", "Fixed lines: %d",
     742           0 :                          poNoncrossingLines->getNumGeometries() -
     743           0 :                              poLines->getNumGeometries());
     744             :             }
     745             :             else
     746             :             {
     747           0 :                 delete poUnion;
     748             :             }
     749             :         }
     750             :     }
     751             : 
     752           2 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     753             : 
     754             :     ahInGeoms = static_cast<GEOSGeom *>(
     755           2 :         CPLCalloc(sizeof(void *), poNoncrossingLines->getNumGeometries()));
     756          10 :     for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
     757          16 :         ahInGeoms[i] =
     758           8 :             poNoncrossingLines->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
     759             : 
     760           4 :     hResultGeom = GEOSPolygonize_r(hGEOSCtxt, ahInGeoms,
     761           2 :                                    poNoncrossingLines->getNumGeometries());
     762             : 
     763          10 :     for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
     764           8 :         GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
     765           2 :     CPLFree(ahInGeoms);
     766           2 :     if (poNoncrossingLines != poLines)
     767           0 :         delete poNoncrossingLines;
     768             : 
     769           2 :     if (hResultGeom == nullptr)
     770             :     {
     771           0 :         OGRGeometry::freeGEOSContext(hGEOSCtxt);
     772           0 :         return new OGRMultiPolygon();
     773             :     }
     774             : 
     775           2 :     poMP = OGRGeometryFactory::createFromGEOS(hGEOSCtxt, hResultGeom);
     776             : 
     777           2 :     GEOSGeom_destroy_r(hGEOSCtxt, hResultGeom);
     778           2 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     779             : 
     780           2 :     poMP = OGRGeometryFactory::forceToMultiPolygon(poMP);
     781           2 :     if (poMP && wkbFlatten(poMP->getGeometryType()) == wkbMultiPolygon)
     782           2 :         return dynamic_cast<OGRMultiPolygon *>(poMP);
     783             : 
     784           0 :     delete poMP;
     785           0 :     return new OGRMultiPolygon();
     786             : 
     787             : #else
     788             :     return new OGRMultiPolygon();
     789             : #endif
     790             : }
     791             : 
     792           2 : void OGRILI1Layer::PolygonizeAreaLayer(OGRILI1Layer *poAreaLineLayer,
     793             :                                        int
     794             : #if defined(HAVE_GEOS)
     795             :                                            nAreaFieldIndex
     796             : #endif
     797             :                                        ,
     798             :                                        int
     799             : #if defined(HAVE_GEOS)
     800             :                                            nPointFieldIndex
     801             : #endif
     802             : )
     803             : {
     804             :     // add all lines from poAreaLineLayer to collection
     805           2 :     OGRGeometryCollection *gc = new OGRGeometryCollection();
     806           2 :     poAreaLineLayer->ResetReading();
     807          10 :     while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
     808           8 :         gc->addGeometry(feature->GetGeometryRef());
     809             : 
     810             :     // polygonize lines
     811           4 :     CPLDebug("OGR_ILI", "Polygonizing layer %s with %d multilines",
     812           2 :              poAreaLineLayer->GetLayerDefn()->GetName(),
     813             :              gc->getNumGeometries());
     814           2 :     poAreaLineLayer = nullptr;
     815           2 :     OGRMultiPolygon *polys = Polygonize(gc, false);
     816           2 :     CPLDebug("OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
     817           2 :     if (polys->getNumGeometries() != GetFeatureCount())
     818             :     {
     819           0 :         CPLDebug("OGR_ILI", "Feature count of layer %s: " CPL_FRMT_GIB,
     820           0 :                  GetLayerDefn()->GetName(), GetFeatureCount());
     821           0 :         CPLDebug("OGR_ILI", "Polygonizing again with crossing line fix");
     822           0 :         delete polys;
     823           0 :         polys = Polygonize(gc, true);  // try again with crossing line fix
     824           0 :         CPLDebug("OGR_ILI", "Resulting polygons: %d",
     825             :                  polys->getNumGeometries());
     826             :     }
     827           2 :     delete gc;
     828             : 
     829             :     // associate polygon feature with data row according to centroid
     830             : #if defined(HAVE_GEOS)
     831           4 :     OGRPolygon emptyPoly;
     832             : 
     833           2 :     CPLDebug("OGR_ILI", "Associating layer %s with area polygons",
     834           2 :              GetLayerDefn()->GetName());
     835             :     GEOSGeom *ahInGeoms = static_cast<GEOSGeom *>(
     836           2 :         CPLCalloc(sizeof(void *), polys->getNumGeometries()));
     837           2 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     838           8 :     for (int i = 0; i < polys->getNumGeometries(); i++)
     839             :     {
     840           6 :         ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
     841           6 :         if (!GEOSisValid_r(hGEOSCtxt, ahInGeoms[i]))
     842           0 :             ahInGeoms[i] = nullptr;
     843             :     }
     844           8 :     for (int nFidx = 0; nFidx < nFeatures; nFidx++)
     845             :     {
     846           6 :         OGRFeature *feature = papoFeatures[nFidx];
     847           6 :         OGRGeometry *geomRef = feature->GetGeomFieldRef(nPointFieldIndex);
     848           6 :         if (!geomRef)
     849             :         {
     850           0 :             continue;
     851             :         }
     852             :         GEOSGeom point =
     853           6 :             reinterpret_cast<GEOSGeom>(geomRef->exportToGEOS(hGEOSCtxt));
     854             : 
     855           6 :         int i = 0;
     856          12 :         for (; i < polys->getNumGeometries(); i++)
     857             :         {
     858          12 :             if (ahInGeoms[i] && GEOSWithin_r(hGEOSCtxt, point, ahInGeoms[i]))
     859             :             {
     860           6 :                 feature->SetGeomField(nAreaFieldIndex,
     861           6 :                                       polys->getGeometryRef(i));
     862           6 :                 break;
     863             :             }
     864             :         }
     865           6 :         if (i == polys->getNumGeometries())
     866             :         {
     867           0 :             CPLDebug("OGR_ILI", "Association between area and point failed.");
     868           0 :             feature->SetGeometry(&emptyPoly);
     869             :         }
     870           6 :         GEOSGeom_destroy_r(hGEOSCtxt, point);
     871             :     }
     872           8 :     for (int i = 0; i < polys->getNumGeometries(); i++)
     873           6 :         GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
     874           2 :     CPLFree(ahInGeoms);
     875           2 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     876             : #endif
     877           2 :     delete polys;
     878           2 : }
     879             : 
     880             : /************************************************************************/
     881             : /*                             GetDataset()                             */
     882             : /************************************************************************/
     883             : 
     884          17 : GDALDataset *OGRILI1Layer::GetDataset()
     885             : {
     886          17 :     return poDS;
     887             : }