diff --git a/autotest/ogr/ogr_vfk.py b/autotest/ogr/ogr_vfk.py index f7d316cbbbcd..d7f510883270 100755 --- a/autotest/ogr/ogr_vfk.py +++ b/autotest/ogr/ogr_vfk.py @@ -70,6 +70,7 @@ def test_ogr_vfk_1(vfk_ds): # Read the first feature from layer 'PAR', check envelope +@pytest.mark.require_geos def test_ogr_vfk_2(vfk_ds): vfk_layer_par = vfk_ds.GetLayer(0) diff --git a/ogr/ogr_api.h b/ogr/ogr_api.h index 94d1f2c6825d..1a854dc2c6f9 100644 --- a/ogr/ogr_api.h +++ b/ogr/ogr_api.h @@ -258,6 +258,7 @@ OGRGeometryH CPL_DLL OGR_G_SetPrecision(OGRGeometryH, double dfGridSize, int nFlags) CPL_WARN_UNUSED_RESULT; OGRGeometryH CPL_DLL OGR_G_Polygonize(OGRGeometryH) CPL_WARN_UNUSED_RESULT; +OGRGeometryH CPL_DLL OGR_G_BuildArea(OGRGeometryH) CPL_WARN_UNUSED_RESULT; /*! @cond Doxygen_Suppress */ /* backward compatibility (non-standard methods) */ diff --git a/ogr/ogr_geometry.h b/ogr/ogr_geometry.h index b64c98f6ff88..e2227f3443f0 100644 --- a/ogr/ogr_geometry.h +++ b/ogr/ogr_geometry.h @@ -589,6 +589,7 @@ class CPL_DLL OGRGeometry int bOnlyEdges) const CPL_WARN_UNUSED_RESULT; virtual OGRGeometry *Polygonize() const CPL_WARN_UNUSED_RESULT; + virtual OGRGeometry *BuildArea() const CPL_WARN_UNUSED_RESULT; virtual double Distance3D(const OGRGeometry *poOtherGeom) const; diff --git a/ogr/ogrgeometry.cpp b/ogr/ogrgeometry.cpp index f429bee27211..f209396e5467 100644 --- a/ogr/ogrgeometry.cpp +++ b/ogr/ogrgeometry.cpp @@ -6885,6 +6885,100 @@ OGRGeometryH OGR_G_Polygonize(OGRGeometryH hTarget) OGRGeometry::FromHandle(hTarget)->Polygonize()); } +/************************************************************************/ +/* BuildArea() */ +/************************************************************************/ + +/** + * \brief Polygonize a linework assuming inner polygons are holes. + * + * This method is the same as the C function OGR_G_BuildArea(). + * + * Polygonization is performed similarly to OGRGeometry::Polygonize(). + * Additionally, holes are dropped and the result is unified producing + * a single Polygon or a MultiPolygon. + * + * A new geometry object is created and returned: NULL on failure, + * empty GeometryCollection if the input geometry cannot be polygonized, + * Polygon or MultiPolygon on success. + * + * This method is built on the GEOSBuildArea_r() function of the GEOS + * library, check it for the definition of the geometry operation. + * If OGR is built without the GEOS library, this method will always fail, + * issuing a CPLE_NotSupported error. + * + * @return a newly allocated geometry now owned by the caller, + * or NULL on failure. + * + * @since OGR 3.11 + */ + +OGRGeometry *OGRGeometry::BuildArea() const + +{ +#ifndef HAVE_GEOS + + CPLError(CE_Failure, CPLE_NotSupported, "GEOS support not enabled."); + return nullptr; + +#else + + OGRGeometry *poPolygsOGRGeom = nullptr; + + GEOSContextHandle_t hGEOSCtxt = createGEOSContext(); + GEOSGeom hThisGeosGeom = exportToGEOS(hGEOSCtxt); + if (hThisGeosGeom != nullptr) + { + GEOSGeom hGeosPolygs = GEOSBuildArea_r(hGEOSCtxt, hThisGeosGeom); + poPolygsOGRGeom = + BuildGeometryFromGEOS(hGEOSCtxt, hGeosPolygs, this, nullptr); + GEOSGeom_destroy_r(hGEOSCtxt, hThisGeosGeom); + } + freeGEOSContext(hGEOSCtxt); + + return poPolygsOGRGeom; + +#endif // HAVE_GEOS +} + +/************************************************************************/ +/* OGR_G_BuildArea() */ +/************************************************************************/ + +/** + * \brief Polygonize a linework assuming inner polygons are holes. + * + * This function is the same as the C++ method OGRGeometry::BuildArea(). + * + * Polygonization is performed similarly to OGR_G_Polygonize(). + * Additionally, holes are dropped and the result is unified producing + * a single Polygon or a MultiPolygon. + * + * A new geometry object is created and returned: NULL on failure, + * empty GeometryCollection if the input geometry cannot be polygonized, + * Polygon or MultiPolygon on success. + * + * This function is built on the GEOSBuildArea_r() function of the GEOS + * library, check it for the definition of the geometry operation. + * If OGR is built without the GEOS library, this function will always fail, + * issuing a CPLE_NotSupported error. + * + * @param hGeom handle on the geometry to polygonize. + * + * @return a handle on newly allocated geometry now owned by the caller, + * or NULL on failure. + * + * @since OGR 3.11 + */ + +OGRGeometryH OGR_G_BuildArea(OGRGeometryH hGeom) + +{ + VALIDATE_POINTER1(hGeom, "OGR_G_BuildArea", nullptr); + + return OGRGeometry::ToHandle(OGRGeometry::FromHandle(hGeom)->BuildArea()); +} + /************************************************************************/ /* swapXY() */ /************************************************************************/ diff --git a/ogr/ogrsf_frmts/vfk/CMakeLists.txt b/ogr/ogrsf_frmts/vfk/CMakeLists.txt index ff3181f0a0d6..2d1045a78781 100644 --- a/ogr/ogrsf_frmts/vfk/CMakeLists.txt +++ b/ogr/ogrsf_frmts/vfk/CMakeLists.txt @@ -29,3 +29,7 @@ if (GDAL_USE_SQLITE3) gdal_target_link_libraries(ogr_VFK PRIVATE SQLite::SQLite3) target_compile_definitions(ogr_VFK PRIVATE -DHAVE_SQLITE) endif () +if (GDAL_USE_GEOS) + target_compile_definitions(ogr_VFK PRIVATE -DHAVE_GEOS=1) + gdal_target_link_libraries(ogr_VFK PRIVATE ${GEOS_TARGET}) +endif () diff --git a/ogr/ogrsf_frmts/vfk/vfkdatablocksqlite.cpp b/ogr/ogrsf_frmts/vfk/vfkdatablocksqlite.cpp index 479a24a197d7..ed7131e7cc08 100644 --- a/ogr/ogrsf_frmts/vfk/vfkdatablocksqlite.cpp +++ b/ogr/ogrsf_frmts/vfk/vfkdatablocksqlite.cpp @@ -452,6 +452,15 @@ int VFKDataBlockSQLite::LoadGeometryLineStringHP() */ int VFKDataBlockSQLite::LoadGeometryPolygon() { +#ifndef HAVE_GEOS + + CPLError(CE_Warning, CPLE_NotSupported, + "GEOS support not enabled. Unable to build geometry for %s.", + m_pszName); + return -1; + +#else + VFKReaderSQLite *poReader = (VFKReaderSQLite *)m_poReader; VFKDataBlockSQLite *poDataBlockLines1 = nullptr; @@ -514,11 +523,6 @@ int VFKDataBlockSQLite::LoadGeometryPolygon() if (poReader->IsSpatial()) poReader->ExecuteSQL("BEGIN"); - VFKFeatureSQLiteList poLineList; - /* first is to be considered as exterior */ - PointListArray poRingList; - std::vector poLinearRingList; - OGRPolygon ogrPolygon; int nInvalidNoLines = 0; int nInvalidNoRings = 0; int nGeometries = 0; @@ -537,6 +541,8 @@ int VFKDataBlockSQLite::LoadGeometryPolygon() continue; } + /* collect boundary lines */ + VFKFeatureSQLiteList poLineList; if (bIsPar) { vrValue[0] = vrValue[1] = id; @@ -544,8 +550,6 @@ int VFKDataBlockSQLite::LoadGeometryPolygon() } else { - // std::vector poLineListOb; - osSQL.Printf("SELECT ID FROM %s WHERE BUD_ID = " CPL_FRMT_GUIB, poDataBlockLines1->GetName(), id); if (poReader->IsSpatial()) @@ -567,236 +571,40 @@ int VFKDataBlockSQLite::LoadGeometryPolygon() poLineList.push_back(poLineSbp); } } - size_t nLines = poLineList.size(); - if (nLines < 1) + if (poLineList.empty()) { - CPLDebug( - "OGR-VFK", - "%s: unable to collect rings for polygon fid = %ld (no lines)", - m_pszName, iFID); + CPLDebug("OGR-VFK", "%s: no lines to polygonize (fid = %ld)", + m_pszName, iFID); nInvalidNoLines++; continue; } - /* clear */ - ogrPolygon.empty(); - - /* free ring list */ - for (PointListArray::iterator iRing = poRingList.begin(), - eRing = poRingList.end(); - iRing != eRing; ++iRing) - { - delete (*iRing); - *iRing = nullptr; - } - poRingList.clear(); - - /* collect rings from lines */ -#if 1 - // Fast version using a map to identify quickly a ring from its end - // point. - std::map, PointList *> oMapEndRing; - while (!poLineList.empty()) - { - auto pGeom = poLineList.front()->GetGeometry(); - if (pGeom) - { - auto poLine = pGeom->toLineString(); - if (poLine == nullptr || poLine->getNumPoints() < 2) - continue; - poLineList.erase(poLineList.begin()); - PointList *poList = new PointList(); - FillPointList(poList, poLine); - poRingList.emplace_back(poList); - OGRPoint oFirst, oLast; - poLine->StartPoint(&oFirst); - poLine->EndPoint(&oLast); - oMapEndRing[std::pair(oLast.getX(), - oLast.getY())] = poList; - - bool bWorkDone = true; - while (bWorkDone && (*poList).front() != (*poList).back()) - { - bWorkDone = false; - for (auto oIter = poLineList.begin(); - oIter != poLineList.end(); ++oIter) - { - const auto &oCandidate = *oIter; - auto poCandidateGeom = oCandidate->GetGeometry(); - if (poCandidateGeom == nullptr) - continue; - poLine = poCandidateGeom->toLineString(); - if (poLine == nullptr || poLine->getNumPoints() < 2) - continue; - poLine->StartPoint(&oFirst); - poLine->EndPoint(&oLast); - // MER = MapEndRing - auto oIterMER = - oMapEndRing.find(std::pair( - oFirst.getX(), oFirst.getY())); - if (oIterMER != oMapEndRing.end()) - { - auto ring = oIterMER->second; - PointList oList; - FillPointList(&oList, poLine); - /* forward, skip first point */ - ring->insert(ring->end(), oList.begin() + 1, - oList.end()); - poLineList.erase(oIter); - oMapEndRing.erase(oIterMER); - oMapEndRing[std::pair( - oLast.getX(), oLast.getY())] = poList; - bWorkDone = true; - break; - } - oIterMER = oMapEndRing.find(std::pair( - oLast.getX(), oLast.getY())); - if (oIterMER != oMapEndRing.end()) - { - auto ring = oIterMER->second; - PointList oList; - FillPointList(&oList, poLine); - /* backward, skip first point */ - ring->insert(ring->end(), oList.rbegin() + 1, - oList.rend()); - poLineList.erase(oIter); - oMapEndRing.erase(oIterMER); - oMapEndRing[std::pair( - oFirst.getX(), oFirst.getY())] = ring; - bWorkDone = true; - break; - } - } - } - } - } -#else - bool bFound = false; - int nCount = 0; - const int nCountMax = static_cast(nLines) * 2; - while (!poLineList.empty() && nCount < nCountMax) + OGRMultiLineString oMultiLine; + for (VFKFeatureSQLite *poLineFeature : poLineList) { - bool bNewRing = !bFound; - bFound = false; - int i = 1; - for (VFKFeatureSQLiteList::iterator iHp = poLineList.begin(), - eHp = poLineList.end(); - iHp != eHp; ++iHp, ++i) + const OGRGeometry *poLineGeom = poLineFeature->GetGeometry(); + if (poLineGeom) { - auto pGeom = (*iHp)->GetGeometry(); - if (pGeom && AppendLineToRing(&poRingList, - pGeom->toLineString(), bNewRing)) - { - bFound = true; - poLineList.erase(iHp); - break; - } + oMultiLine.addGeometry(poLineGeom); } - nCount++; } -#endif - CPLDebug("OGR-VFK", "%s: fid = %ld nlines = %d -> nrings = %d", - m_pszName, iFID, (int)nLines, (int)poRingList.size()); - if (!poLineList.empty()) + /* polygonize using GEOSBuildArea() */ + auto poPolygonGeom = + std::unique_ptr(oMultiLine.BuildArea()); + /* only Polygons are allowed in VFK, in particular, no MultiPolygons */ + if (!poPolygonGeom || poPolygonGeom->IsEmpty() || + wkbFlatten(poPolygonGeom->getGeometryType()) != wkbPolygon) { - CPLDebug("OGR-VFK", - "%s: unable to collect rings for polygon fid = %ld", + CPLDebug("OGR-VFK", "%s: unable to polygonize (fid = %ld)", m_pszName, iFID); nInvalidNoRings++; continue; } - /* build rings */ - poLinearRingList.clear(); - OGRLinearRing *poOgrRing = nullptr; - for (PointListArray::const_iterator iRing = poRingList.begin(), - eRing = poRingList.end(); - iRing != eRing; ++iRing) - { - PointList *poList = *iRing; - - poLinearRingList.push_back(new OGRLinearRing()); - poOgrRing = poLinearRingList.back(); - CPLAssert(nullptr != poOgrRing); - - for (PointList::iterator iPoint = poList->begin(), - ePoint = poList->end(); - iPoint != ePoint; ++iPoint) - { - OGRPoint *poPoint = &(*iPoint); - poOgrRing->addPoint(poPoint); - } - } - - /* find exterior ring */ - if (poLinearRingList.size() > 1) - { - std::vector::iterator exteriorRing; - - exteriorRing = poLinearRingList.begin(); - double dMaxArea = -1.0; - for (std::vector::iterator - iRing = poLinearRingList.begin(), - eRing = poLinearRingList.end(); - iRing != eRing; ++iRing) - { - poOgrRing = *iRing; - if (!IsRingClosed(poOgrRing)) - continue; /* skip unclosed rings */ - - const double dArea = poOgrRing->get_Area(); - if (dArea > dMaxArea) - { - dMaxArea = dArea; - exteriorRing = iRing; - } - } - if (exteriorRing != poLinearRingList.begin()) - { - std::swap(*poLinearRingList.begin(), *exteriorRing); - } - } - - /* build polygon from rings */ - int nBridges = 0; - for (std::vector::iterator - iRing = poLinearRingList.begin(), - eRing = poLinearRingList.end(); - iRing != eRing; ++iRing) - { - poOgrRing = *iRing; - - /* check if ring is closed */ - if (IsRingClosed(poOgrRing)) - { - ogrPolygon.addRing(poOgrRing); - } - else - { - if (poOgrRing->getNumPoints() == 2) - { - CPLDebug("OGR-VFK", - "%s: Polygon (fid = %ld) bridge removed", - m_pszName, iFID); - nBridges++; - } - else - { - CPLDebug("OGR-VFK", - "%s: Polygon (fid = %ld) unclosed ring skipped", - m_pszName, iFID); - } - } - delete poOgrRing; - *iRing = nullptr; - } - /* set polygon */ - ogrPolygon.setCoordinateDimension(2); /* force 2D */ - if (ogrPolygon.getNumInteriorRings() + nBridges != - (int)poLinearRingList.size() - 1 || - !poFeature->SetGeometry(&ogrPolygon)) + poPolygonGeom->setCoordinateDimension(2); /* force 2D */ + if (!poFeature->SetGeometry(poPolygonGeom.get())) { nInvalidNoRings++; continue; @@ -804,19 +612,10 @@ int VFKDataBlockSQLite::LoadGeometryPolygon() /* store also geometry in DB */ if (poReader->IsSpatial() && - SaveGeometryToDB(&ogrPolygon, rowId) != OGRERR_FAILURE) + SaveGeometryToDB(poPolygonGeom.get(), rowId) != OGRERR_FAILURE) nGeometries++; } - /* free ring list */ - for (PointListArray::iterator iRing = poRingList.begin(), - eRing = poRingList.end(); - iRing != eRing; ++iRing) - { - delete (*iRing); - *iRing = nullptr; - } - CPLDebug("OGR-VFK", "%s: nolines = %d norings = %d", m_pszName, nInvalidNoLines, nInvalidNoRings); @@ -827,6 +626,7 @@ int VFKDataBlockSQLite::LoadGeometryPolygon() poReader->ExecuteSQL("COMMIT"); return nInvalidNoLines + nInvalidNoRings; +#endif // HAVE_GEOS } /*! diff --git a/swig/include/ogr.i b/swig/include/ogr.i index 32ec24a34393..5054c3128d37 100644 --- a/swig/include/ogr.i +++ b/swig/include/ogr.i @@ -3758,6 +3758,11 @@ public: return (OGRGeometryShadow*) OGR_G_Polygonize(self); } + %newobject BuildArea; + OGRGeometryShadow* BuildArea() { + return (OGRGeometryShadow*) OGR_G_BuildArea(self); + } + %newobject Boundary; OGRGeometryShadow* Boundary() { return (OGRGeometryShadow*) OGR_G_Boundary(self); diff --git a/swig/include/python/docs/ogr_geometry_docs.i b/swig/include/python/docs/ogr_geometry_docs.i index 00a9dff17493..c5c923b67571 100644 --- a/swig/include/python/docs/ogr_geometry_docs.i +++ b/swig/include/python/docs/ogr_geometry_docs.i @@ -874,6 +874,17 @@ Geometry: A new geometry or None on failure. "; +%feature("docstring") BuildArea " +Polygonize a linework assuming inner polygons are holes. + +For more details: :cpp:func:`OGR_G_BuildArea` + +Returns +-------- +Geometry: + A new geometry or None on failure. +"; + %feature("docstring") SwapXY " Swap x and y coordinates.