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          46 : OGRILI1Layer::OGRILI1Layer(OGRFeatureDefn *poFeatureDefnIn,
      27             :                            const GeomFieldInfos &oGeomFieldInfosIn,
      28          46 :                            OGRILI1DataSource *poDSIn)
      29             :     : poFeatureDefn(poFeatureDefnIn), oGeomFieldInfos(oGeomFieldInfosIn),
      30             :       nFeatures(0), papoFeatures(nullptr), nFeatureIdx(0), bGeomsJoined(FALSE),
      31          46 :       poDS(poDSIn)
      32             : {
      33          46 :     SetDescription(poFeatureDefn->GetName());
      34          46 :     poFeatureDefn->Reference();
      35          46 : }
      36             : 
      37             : /************************************************************************/
      38             : /*                           ~OGRILI1Layer()                           */
      39             : /************************************************************************/
      40             : 
      41          92 : OGRILI1Layer::~OGRILI1Layer()
      42             : {
      43         146 :     for (int i = 0; i < nFeatures; i++)
      44             :     {
      45         100 :         delete papoFeatures[i];
      46             :     }
      47          46 :     CPLFree(papoFeatures);
      48             : 
      49          46 :     if (poFeatureDefn)
      50          46 :         poFeatureDefn->Release();
      51          92 : }
      52             : 
      53         100 : OGRErr OGRILI1Layer::AddFeature(OGRFeature *poFeature)
      54             : {
      55         100 :     nFeatures++;
      56             : 
      57         100 :     papoFeatures = static_cast<OGRFeature **>(
      58         100 :         CPLRealloc(papoFeatures, sizeof(void *) * nFeatures));
      59             : 
      60         100 :     papoFeatures[nFeatures - 1] = poFeature;
      61             : 
      62         100 :     return OGRERR_NONE;
      63             : }
      64             : 
      65             : /************************************************************************/
      66             : /*                            ResetReading()                            */
      67             : /************************************************************************/
      68             : 
      69          33 : void OGRILI1Layer::ResetReading()
      70             : {
      71          33 :     nFeatureIdx = 0;
      72          33 : }
      73             : 
      74             : /************************************************************************/
      75             : /*                           GetNextFeature()                           */
      76             : /************************************************************************/
      77             : 
      78          32 : OGRFeature *OGRILI1Layer::GetNextFeature()
      79             : {
      80          32 :     if (!bGeomsJoined)
      81          17 :         JoinGeomLayers();
      82             : 
      83          32 :     while (nFeatureIdx < nFeatures)
      84             :     {
      85          31 :         OGRFeature *poFeature = GetNextFeatureRef();
      86          31 :         if (poFeature)
      87          31 :             return poFeature->Clone();
      88             :     }
      89           1 :     return nullptr;
      90             : }
      91             : 
      92         106 : OGRFeature *OGRILI1Layer::GetNextFeatureRef()
      93             : {
      94         106 :     OGRFeature *poFeature = nullptr;
      95         106 :     if (nFeatureIdx < nFeatures)
      96             :     {
      97          99 :         poFeature = papoFeatures[nFeatureIdx++];
      98             :         // apply filters
      99         198 :         if ((m_poFilterGeom == nullptr ||
     100         198 :              FilterGeometry(poFeature->GetGeometryRef())) &&
     101          99 :             (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
     102          99 :             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          58 : GIntBig OGRILI1Layer::GetFeatureCount(int bForce)
     146             : {
     147          58 :     if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr
     148             :         /* && poAreaLineLayer == NULL*/)
     149             :     {
     150          58 :         return nFeatures;
     151             :     }
     152             :     else
     153             :     {
     154           0 :         return OGRLayer::GetFeatureCount(bForce);
     155             :     }
     156             : }
     157             : 
     158             : /************************************************************************/
     159             : /*                           TestCapability()                           */
     160             : /************************************************************************/
     161             : 
     162           2 : int OGRILI1Layer::TestCapability(CPL_UNUSED const char *pszCap)
     163             : {
     164           2 :     if (EQUAL(pszCap, OLCCurveGeometries))
     165           0 :         return TRUE;
     166           2 :     if (EQUAL(pszCap, OLCZGeometries))
     167           0 :         return TRUE;
     168             : 
     169           2 :     return FALSE;
     170             : }
     171             : 
     172             : /************************************************************************/
     173             : /*                         Internal routines                            */
     174             : /************************************************************************/
     175             : 
     176          17 : void OGRILI1Layer::JoinGeomLayers()
     177             : {
     178          17 :     bGeomsJoined = true;
     179             : 
     180             :     CPLConfigOptionSetter oSetter("OGR_ARC_STEPSIZE", "0.96",
     181          34 :                                   /* bSetOnlyIfUndefined = */ true);
     182             : 
     183          29 :     for (GeomFieldInfos::const_iterator it = oGeomFieldInfos.begin();
     184          41 :          it != oGeomFieldInfos.end(); ++it)
     185             :     {
     186          12 :         OGRFeatureDefn *geomFeatureDefn = it->second.GetGeomTableDefnRef();
     187          12 :         if (geomFeatureDefn)
     188             :         {
     189          14 :             CPLDebug("OGR_ILI", "Join geometry table %s of field '%s'",
     190          14 :                      geomFeatureDefn->GetName(), it->first.c_str());
     191             :             OGRILI1Layer *poGeomLayer =
     192           7 :                 poDS->GetLayerByName(geomFeatureDefn->GetName());
     193             :             const int nGeomFieldIndex =
     194           7 :                 GetLayerDefn()->GetGeomFieldIndex(it->first.c_str());
     195           7 :             if (it->second.iliGeomType == "Surface")
     196             :             {
     197           5 :                 JoinSurfaceLayer(poGeomLayer, nGeomFieldIndex);
     198             :             }
     199           2 :             else if (it->second.iliGeomType == "Area")
     200             :             {
     201           4 :                 CPLString pointField = it->first + "__Point";
     202             :                 const int nPointFieldIndex =
     203           2 :                     GetLayerDefn()->GetGeomFieldIndex(pointField.c_str());
     204           2 :                 PolygonizeAreaLayer(poGeomLayer, nGeomFieldIndex,
     205             :                                     nPointFieldIndex);
     206             :             }
     207             :         }
     208             :     }
     209          17 : }
     210             : 
     211           5 : void OGRILI1Layer::JoinSurfaceLayer(OGRILI1Layer *poSurfaceLineLayer,
     212             :                                     int nSurfaceFieldIndex)
     213             : {
     214           5 :     CPLDebug("OGR_ILI", "Joining surface layer %s with geometries",
     215           5 :              GetLayerDefn()->GetName());
     216             :     OGRwkbGeometryType geomType =
     217           5 :         GetLayerDefn()->GetGeomFieldDefn(nSurfaceFieldIndex)->GetType();
     218             : 
     219          10 :     std::map<OGRFeature *, std::vector<OGRCurve *>> oMapFeatureToGeomSet;
     220             : 
     221           5 :     poSurfaceLineLayer->ResetReading();
     222             : 
     223             :     // First map: for each target curvepolygon, find all belonging curves
     224          24 :     while (OGRFeature *linefeature = poSurfaceLineLayer->GetNextFeatureRef())
     225             :     {
     226             :         // OBJE entries with same _RefTID are polygon rings of same feature
     227          19 :         OGRFeature *feature = nullptr;
     228          19 :         if (poFeatureDefn->GetFieldDefn(0)->GetType() == OFTString)
     229             :         {
     230          19 :             feature = GetFeatureRef(linefeature->GetFieldAsString(1));
     231             :         }
     232             :         else
     233             :         {
     234           0 :             GIntBig reftid = linefeature->GetFieldAsInteger64(1);
     235           0 :             feature = GetFeatureRef(reftid);
     236             :         }
     237          19 :         if (feature)
     238             :         {
     239          19 :             OGRGeometry *poGeom = linefeature->GetGeomFieldRef(0);
     240          19 :             OGRMultiCurve *curves = dynamic_cast<OGRMultiCurve *>(poGeom);
     241          19 :             if (curves)
     242             :             {
     243          38 :                 for (auto &&curve : curves)
     244             :                 {
     245          19 :                     if (!curve->IsEmpty())
     246          19 :                         oMapFeatureToGeomSet[feature].push_back(curve);
     247             :                 }
     248             :             }
     249             :         }
     250             :         else
     251             :         {
     252           0 :             CPLError(CE_Warning, CPLE_AppDefined,
     253             :                      "Couldn't join feature FID " CPL_FRMT_GIB,
     254             :                      linefeature->GetFieldAsInteger64(1));
     255             :         }
     256          19 :     }
     257             : 
     258             :     // Now for each target polygon, assemble the curves together.
     259             :     std::map<OGRFeature *, std::vector<OGRCurve *>>::const_iterator oIter =
     260           5 :         oMapFeatureToGeomSet.begin();
     261          18 :     for (; oIter != oMapFeatureToGeomSet.end(); ++oIter)
     262             :     {
     263          13 :         OGRFeature *feature = oIter->first;
     264          26 :         std::vector<OGRCurve *> oCurves = oIter->second;
     265             : 
     266          26 :         std::vector<OGRCurve *> oSetDestCurves;
     267          13 :         double dfLargestArea = 0.0;
     268          13 :         OGRCurve *poLargestCurve = nullptr;
     269             :         while (true)
     270             :         {
     271          28 :             std::vector<OGRCurve *>::iterator oIterCurves = oCurves.begin();
     272          28 :             if (oIterCurves == oCurves.end())
     273          13 :                 break;
     274             : 
     275          30 :             OGRPoint endPointCC;
     276          15 :             OGRCompoundCurve *poCC = new OGRCompoundCurve();
     277             : 
     278          15 :             bool bFirst = true;
     279             :             while (true)
     280             :             {
     281          19 :                 bool bNewCurveAdded = false;
     282          19 :                 const double dfEps = 1e-14;
     283          20 :                 for (oIterCurves = oCurves.begin();
     284          21 :                      oIterCurves != oCurves.end(); ++oIterCurves)
     285             :                 {
     286          20 :                     OGRCurve *curve = *oIterCurves;
     287          20 :                     OGRPoint startPoint;
     288          20 :                     OGRPoint endPoint;
     289          20 :                     curve->StartPoint(&startPoint);
     290          20 :                     curve->EndPoint(&endPoint);
     291          25 :                     if (bFirst ||
     292           5 :                         (fabs(startPoint.getX() - endPointCC.getX()) < dfEps &&
     293           4 :                          fabs(startPoint.getY() - endPointCC.getY()) < dfEps))
     294             :                     {
     295          19 :                         bFirst = false;
     296             : 
     297          19 :                         curve->EndPoint(&endPointCC);
     298             : 
     299             :                         const OGRwkbGeometryType eCurveType =
     300          19 :                             wkbFlatten(curve->getGeometryType());
     301          19 :                         if (eCurveType == wkbCompoundCurve)
     302             :                         {
     303             :                             OGRCompoundCurve *poCCSub =
     304          11 :                                 curve->toCompoundCurve();
     305          24 :                             for (auto &&subCurve : poCCSub)
     306             :                             {
     307          13 :                                 poCC->addCurve(subCurve);
     308             :                             }
     309             :                         }
     310             :                         else
     311             :                         {
     312           8 :                             poCC->addCurve(curve);
     313             :                         }
     314          19 :                         oCurves.erase(oIterCurves);
     315          19 :                         bNewCurveAdded = true;
     316          19 :                         break;
     317             :                     }
     318           1 :                     else if (fabs(endPoint.getX() - endPointCC.getX()) <
     319           1 :                                  dfEps &&
     320           0 :                              fabs(endPoint.getY() - endPointCC.getY()) < dfEps)
     321             :                     {
     322           0 :                         curve->StartPoint(&endPointCC);
     323             : 
     324             :                         const OGRwkbGeometryType eCurveType =
     325           0 :                             wkbFlatten(curve->getGeometryType());
     326           0 :                         if (eCurveType == wkbLineString ||
     327             :                             eCurveType == wkbCircularString)
     328             :                         {
     329             :                             OGRSimpleCurve *poSC =
     330           0 :                                 curve->toSimpleCurve()->clone();
     331           0 :                             poSC->reversePoints();
     332           0 :                             poCC->addCurveDirectly(poSC);
     333             :                         }
     334           0 :                         else if (eCurveType == wkbCompoundCurve)
     335             :                         {
     336             :                             // Reverse the order of the elements of the
     337             :                             // compound curve
     338             :                             OGRCompoundCurve *poCCSub =
     339           0 :                                 curve->toCompoundCurve();
     340           0 :                             for (int i = poCCSub->getNumCurves() - 1; i >= 0;
     341             :                                  --i)
     342             :                             {
     343             :                                 OGRSimpleCurve *poSC = poCCSub->getCurve(i)
     344           0 :                                                            ->toSimpleCurve()
     345           0 :                                                            ->clone();
     346           0 :                                 poSC->reversePoints();
     347           0 :                                 poCC->addCurveDirectly(poSC);
     348             :                             }
     349             :                         }
     350             : 
     351           0 :                         oCurves.erase(oIterCurves);
     352           0 :                         bNewCurveAdded = true;
     353           0 :                         break;
     354             :                     }
     355             :                 }
     356          19 :                 if (!bNewCurveAdded || oCurves.empty() || poCC->get_IsClosed())
     357          15 :                     break;
     358           4 :             }
     359             : 
     360          15 :             if (!poCC->get_IsClosed())
     361             :             {
     362           0 :                 char *pszJSon = poCC->exportToJson();
     363           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     364             :                          "A ring %s for feature " CPL_FRMT_GIB " in layer %s "
     365             :                          "was not closed. Dropping it",
     366           0 :                          pszJSon, feature->GetFID(), GetName());
     367           0 :                 delete poCC;
     368           0 :                 CPLFree(pszJSon);
     369             :             }
     370             :             else
     371             :             {
     372          15 :                 double dfArea = poCC->get_Area();
     373          15 :                 if (dfArea >= dfLargestArea)
     374             :                 {
     375          15 :                     dfLargestArea = dfArea;
     376          15 :                     poLargestCurve = poCC;
     377             :                 }
     378          15 :                 oSetDestCurves.push_back(poCC);
     379             :             }
     380          15 :         }
     381             : 
     382             :         // Now build the final polygon by first inserting the largest ring.
     383             :         OGRCurvePolygon *poPoly =
     384          13 :             (geomType == wkbPolygon) ? new OGRPolygon() : new OGRCurvePolygon();
     385          13 :         if (poLargestCurve)
     386             :         {
     387             :             std::vector<OGRCurve *>::iterator oIterCurves =
     388          13 :                 oSetDestCurves.begin();
     389          15 :             for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
     390             :             {
     391          15 :                 OGRCurve *poCurve = *oIterCurves;
     392          15 :                 if (poCurve == poLargestCurve)
     393             :                 {
     394          13 :                     oSetDestCurves.erase(oIterCurves);
     395          13 :                     break;
     396             :                 }
     397             :             }
     398             : 
     399          13 :             if (geomType == wkbPolygon)
     400             :             {
     401           6 :                 poLargestCurve = OGRCurve::CastToLinearRing(poLargestCurve);
     402             :             }
     403          13 :             OGRErr error = poPoly->addRingDirectly(poLargestCurve);
     404          13 :             if (error != OGRERR_NONE)
     405             :             {
     406           0 :                 char *pszJSon = poLargestCurve->exportToJson();
     407           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     408             :                          "Cannot add ring %s to feature " CPL_FRMT_GIB
     409             :                          " in layer %s",
     410           0 :                          pszJSon, feature->GetFID(), GetName());
     411           0 :                 CPLFree(pszJSon);
     412             :             }
     413             : 
     414          13 :             oIterCurves = oSetDestCurves.begin();
     415          15 :             for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
     416             :             {
     417           2 :                 OGRCurve *poCurve = *oIterCurves;
     418           2 :                 if (geomType == wkbPolygon)
     419             :                 {
     420           1 :                     poCurve = OGRCurve::CastToLinearRing(poCurve);
     421             :                 }
     422           2 :                 error = poPoly->addRingDirectly(poCurve);
     423           2 :                 if (error != OGRERR_NONE)
     424             :                 {
     425           0 :                     char *pszJSon = poCurve->exportToJson();
     426           0 :                     CPLError(CE_Warning, CPLE_AppDefined,
     427             :                              "Cannot add ring %s to feature " CPL_FRMT_GIB
     428             :                              " in layer %s",
     429           0 :                              pszJSon, feature->GetFID(), GetName());
     430           0 :                     CPLFree(pszJSon);
     431             :                 }
     432             :             }
     433             :         }
     434             : 
     435          13 :         feature->SetGeomFieldDirectly(nSurfaceFieldIndex, poPoly);
     436             :     }
     437             : 
     438           5 :     ResetReading();
     439           5 : }
     440             : 
     441           2 : OGRMultiPolygon *OGRILI1Layer::Polygonize(OGRGeometryCollection *poLines,
     442             :                                           bool
     443             : #if defined(HAVE_GEOS)
     444             :                                               fix_crossing_lines
     445             : #endif
     446             : )
     447             : {
     448           2 :     if (poLines->getNumGeometries() == 0)
     449             :     {
     450           0 :         return new OGRMultiPolygon();
     451             :     }
     452             : 
     453             : #if defined(HAVE_GEOS)
     454           2 :     GEOSGeom *ahInGeoms = nullptr;
     455           2 :     OGRGeometryCollection *poNoncrossingLines = poLines;
     456           2 :     GEOSGeom hResultGeom = nullptr;
     457           2 :     OGRGeometry *poMP = nullptr;
     458             : 
     459           2 :     if (fix_crossing_lines && poLines->getNumGeometries() > 0)
     460             :     {
     461           0 :         CPLDebug("OGR_ILI", "Fixing crossing lines");
     462             :         // A union of the geometry collection with one line fixes
     463             :         // invalid geometries.
     464           0 :         OGRGeometry *poUnion = poLines->Union(poLines->getGeometryRef(0));
     465           0 :         if (poUnion != nullptr)
     466             :         {
     467           0 :             if (wkbFlatten(poUnion->getGeometryType()) ==
     468           0 :                     wkbGeometryCollection ||
     469           0 :                 wkbFlatten(poUnion->getGeometryType()) == wkbMultiLineString)
     470             :             {
     471           0 :                 poNoncrossingLines =
     472           0 :                     dynamic_cast<OGRGeometryCollection *>(poUnion);
     473           0 :                 CPLDebug("OGR_ILI", "Fixed lines: %d",
     474           0 :                          poNoncrossingLines->getNumGeometries() -
     475           0 :                              poLines->getNumGeometries());
     476             :             }
     477             :             else
     478             :             {
     479           0 :                 delete poUnion;
     480             :             }
     481             :         }
     482             :     }
     483             : 
     484           2 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     485             : 
     486             :     ahInGeoms = static_cast<GEOSGeom *>(
     487           2 :         CPLCalloc(sizeof(void *), poNoncrossingLines->getNumGeometries()));
     488          10 :     for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
     489          16 :         ahInGeoms[i] =
     490           8 :             poNoncrossingLines->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
     491             : 
     492           4 :     hResultGeom = GEOSPolygonize_r(hGEOSCtxt, ahInGeoms,
     493           2 :                                    poNoncrossingLines->getNumGeometries());
     494             : 
     495          10 :     for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
     496           8 :         GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
     497           2 :     CPLFree(ahInGeoms);
     498           2 :     if (poNoncrossingLines != poLines)
     499           0 :         delete poNoncrossingLines;
     500             : 
     501           2 :     if (hResultGeom == nullptr)
     502             :     {
     503           0 :         OGRGeometry::freeGEOSContext(hGEOSCtxt);
     504           0 :         return new OGRMultiPolygon();
     505             :     }
     506             : 
     507           2 :     poMP = OGRGeometryFactory::createFromGEOS(hGEOSCtxt, hResultGeom);
     508             : 
     509           2 :     GEOSGeom_destroy_r(hGEOSCtxt, hResultGeom);
     510           2 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     511             : 
     512           2 :     poMP = OGRGeometryFactory::forceToMultiPolygon(poMP);
     513           2 :     if (poMP && wkbFlatten(poMP->getGeometryType()) == wkbMultiPolygon)
     514           2 :         return dynamic_cast<OGRMultiPolygon *>(poMP);
     515             : 
     516           0 :     delete poMP;
     517           0 :     return new OGRMultiPolygon();
     518             : 
     519             : #else
     520             :     return new OGRMultiPolygon();
     521             : #endif
     522             : }
     523             : 
     524           2 : void OGRILI1Layer::PolygonizeAreaLayer(OGRILI1Layer *poAreaLineLayer,
     525             :                                        int
     526             : #if defined(HAVE_GEOS)
     527             :                                            nAreaFieldIndex
     528             : #endif
     529             :                                        ,
     530             :                                        int
     531             : #if defined(HAVE_GEOS)
     532             :                                            nPointFieldIndex
     533             : #endif
     534             : )
     535             : {
     536             :     // add all lines from poAreaLineLayer to collection
     537           2 :     OGRGeometryCollection *gc = new OGRGeometryCollection();
     538           2 :     poAreaLineLayer->ResetReading();
     539          10 :     while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
     540           8 :         gc->addGeometry(feature->GetGeometryRef());
     541             : 
     542             :     // polygonize lines
     543           4 :     CPLDebug("OGR_ILI", "Polygonizing layer %s with %d multilines",
     544           2 :              poAreaLineLayer->GetLayerDefn()->GetName(),
     545             :              gc->getNumGeometries());
     546           2 :     poAreaLineLayer = nullptr;
     547           2 :     OGRMultiPolygon *polys = Polygonize(gc, false);
     548           2 :     CPLDebug("OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
     549           2 :     if (polys->getNumGeometries() != GetFeatureCount())
     550             :     {
     551           0 :         CPLDebug("OGR_ILI", "Feature count of layer %s: " CPL_FRMT_GIB,
     552           0 :                  GetLayerDefn()->GetName(), GetFeatureCount());
     553           0 :         CPLDebug("OGR_ILI", "Polygonizing again with crossing line fix");
     554           0 :         delete polys;
     555           0 :         polys = Polygonize(gc, true);  // try again with crossing line fix
     556           0 :         CPLDebug("OGR_ILI", "Resulting polygons: %d",
     557             :                  polys->getNumGeometries());
     558             :     }
     559           2 :     delete gc;
     560             : 
     561             :     // associate polygon feature with data row according to centroid
     562             : #if defined(HAVE_GEOS)
     563           4 :     OGRPolygon emptyPoly;
     564             : 
     565           2 :     CPLDebug("OGR_ILI", "Associating layer %s with area polygons",
     566           2 :              GetLayerDefn()->GetName());
     567             :     GEOSGeom *ahInGeoms = static_cast<GEOSGeom *>(
     568           2 :         CPLCalloc(sizeof(void *), polys->getNumGeometries()));
     569           2 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     570           8 :     for (int i = 0; i < polys->getNumGeometries(); i++)
     571             :     {
     572           6 :         ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
     573           6 :         if (!GEOSisValid_r(hGEOSCtxt, ahInGeoms[i]))
     574           0 :             ahInGeoms[i] = nullptr;
     575             :     }
     576           8 :     for (int nFidx = 0; nFidx < nFeatures; nFidx++)
     577             :     {
     578           6 :         OGRFeature *feature = papoFeatures[nFidx];
     579           6 :         OGRGeometry *geomRef = feature->GetGeomFieldRef(nPointFieldIndex);
     580           6 :         if (!geomRef)
     581             :         {
     582           0 :             continue;
     583             :         }
     584             :         GEOSGeom point =
     585           6 :             reinterpret_cast<GEOSGeom>(geomRef->exportToGEOS(hGEOSCtxt));
     586             : 
     587           6 :         int i = 0;
     588          12 :         for (; i < polys->getNumGeometries(); i++)
     589             :         {
     590          12 :             if (ahInGeoms[i] && GEOSWithin_r(hGEOSCtxt, point, ahInGeoms[i]))
     591             :             {
     592           6 :                 feature->SetGeomField(nAreaFieldIndex,
     593           6 :                                       polys->getGeometryRef(i));
     594           6 :                 break;
     595             :             }
     596             :         }
     597           6 :         if (i == polys->getNumGeometries())
     598             :         {
     599           0 :             CPLDebug("OGR_ILI", "Association between area and point failed.");
     600           0 :             feature->SetGeometry(&emptyPoly);
     601             :         }
     602           6 :         GEOSGeom_destroy_r(hGEOSCtxt, point);
     603             :     }
     604           8 :     for (int i = 0; i < polys->getNumGeometries(); i++)
     605           6 :         GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
     606           2 :     CPLFree(ahInGeoms);
     607           2 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     608             : #endif
     609           2 :     delete polys;
     610           2 : }
     611             : 
     612             : /************************************************************************/
     613             : /*                             GetDataset()                             */
     614             : /************************************************************************/
     615             : 
     616           1 : GDALDataset *OGRILI1Layer::GetDataset()
     617             : {
     618           1 :     return poDS;
     619             : }