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          36 :     bool bResetConfigOption = false;
     448          36 :     if (EQUAL(CPLGetConfigOption("OGR_ARC_STEPSIZE", ""), ""))
     449             :     {
     450          36 :         bResetConfigOption = true;
     451          36 :         CPLSetThreadLocalConfigOption("OGR_ARC_STEPSIZE", "0.96");
     452             :     }
     453             : 
     454          49 :     for (GeomFieldInfos::const_iterator it = oGeomFieldInfos.begin();
     455          62 :          it != oGeomFieldInfos.end(); ++it)
     456             :     {
     457          13 :         OGRFeatureDefn *geomFeatureDefn = it->second.GetGeomTableDefnRef();
     458          13 :         if (geomFeatureDefn)
     459             :         {
     460          14 :             CPLDebug("OGR_ILI", "Join geometry table %s of field '%s'",
     461          14 :                      geomFeatureDefn->GetName(), it->first.c_str());
     462             :             OGRILI1Layer *poGeomLayer =
     463           7 :                 poDS->GetLayerByName(geomFeatureDefn->GetName());
     464             :             const int nGeomFieldIndex =
     465           7 :                 GetLayerDefn()->GetGeomFieldIndex(it->first.c_str());
     466           7 :             if (it->second.iliGeomType == "Surface")
     467             :             {
     468           5 :                 JoinSurfaceLayer(poGeomLayer, nGeomFieldIndex);
     469             :             }
     470           2 :             else if (it->second.iliGeomType == "Area")
     471             :             {
     472           4 :                 CPLString pointField = it->first + "__Point";
     473             :                 const int nPointFieldIndex =
     474           2 :                     GetLayerDefn()->GetGeomFieldIndex(pointField.c_str());
     475           2 :                 PolygonizeAreaLayer(poGeomLayer, nGeomFieldIndex,
     476             :                                     nPointFieldIndex);
     477             :             }
     478             :         }
     479             :     }
     480             : 
     481          36 :     if (bResetConfigOption)
     482          36 :         CPLSetThreadLocalConfigOption("OGR_ARC_STEPSIZE", nullptr);
     483          36 : }
     484             : 
     485           5 : void OGRILI1Layer::JoinSurfaceLayer(OGRILI1Layer *poSurfaceLineLayer,
     486             :                                     int nSurfaceFieldIndex)
     487             : {
     488           5 :     CPLDebug("OGR_ILI", "Joining surface layer %s with geometries",
     489           5 :              GetLayerDefn()->GetName());
     490             :     OGRwkbGeometryType geomType =
     491           5 :         GetLayerDefn()->GetGeomFieldDefn(nSurfaceFieldIndex)->GetType();
     492             : 
     493          10 :     std::map<OGRFeature *, std::vector<OGRCurve *>> oMapFeatureToGeomSet;
     494             : 
     495           5 :     poSurfaceLineLayer->ResetReading();
     496             : 
     497             :     // First map: for each target curvepolygon, find all belonging curves
     498          24 :     while (OGRFeature *linefeature = poSurfaceLineLayer->GetNextFeatureRef())
     499             :     {
     500             :         // OBJE entries with same _RefTID are polygon rings of same feature
     501          19 :         OGRFeature *feature = nullptr;
     502          19 :         if (poFeatureDefn->GetFieldDefn(0)->GetType() == OFTString)
     503             :         {
     504          19 :             feature = GetFeatureRef(linefeature->GetFieldAsString(1));
     505             :         }
     506             :         else
     507             :         {
     508           0 :             GIntBig reftid = linefeature->GetFieldAsInteger64(1);
     509           0 :             feature = GetFeatureRef(reftid);
     510             :         }
     511          19 :         if (feature)
     512             :         {
     513          19 :             OGRGeometry *poGeom = linefeature->GetGeomFieldRef(0);
     514          19 :             OGRMultiCurve *curves = dynamic_cast<OGRMultiCurve *>(poGeom);
     515          19 :             if (curves)
     516             :             {
     517          38 :                 for (auto &&curve : curves)
     518             :                 {
     519          19 :                     if (!curve->IsEmpty())
     520          19 :                         oMapFeatureToGeomSet[feature].push_back(curve);
     521             :                 }
     522             :             }
     523             :         }
     524             :         else
     525             :         {
     526           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     527             :                      "Couldn't join feature FID " CPL_FRMT_GIB,
     528             :                      linefeature->GetFieldAsInteger64(1));
     529             :         }
     530          19 :     }
     531             : 
     532             :     // Now for each target polygon, assemble the curves together.
     533             :     std::map<OGRFeature *, std::vector<OGRCurve *>>::const_iterator oIter =
     534           5 :         oMapFeatureToGeomSet.begin();
     535          18 :     for (; oIter != oMapFeatureToGeomSet.end(); ++oIter)
     536             :     {
     537          13 :         OGRFeature *feature = oIter->first;
     538          26 :         std::vector<OGRCurve *> oCurves = oIter->second;
     539             : 
     540          26 :         std::vector<OGRCurve *> oSetDestCurves;
     541          13 :         double dfLargestArea = 0.0;
     542          13 :         OGRCurve *poLargestCurve = nullptr;
     543             :         while (true)
     544             :         {
     545          28 :             std::vector<OGRCurve *>::iterator oIterCurves = oCurves.begin();
     546          28 :             if (oIterCurves == oCurves.end())
     547          13 :                 break;
     548             : 
     549          30 :             OGRPoint endPointCC;
     550          15 :             OGRCompoundCurve *poCC = new OGRCompoundCurve();
     551             : 
     552          15 :             bool bFirst = true;
     553             :             while (true)
     554             :             {
     555          19 :                 bool bNewCurveAdded = false;
     556          19 :                 const double dfEps = 1e-14;
     557          20 :                 for (oIterCurves = oCurves.begin();
     558          21 :                      oIterCurves != oCurves.end(); ++oIterCurves)
     559             :                 {
     560          20 :                     OGRCurve *curve = *oIterCurves;
     561          20 :                     OGRPoint startPoint;
     562          20 :                     OGRPoint endPoint;
     563          20 :                     curve->StartPoint(&startPoint);
     564          20 :                     curve->EndPoint(&endPoint);
     565          25 :                     if (bFirst ||
     566           5 :                         (fabs(startPoint.getX() - endPointCC.getX()) < dfEps &&
     567           4 :                          fabs(startPoint.getY() - endPointCC.getY()) < dfEps))
     568             :                     {
     569          19 :                         bFirst = false;
     570             : 
     571          19 :                         curve->EndPoint(&endPointCC);
     572             : 
     573             :                         const OGRwkbGeometryType eCurveType =
     574          19 :                             wkbFlatten(curve->getGeometryType());
     575          19 :                         if (eCurveType == wkbCompoundCurve)
     576             :                         {
     577             :                             OGRCompoundCurve *poCCSub =
     578          11 :                                 curve->toCompoundCurve();
     579          24 :                             for (auto &&subCurve : poCCSub)
     580             :                             {
     581          13 :                                 poCC->addCurve(subCurve);
     582             :                             }
     583             :                         }
     584             :                         else
     585             :                         {
     586           8 :                             poCC->addCurve(curve);
     587             :                         }
     588          19 :                         oCurves.erase(oIterCurves);
     589          19 :                         bNewCurveAdded = true;
     590          19 :                         break;
     591             :                     }
     592           1 :                     else if (fabs(endPoint.getX() - endPointCC.getX()) <
     593           1 :                                  dfEps &&
     594           0 :                              fabs(endPoint.getY() - endPointCC.getY()) < dfEps)
     595             :                     {
     596           0 :                         curve->StartPoint(&endPointCC);
     597             : 
     598             :                         const OGRwkbGeometryType eCurveType =
     599           0 :                             wkbFlatten(curve->getGeometryType());
     600           0 :                         if (eCurveType == wkbLineString ||
     601             :                             eCurveType == wkbCircularString)
     602             :                         {
     603             :                             OGRSimpleCurve *poSC =
     604           0 :                                 curve->toSimpleCurve()->clone();
     605           0 :                             poSC->reversePoints();
     606           0 :                             poCC->addCurveDirectly(poSC);
     607             :                         }
     608           0 :                         else if (eCurveType == wkbCompoundCurve)
     609             :                         {
     610             :                             // Reverse the order of the elements of the
     611             :                             // compound curve
     612             :                             OGRCompoundCurve *poCCSub =
     613           0 :                                 curve->toCompoundCurve();
     614           0 :                             for (int i = poCCSub->getNumCurves() - 1; i >= 0;
     615             :                                  --i)
     616             :                             {
     617             :                                 OGRSimpleCurve *poSC = poCCSub->getCurve(i)
     618           0 :                                                            ->toSimpleCurve()
     619           0 :                                                            ->clone();
     620           0 :                                 poSC->reversePoints();
     621           0 :                                 poCC->addCurveDirectly(poSC);
     622             :                             }
     623             :                         }
     624             : 
     625           0 :                         oCurves.erase(oIterCurves);
     626           0 :                         bNewCurveAdded = true;
     627           0 :                         break;
     628             :                     }
     629             :                 }
     630          19 :                 if (!bNewCurveAdded || oCurves.empty() || poCC->get_IsClosed())
     631          15 :                     break;
     632           4 :             }
     633             : 
     634          15 :             if (!poCC->get_IsClosed())
     635             :             {
     636           0 :                 char *pszJSon = poCC->exportToJson();
     637           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     638             :                          "A ring %s for feature " CPL_FRMT_GIB " in layer %s "
     639             :                          "was not closed. Dropping it",
     640           0 :                          pszJSon, feature->GetFID(), GetName());
     641           0 :                 delete poCC;
     642           0 :                 CPLFree(pszJSon);
     643             :             }
     644             :             else
     645             :             {
     646          15 :                 double dfArea = poCC->get_Area();
     647          15 :                 if (dfArea >= dfLargestArea)
     648             :                 {
     649          15 :                     dfLargestArea = dfArea;
     650          15 :                     poLargestCurve = poCC;
     651             :                 }
     652          15 :                 oSetDestCurves.push_back(poCC);
     653             :             }
     654          15 :         }
     655             : 
     656             :         // Now build the final polygon by first inserting the largest ring.
     657             :         OGRCurvePolygon *poPoly =
     658          13 :             (geomType == wkbPolygon) ? new OGRPolygon() : new OGRCurvePolygon();
     659          13 :         if (poLargestCurve)
     660             :         {
     661             :             std::vector<OGRCurve *>::iterator oIterCurves =
     662          13 :                 oSetDestCurves.begin();
     663          15 :             for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
     664             :             {
     665          15 :                 OGRCurve *poCurve = *oIterCurves;
     666          15 :                 if (poCurve == poLargestCurve)
     667             :                 {
     668          13 :                     oSetDestCurves.erase(oIterCurves);
     669          13 :                     break;
     670             :                 }
     671             :             }
     672             : 
     673          13 :             if (geomType == wkbPolygon)
     674             :             {
     675           6 :                 poLargestCurve = OGRCurve::CastToLinearRing(poLargestCurve);
     676             :             }
     677          13 :             OGRErr error = poPoly->addRingDirectly(poLargestCurve);
     678          13 :             if (error != OGRERR_NONE)
     679             :             {
     680           0 :                 char *pszJSon = poLargestCurve->exportToJson();
     681           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     682             :                          "Cannot add ring %s to feature " CPL_FRMT_GIB
     683             :                          " in layer %s",
     684           0 :                          pszJSon, feature->GetFID(), GetName());
     685           0 :                 CPLFree(pszJSon);
     686             :             }
     687             : 
     688          13 :             oIterCurves = oSetDestCurves.begin();
     689          15 :             for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
     690             :             {
     691           2 :                 OGRCurve *poCurve = *oIterCurves;
     692           2 :                 if (geomType == wkbPolygon)
     693             :                 {
     694           1 :                     poCurve = OGRCurve::CastToLinearRing(poCurve);
     695             :                 }
     696           2 :                 error = poPoly->addRingDirectly(poCurve);
     697           2 :                 if (error != OGRERR_NONE)
     698             :                 {
     699           0 :                     char *pszJSon = poCurve->exportToJson();
     700           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
     701             :                              "Cannot add ring %s to feature " CPL_FRMT_GIB
     702             :                              " in layer %s",
     703           0 :                              pszJSon, feature->GetFID(), GetName());
     704           0 :                     CPLFree(pszJSon);
     705             :                 }
     706             :             }
     707             :         }
     708             : 
     709          13 :         feature->SetGeomFieldDirectly(nSurfaceFieldIndex, poPoly);
     710             :     }
     711             : 
     712           5 :     ResetReading();
     713           5 : }
     714             : 
     715           2 : OGRMultiPolygon *OGRILI1Layer::Polygonize(OGRGeometryCollection *poLines,
     716             :                                           bool
     717             : #if defined(HAVE_GEOS)
     718             :                                               fix_crossing_lines
     719             : #endif
     720             : )
     721             : {
     722           2 :     if (poLines->getNumGeometries() == 0)
     723             :     {
     724           0 :         return new OGRMultiPolygon();
     725             :     }
     726             : 
     727             : #if defined(HAVE_GEOS)
     728           2 :     GEOSGeom *ahInGeoms = nullptr;
     729           2 :     OGRGeometryCollection *poNoncrossingLines = poLines;
     730           2 :     GEOSGeom hResultGeom = nullptr;
     731           2 :     OGRGeometry *poMP = nullptr;
     732             : 
     733           2 :     if (fix_crossing_lines && poLines->getNumGeometries() > 0)
     734             :     {
     735           0 :         CPLDebug("OGR_ILI", "Fixing crossing lines");
     736             :         // A union of the geometry collection with one line fixes
     737             :         // invalid geometries.
     738           0 :         OGRGeometry *poUnion = poLines->Union(poLines->getGeometryRef(0));
     739           0 :         if (poUnion != nullptr)
     740             :         {
     741           0 :             if (wkbFlatten(poUnion->getGeometryType()) ==
     742           0 :                     wkbGeometryCollection ||
     743           0 :                 wkbFlatten(poUnion->getGeometryType()) == wkbMultiLineString)
     744             :             {
     745           0 :                 poNoncrossingLines =
     746           0 :                     dynamic_cast<OGRGeometryCollection *>(poUnion);
     747           0 :                 CPLDebug("OGR_ILI", "Fixed lines: %d",
     748           0 :                          poNoncrossingLines->getNumGeometries() -
     749           0 :                              poLines->getNumGeometries());
     750             :             }
     751             :             else
     752             :             {
     753           0 :                 delete poUnion;
     754             :             }
     755             :         }
     756             :     }
     757             : 
     758           2 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     759             : 
     760             :     ahInGeoms = static_cast<GEOSGeom *>(
     761           2 :         CPLCalloc(sizeof(void *), poNoncrossingLines->getNumGeometries()));
     762          10 :     for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
     763          16 :         ahInGeoms[i] =
     764           8 :             poNoncrossingLines->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
     765             : 
     766           4 :     hResultGeom = GEOSPolygonize_r(hGEOSCtxt, ahInGeoms,
     767           2 :                                    poNoncrossingLines->getNumGeometries());
     768             : 
     769          10 :     for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
     770           8 :         GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
     771           2 :     CPLFree(ahInGeoms);
     772           2 :     if (poNoncrossingLines != poLines)
     773           0 :         delete poNoncrossingLines;
     774             : 
     775           2 :     if (hResultGeom == nullptr)
     776             :     {
     777           0 :         OGRGeometry::freeGEOSContext(hGEOSCtxt);
     778           0 :         return new OGRMultiPolygon();
     779             :     }
     780             : 
     781           2 :     poMP = OGRGeometryFactory::createFromGEOS(hGEOSCtxt, hResultGeom);
     782             : 
     783           2 :     GEOSGeom_destroy_r(hGEOSCtxt, hResultGeom);
     784           2 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     785             : 
     786           2 :     poMP = OGRGeometryFactory::forceToMultiPolygon(poMP);
     787           2 :     if (poMP && wkbFlatten(poMP->getGeometryType()) == wkbMultiPolygon)
     788           2 :         return dynamic_cast<OGRMultiPolygon *>(poMP);
     789             : 
     790           0 :     delete poMP;
     791           0 :     return new OGRMultiPolygon();
     792             : 
     793             : #else
     794             :     return new OGRMultiPolygon();
     795             : #endif
     796             : }
     797             : 
     798           2 : void OGRILI1Layer::PolygonizeAreaLayer(OGRILI1Layer *poAreaLineLayer,
     799             :                                        int
     800             : #if defined(HAVE_GEOS)
     801             :                                            nAreaFieldIndex
     802             : #endif
     803             :                                        ,
     804             :                                        int
     805             : #if defined(HAVE_GEOS)
     806             :                                            nPointFieldIndex
     807             : #endif
     808             : )
     809             : {
     810             :     // add all lines from poAreaLineLayer to collection
     811           2 :     OGRGeometryCollection *gc = new OGRGeometryCollection();
     812           2 :     poAreaLineLayer->ResetReading();
     813          10 :     while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
     814           8 :         gc->addGeometry(feature->GetGeometryRef());
     815             : 
     816             :     // polygonize lines
     817           4 :     CPLDebug("OGR_ILI", "Polygonizing layer %s with %d multilines",
     818           2 :              poAreaLineLayer->GetLayerDefn()->GetName(),
     819             :              gc->getNumGeometries());
     820           2 :     poAreaLineLayer = nullptr;
     821           2 :     OGRMultiPolygon *polys = Polygonize(gc, false);
     822           2 :     CPLDebug("OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
     823           2 :     if (polys->getNumGeometries() != GetFeatureCount())
     824             :     {
     825           0 :         CPLDebug("OGR_ILI", "Feature count of layer %s: " CPL_FRMT_GIB,
     826           0 :                  GetLayerDefn()->GetName(), GetFeatureCount());
     827           0 :         CPLDebug("OGR_ILI", "Polygonizing again with crossing line fix");
     828           0 :         delete polys;
     829           0 :         polys = Polygonize(gc, true);  // try again with crossing line fix
     830           0 :         CPLDebug("OGR_ILI", "Resulting polygons: %d",
     831             :                  polys->getNumGeometries());
     832             :     }
     833           2 :     delete gc;
     834             : 
     835             :     // associate polygon feature with data row according to centroid
     836             : #if defined(HAVE_GEOS)
     837           4 :     OGRPolygon emptyPoly;
     838             : 
     839           2 :     CPLDebug("OGR_ILI", "Associating layer %s with area polygons",
     840           2 :              GetLayerDefn()->GetName());
     841             :     GEOSGeom *ahInGeoms = static_cast<GEOSGeom *>(
     842           2 :         CPLCalloc(sizeof(void *), polys->getNumGeometries()));
     843           2 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     844           8 :     for (int i = 0; i < polys->getNumGeometries(); i++)
     845             :     {
     846           6 :         ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
     847           6 :         if (!GEOSisValid_r(hGEOSCtxt, ahInGeoms[i]))
     848           0 :             ahInGeoms[i] = nullptr;
     849             :     }
     850           8 :     for (int nFidx = 0; nFidx < nFeatures; nFidx++)
     851             :     {
     852           6 :         OGRFeature *feature = papoFeatures[nFidx];
     853           6 :         OGRGeometry *geomRef = feature->GetGeomFieldRef(nPointFieldIndex);
     854           6 :         if (!geomRef)
     855             :         {
     856           0 :             continue;
     857             :         }
     858             :         GEOSGeom point =
     859           6 :             reinterpret_cast<GEOSGeom>(geomRef->exportToGEOS(hGEOSCtxt));
     860             : 
     861           6 :         int i = 0;
     862          12 :         for (; i < polys->getNumGeometries(); i++)
     863             :         {
     864          12 :             if (ahInGeoms[i] && GEOSWithin_r(hGEOSCtxt, point, ahInGeoms[i]))
     865             :             {
     866           6 :                 feature->SetGeomField(nAreaFieldIndex,
     867           6 :                                       polys->getGeometryRef(i));
     868           6 :                 break;
     869             :             }
     870             :         }
     871           6 :         if (i == polys->getNumGeometries())
     872             :         {
     873           0 :             CPLDebug("OGR_ILI", "Association between area and point failed.");
     874           0 :             feature->SetGeometry(&emptyPoly);
     875             :         }
     876           6 :         GEOSGeom_destroy_r(hGEOSCtxt, point);
     877             :     }
     878           8 :     for (int i = 0; i < polys->getNumGeometries(); i++)
     879           6 :         GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
     880           2 :     CPLFree(ahInGeoms);
     881           2 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     882             : #endif
     883           2 :     delete polys;
     884           2 : }
     885             : 
     886             : /************************************************************************/
     887             : /*                             GetDataset()                             */
     888             : /************************************************************************/
     889             : 
     890          17 : GDALDataset *OGRILI1Layer::GetDataset()
     891             : {
     892          17 :     return poDS;
     893             : }