From fa4c932d2a14ab861b1ab15ad7ea6cdaf7d79bf0 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 16:51:12 +0200 Subject: [PATCH 01/10] ogr2ogr: implement GetInverse() in CompositeCT and AxisMappingCoordinateTransformation --- apps/ogr2ogr_lib.cpp | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index 16b3539ddd42..42dea3bce9dc 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -1205,6 +1206,14 @@ class GCPCoordTransformation : public OGRCoordinateTransformation virtual OGRCoordinateTransformation *GetInverse() const override { + static std::once_flag flag; + std::call_once(flag, + []() + { + CPLDebug("OGR2OGR", + "GCPCoordTransformation::GetInverse() " + "called, but not implemented"); + }); return nullptr; } }; @@ -1292,7 +1301,21 @@ class CompositeCT : public OGRCoordinateTransformation virtual OGRCoordinateTransformation *GetInverse() const override { - return nullptr; + if (!poCT1 && !poCT2) + return nullptr; + if (!poCT2) + return poCT1->GetInverse(); + if (!poCT1) + return poCT2->GetInverse(); + auto poInvCT1 = + std::unique_ptr(poCT1->GetInverse()); + auto poInvCT2 = + std::unique_ptr(poCT2->GetInverse()); + if (!poInvCT1 || !poInvCT2) + return nullptr; + return std::make_unique(poInvCT2.release(), true, + poInvCT1.release(), true) + .release(); } }; @@ -1302,9 +1325,14 @@ class CompositeCT : public OGRCoordinateTransformation class AxisMappingCoordinateTransformation : public OGRCoordinateTransformation { - public: bool bSwapXY = false; + explicit AxisMappingCoordinateTransformation(bool bSwapXYIn) + : bSwapXY(bSwapXYIn) + { + } + + public: AxisMappingCoordinateTransformation(const std::vector &mappingIn, const std::vector &mappingOut) { @@ -1360,7 +1388,7 @@ class AxisMappingCoordinateTransformation : public OGRCoordinateTransformation virtual OGRCoordinateTransformation *GetInverse() const override { - return nullptr; + return new AxisMappingCoordinateTransformation(bSwapXY); } }; From ffb6004c25503d80d00575524c44c4adb8bccadd Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 15:40:08 +0200 Subject: [PATCH 02/10] OGRGeometryFactory::transformWithOptions(): deal with polar or anti-meridian discontinuities when going from projected to (any) geographic CRS The logic already existed but was restricted to WGS 84 without any good reason. --- autotest/cpp/test_ogr.cpp | 69 ++++++++++++++++++++ ogr/ogrgeometryfactory.cpp | 125 ++++++++++++++++++++++++------------- 2 files changed, 150 insertions(+), 44 deletions(-) diff --git a/autotest/cpp/test_ogr.cpp b/autotest/cpp/test_ogr.cpp index e75d0b6253b9..8fe4e7469977 100644 --- a/autotest/cpp/test_ogr.cpp +++ b/autotest/cpp/test_ogr.cpp @@ -4201,4 +4201,73 @@ TEST_F(test_ogr, OGRCurve_reversePoints) } } +// Test OGRGeometryFactory::transformWithOptions() +TEST_F(test_ogr, transformWithOptions) +{ + // Projected CRS to national geographic CRS (not including poles or antimeridian) + { + OGRGeometry *poGeom = nullptr; + OGRGeometryFactory::createFromWkt( + "LINESTRING(700000 6600000, 700001 6600001)", nullptr, &poGeom); + ASSERT_NE(poGeom, nullptr); + + OGRSpatialReference oEPSG_2154; + oEPSG_2154.importFromEPSG(2154); // "RGF93 v1 / Lambert-93" + OGRSpatialReference oEPSG_4171; + oEPSG_4171.importFromEPSG(4171); // "RGF93 v1" + oEPSG_4171.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + auto poCT = std::unique_ptr( + OGRCreateCoordinateTransformation(&oEPSG_2154, &oEPSG_4171)); + OGRGeometryFactory::TransformWithOptionsCache oCache; + poGeom = OGRGeometryFactory::transformWithOptions(poGeom, poCT.get(), + nullptr, oCache); + EXPECT_NEAR(poGeom->toLineString()->getX(0), 3, 1e-8); + EXPECT_NEAR(poGeom->toLineString()->getY(0), 46.5, 1e-8); + + delete poGeom; + } + +#ifdef HAVE_GEOS + // Projected CRS to national geographic CRS including antimeridian + { + OGRGeometry *poGeom = nullptr; + OGRGeometryFactory::createFromWkt( + "LINESTRING(657630.64 4984896.17,815261.43 4990738.26)", nullptr, + &poGeom); + ASSERT_NE(poGeom, nullptr); + + OGRSpatialReference oEPSG_6329; + oEPSG_6329.importFromEPSG(6329); // "NAD83(2011) / UTM zone 60N" + OGRSpatialReference oEPSG_6318; + oEPSG_6318.importFromEPSG(6318); // "NAD83(2011)" + oEPSG_6318.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + auto poCT = std::unique_ptr( + OGRCreateCoordinateTransformation(&oEPSG_6329, &oEPSG_6318)); + OGRGeometryFactory::TransformWithOptionsCache oCache; + poGeom = OGRGeometryFactory::transformWithOptions(poGeom, poCT.get(), + nullptr, oCache); + EXPECT_EQ(poGeom->getGeometryType(), wkbMultiLineString); + if (poGeom->getGeometryType() == wkbMultiLineString) + { + const auto poMLS = poGeom->toMultiLineString(); + EXPECT_EQ(poMLS->getNumGeometries(), 2); + if (poMLS->getNumGeometries() == 2) + { + const auto poLS = poMLS->getGeometryRef(0); + EXPECT_EQ(poLS->getNumPoints(), 2); + if (poLS->getNumPoints() == 2) + { + EXPECT_NEAR(poLS->getX(0), 179, 1e-6); + EXPECT_NEAR(poLS->getY(0), 45, 1e-6); + EXPECT_NEAR(poLS->getX(1), 180, 1e-6); + EXPECT_NEAR(poLS->getY(1), 45.004384301691303, 1e-6); + } + } + } + + delete poGeom; + } +#endif +} + } // namespace diff --git a/ogr/ogrgeometryfactory.cpp b/ogr/ogrgeometryfactory.cpp index 90583391997e..532028a25507 100644 --- a/ogr/ogrgeometryfactory.cpp +++ b/ogr/ogrgeometryfactory.cpp @@ -3302,15 +3302,15 @@ static void AlterPole(OGRGeometry *poGeom, OGRPoint *poPole, } /************************************************************************/ -/* IsPolarToWGS84() */ +/* IsPolarToGeographic() */ /* */ /* Returns true if poCT transforms from a projection that includes one */ /* of the pole in a continuous way. */ /************************************************************************/ -static bool IsPolarToWGS84(OGRCoordinateTransformation *poCT, - OGRCoordinateTransformation *poRevCT, - bool &bIsNorthPolarOut) +static bool IsPolarToGeographic(OGRCoordinateTransformation *poCT, + OGRCoordinateTransformation *poRevCT, + bool &bIsNorthPolarOut) { bool bIsNorthPolar = false; bool bIsSouthPolar = false; @@ -3366,14 +3366,14 @@ static bool IsPolarToWGS84(OGRCoordinateTransformation *poCT, } /************************************************************************/ -/* TransformBeforePolarToWGS84() */ +/* TransformBeforePolarToGeographic() */ /* */ /* Transform the geometry (by intersection), so as to cut each geometry */ /* that crosses the pole, in 2 parts. Do also tricks for geometries */ /* that just touch the pole. */ /************************************************************************/ -static std::unique_ptr TransformBeforePolarToWGS84( +static std::unique_ptr TransformBeforePolarToGeographic( OGRCoordinateTransformation *poRevCT, bool bIsNorthPolar, std::unique_ptr poDstGeom, bool &bNeedPostCorrectionOut) { @@ -3449,15 +3449,15 @@ static std::unique_ptr TransformBeforePolarToWGS84( } /************************************************************************/ -/* IsAntimeridianProjToWGS84() */ +/* IsAntimeridianProjToGeographic() */ /* */ /* Returns true if poCT transforms from a projection that includes the */ /* antimeridian in a continuous way. */ /************************************************************************/ -static bool IsAntimeridianProjToWGS84(OGRCoordinateTransformation *poCT, - OGRCoordinateTransformation *poRevCT, - OGRGeometry *poDstGeometry) +static bool IsAntimeridianProjToGeographic(OGRCoordinateTransformation *poCT, + OGRCoordinateTransformation *poRevCT, + OGRGeometry *poDstGeometry) { const bool bBackupEmitErrors = poCT->GetEmitErrors(); poRevCT->SetEmitErrors(false); @@ -3633,13 +3633,13 @@ struct SortPointsByAscendingY }; /************************************************************************/ -/* TransformBeforeAntimeridianToWGS84() */ +/* TransformBeforeAntimeridianToGeographic() */ /* */ /* Transform the geometry (by intersection), so as to cut each geometry */ /* that crosses the antimeridian, in 2 parts. */ /************************************************************************/ -static std::unique_ptr TransformBeforeAntimeridianToWGS84( +static std::unique_ptr TransformBeforeAntimeridianToGeographic( OGRCoordinateTransformation *poCT, OGRCoordinateTransformation *poRevCT, std::unique_ptr poDstGeom, bool &bNeedPostCorrectionOut) { @@ -3820,9 +3820,22 @@ static void SnapCoordsCloseToLatLongBounds(OGRGeometry *poGeom) struct OGRGeometryFactory::TransformWithOptionsCache::Private { + const OGRSpatialReference *poSourceCRS = nullptr; + const OGRSpatialReference *poTargetCRS = nullptr; + const OGRCoordinateTransformation *poCT = nullptr; std::unique_ptr poRevCT{}; bool bIsPolar = false; bool bIsNorthPolar = false; + + void clear() + { + poSourceCRS = nullptr; + poTargetCRS = nullptr; + poCT = nullptr; + poRevCT.reset(); + bIsPolar = false; + bIsNorthPolar = false; + } }; /************************************************************************/ @@ -3858,52 +3871,76 @@ OGRGeometry *OGRGeometryFactory::transformWithOptions( char **papszOptions, CPL_UNUSED const TransformWithOptionsCache &cache) { auto poDstGeom = std::unique_ptr(poSrcGeom->clone()); - if (poCT != nullptr) + if (poCT) { #ifdef HAVE_GEOS bool bNeedPostCorrection = false; - - auto poSourceCRS = poCT->GetSourceCS(); - auto poTargetCRS = poCT->GetTargetCS(); - if (poSourceCRS != nullptr && poTargetCRS != nullptr && - poSourceCRS->IsProjected() && poTargetCRS->IsGeographic()) + const auto poSourceCRS = poCT->GetSourceCS(); + const auto poTargetCRS = poCT->GetTargetCS(); + const auto eSrcGeomType = wkbFlatten(poSrcGeom->getGeometryType()); + // Check if we are transforming from projected coordinates to + // geographic coordinates, with a chance that there might be polar or + // anti-meridian discontinuities. If so, create the inverse transform. + if (eSrcGeomType != wkbPoint && eSrcGeomType != wkbMultiPoint && + (poSourceCRS != cache.d->poSourceCRS || + poTargetCRS != cache.d->poTargetCRS || poCT != cache.d->poCT)) { - OGRSpatialReference oSRSWGS84; - oSRSWGS84.SetWellKnownGeogCS("WGS84"); - oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); - if (poTargetCRS->IsSame(&oSRSWGS84)) + cache.d->clear(); + cache.d->poSourceCRS = poSourceCRS; + cache.d->poTargetCRS = poTargetCRS; + cache.d->poCT = poCT; + if (poSourceCRS && poTargetCRS && poSourceCRS->IsProjected() && + poTargetCRS->IsGeographic() && + poTargetCRS->GetAxisMappingStrategy() == + OAMS_TRADITIONAL_GIS_ORDER && + // check that angular units is degree + std::fabs(poTargetCRS->GetAngularUnits(nullptr) - + CPLAtof(SRS_UA_DEGREE_CONV)) <= + 1e-8 * CPLAtof(SRS_UA_DEGREE_CONV)) { - if (cache.d->poRevCT == nullptr || - !cache.d->poRevCT->GetTargetCS()->IsSame(poSourceCRS)) + double dfWestLong = 0.0; + double dfSouthLat = 0.0; + double dfEastLong = 0.0; + double dfNorthLat = 0.0; + if (poSourceCRS->GetAreaOfUse(&dfWestLong, &dfSouthLat, + &dfEastLong, &dfNorthLat, + nullptr) && + !(dfSouthLat == -90.0 || dfNorthLat == 90.0 || + dfWestLong == -180.0 || dfEastLong == 180.0 || + dfWestLong > dfEastLong)) + { + // Not a global geographic CRS + } + else { cache.d->poRevCT.reset(OGRCreateCoordinateTransformation( - &oSRSWGS84, poSourceCRS)); + poTargetCRS, poSourceCRS)); cache.d->bIsNorthPolar = false; cache.d->bIsPolar = false; + cache.d->poRevCT.reset(poCT->GetInverse()); if (cache.d->poRevCT && - IsPolarToWGS84(poCT, cache.d->poRevCT.get(), - cache.d->bIsNorthPolar)) + IsPolarToGeographic(poCT, cache.d->poRevCT.get(), + cache.d->bIsNorthPolar)) { cache.d->bIsPolar = true; } } - auto poRevCT = cache.d->poRevCT.get(); - if (poRevCT != nullptr) - { - if (cache.d->bIsPolar) - { - poDstGeom = TransformBeforePolarToWGS84( - poRevCT, cache.d->bIsNorthPolar, - std::move(poDstGeom), bNeedPostCorrection); - } - else if (IsAntimeridianProjToWGS84(poCT, poRevCT, - poDstGeom.get())) - { - poDstGeom = TransformBeforeAntimeridianToWGS84( - poCT, poRevCT, std::move(poDstGeom), - bNeedPostCorrection); - } - } + } + } + + if (auto poRevCT = cache.d->poRevCT.get()) + { + if (cache.d->bIsPolar) + { + poDstGeom = TransformBeforePolarToGeographic( + poRevCT, cache.d->bIsNorthPolar, std::move(poDstGeom), + bNeedPostCorrection); + } + else if (IsAntimeridianProjToGeographic(poCT, poRevCT, + poDstGeom.get())) + { + poDstGeom = TransformBeforeAntimeridianToGeographic( + poCT, poRevCT, std::move(poDstGeom), bNeedPostCorrection); } } #endif From 1e6c00c468f8f1134215ac041d535eeec6a0ee86 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 20:21:59 +0200 Subject: [PATCH 03/10] OGRCoordinateTransformation::TransformWithErrorCodes(): optimize when nCount == 1 --- ogr/ogrct.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/ogr/ogrct.cpp b/ogr/ogrct.cpp index 49b83af8b096..4c62b1dc091b 100644 --- a/ogr/ogrct.cpp +++ b/ogr/ogrct.cpp @@ -2232,6 +2232,18 @@ int OGRCoordinateTransformation::TransformWithErrorCodes(size_t nCount, int *panErrorCodes) { + if (nCount == 1) + { + int nSuccess = 0; + const bool bOverallSuccess = + CPL_TO_BOOL(Transform(nCount, x, y, z, t, &nSuccess)); + if (panErrorCodes) + { + panErrorCodes[0] = nSuccess ? 0 : -1; + } + return bOverallSuccess; + } + std::vector abSuccess; try { From cfc5ff0c37fdfc0ce28527f4669b6d89d61f2120 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 20:28:40 +0200 Subject: [PATCH 04/10] OGRProjCT::TransformWithErrorCodes(): speed-up by avoiding OSRGetProjTLSContext() when possible --- ogr/ogrct.cpp | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/ogr/ogrct.cpp b/ogr/ogrct.cpp index 4c62b1dc091b..afc3bbb489ea 100644 --- a/ogr/ogrct.cpp +++ b/ogr/ogrct.cpp @@ -37,6 +37,7 @@ #include #include #include +#include #include "cpl_conv.h" #include "cpl_error.h" @@ -765,6 +766,9 @@ class OGRProjCT : public OGRCoordinateTransformation double dfThreshold = 0.0; + PJ_CONTEXT *m_psLastContext = nullptr; + std::thread::id m_nLastContextThreadId{}; + PjPtr m_pj{}; bool m_bReversePj = false; @@ -1277,9 +1281,10 @@ OGRProjCT::OGRProjCT(const OGRProjCT &other) m_osTargetSRS(other.m_osTargetSRS), bWebMercatorToWGS84LongLat(other.bWebMercatorToWGS84LongLat), nErrorCount(other.nErrorCount), dfThreshold(other.dfThreshold), - m_pj(other.m_pj), m_bReversePj(other.m_bReversePj), - m_bEmitErrors(other.m_bEmitErrors), bNoTransform(other.bNoTransform), - m_eStrategy(other.m_eStrategy), + m_psLastContext(nullptr), + m_nLastContextThreadId(std::this_thread::get_id()), m_pj(other.m_pj), + m_bReversePj(other.m_bReversePj), m_bEmitErrors(other.m_bEmitErrors), + bNoTransform(other.bNoTransform), m_eStrategy(other.m_eStrategy), m_oTransformations(other.m_oTransformations), m_iCurTransformation(other.m_iCurTransformation), m_options(other.m_options) @@ -2544,7 +2549,15 @@ int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y, /* Select dynamically the best transformation for the data, if */ /* needed. */ /* -------------------------------------------------------------------- */ - auto ctx = OSRGetProjTLSContext(); + PJ_CONTEXT *ctx = m_psLastContext; + const auto nThisThreadId = std::this_thread::get_id(); + if (!ctx || nThisThreadId != m_nLastContextThreadId) + { + m_nLastContextThreadId = nThisThreadId; + m_psLastContext = OSRGetProjTLSContext(); + ctx = m_psLastContext; + } + PJ *pj = m_pj; if (!bTransformDone && !pj) { From 4c410fddd70d3388bf1ada7faa523c9d85959d15 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 20:43:21 +0200 Subject: [PATCH 05/10] ogr2ogr: implement CompositeCT::TransformWithErrorCodes() --- apps/ogr2ogr_lib.cpp | 59 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 53 insertions(+), 6 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index 42dea3bce9dc..9d28764d9d98 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -1224,20 +1224,24 @@ class GCPCoordTransformation : public OGRCoordinateTransformation class CompositeCT : public OGRCoordinateTransformation { + OGRCoordinateTransformation *const poCT1; + const bool bOwnCT1; + OGRCoordinateTransformation *const poCT2; + const bool bOwnCT2; + + // Working buffer + std::vector m_anErrorCode{}; + CompositeCT(const CompositeCT &other) : poCT1(other.poCT1 ? other.poCT1->Clone() : nullptr), bOwnCT1(true), - poCT2(other.poCT2 ? other.poCT2->Clone() : nullptr), bOwnCT2(true) + poCT2(other.poCT2 ? other.poCT2->Clone() : nullptr), bOwnCT2(true), + m_anErrorCode({}) { } CompositeCT &operator=(const CompositeCT &) = delete; public: - OGRCoordinateTransformation *poCT1; - bool bOwnCT1; - OGRCoordinateTransformation *poCT2; - bool bOwnCT2; - CompositeCT(OGRCoordinateTransformation *poCT1In, bool bOwnCT1In, OGRCoordinateTransformation *poCT2In, bool bOwnCT2In) : poCT1(poCT1In), bOwnCT1(bOwnCT1In), poCT2(poCT2In), bOwnCT2(bOwnCT2In) @@ -1299,6 +1303,35 @@ class CompositeCT : public OGRCoordinateTransformation return nResult; } + virtual int TransformWithErrorCodes(size_t nCount, double *x, double *y, + double *z, double *t, + int *panErrorCodes) override + { + if (poCT1 && poCT2 && panErrorCodes) + { + m_anErrorCode.resize(nCount); + int nResult = poCT1->TransformWithErrorCodes(nCount, x, y, z, t, + m_anErrorCode.data()); + if (nResult) + nResult = poCT2->TransformWithErrorCodes(nCount, x, y, z, t, + panErrorCodes); + for (size_t i = 0; i < nCount; ++i) + { + if (m_anErrorCode[i]) + panErrorCodes[i] = m_anErrorCode[i]; + } + return nResult; + } + int nResult = TRUE; + if (poCT1) + nResult = poCT1->TransformWithErrorCodes(nCount, x, y, z, t, + panErrorCodes); + if (nResult && poCT2) + nResult = poCT2->TransformWithErrorCodes(nCount, x, y, z, t, + panErrorCodes); + return nResult; + } + virtual OGRCoordinateTransformation *GetInverse() const override { if (!poCT1 && !poCT2) @@ -1386,6 +1419,20 @@ class AxisMappingCoordinateTransformation : public OGRCoordinateTransformation return true; } + virtual int TransformWithErrorCodes(size_t nCount, double *x, double *y, + double * /*z*/, double * /*t*/, + int *panErrorCodes) override + { + for (size_t i = 0; i < nCount; i++) + { + if (panErrorCodes) + panErrorCodes[i] = 0; + if (bSwapXY) + std::swap(x[i], y[i]); + } + return true; + } + virtual OGRCoordinateTransformation *GetInverse() const override { return new AxisMappingCoordinateTransformation(bSwapXY); From 6f88dcf12dc38211306e1e5a6982547cc740b6ad Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 22:02:25 +0200 Subject: [PATCH 06/10] Add OGRWKBTransform() for in-place coordinate transformation of WKB geometries --- autotest/cpp/test_ogr_wkb.cpp | 350 +++++++++++++++++++++++++++++++++ ogr/ogr_wkb.cpp | 351 ++++++++++++++++++++++++++++++++++ ogr/ogr_wkb.h | 42 ++++ 3 files changed, 743 insertions(+) diff --git a/autotest/cpp/test_ogr_wkb.cpp b/autotest/cpp/test_ogr_wkb.cpp index ed86e66fd6a6..dba3eedc46a0 100644 --- a/autotest/cpp/test_ogr_wkb.cpp +++ b/autotest/cpp/test_ogr_wkb.cpp @@ -524,4 +524,354 @@ INSTANTIATE_TEST_SUITE_P( OGRWKBIntersectsPessimisticFixture::ParamType> &l_info) { return std::get<6>(l_info.param); }); +class OGRWKBTransformFixture + : public test_ogr_wkb, + public ::testing::WithParamInterface< + std::tuple> +{ + public: + static std::vector< + std::tuple> + GetTupleValues() + { + return { + std::make_tuple("POINT EMPTY", wkbNDR, "POINT EMPTY", + "POINT_EMPTY_NDR"), + std::make_tuple("POINT EMPTY", wkbXDR, "POINT EMPTY", + "POINT_EMPTY_XDR"), + std::make_tuple("POINT (1 2)", wkbNDR, "POINT (2 4)", "POINT_NDR"), + std::make_tuple("POINT (1 2)", wkbXDR, "POINT (2 4)", "POINT_XDR"), + std::make_tuple("POINT Z EMPTY", wkbNDR, "POINT Z EMPTY", + "POINT_Z_EMPTY_NDR"), + std::make_tuple("POINT Z EMPTY", wkbXDR, "POINT Z EMPTY", + "POINT_Z_EMPTY_XDR"), + std::make_tuple("POINT Z (1 2 3)", wkbNDR, "POINT Z (2 4 6)", + "POINT_Z_NDR"), + std::make_tuple("POINT Z (1 2 3)", wkbXDR, "POINT Z (2 4 6)", + "POINT_Z_XDR"), + std::make_tuple("POINT M EMPTY", wkbNDR, "POINT M EMPTY", + "POINT_M_EMPTY_NDR"), + std::make_tuple("POINT M EMPTY", wkbXDR, "POINT M EMPTY", + "POINT_M_EMPTY_XDR"), + std::make_tuple("POINT M (1 2 -10)", wkbNDR, "POINT M (2 4 -10)", + "POINT_M_NDR"), + std::make_tuple("POINT M (1 2 -10)", wkbXDR, "POINT M (2 4 -10)", + "POINT_M_XDR"), + std::make_tuple("POINT ZM EMPTY", wkbNDR, "POINT ZM EMPTY", + "POINT_ZM_EMPTY_NDR"), + std::make_tuple("POINT ZM EMPTY", wkbXDR, "POINT ZM EMPTY", + "POINT_ZM_EMPTY_XDR"), + std::make_tuple("POINT ZM (1 2 3 10)", wkbNDR, + "POINT ZM (2 4 6 10)", "POINT_ZM_NDR"), + std::make_tuple("POINT ZM (1 2 3 10)", wkbXDR, + "POINT ZM (2 4 6 10)", "POINT_ZM_XDR"), + + std::make_tuple("LINESTRING EMPTY", wkbNDR, "LINESTRING EMPTY", + "LINESTRING_EMPTY"), + std::make_tuple("LINESTRING (1 2,11 12)", wkbNDR, + "LINESTRING (2 4,12 14)", "LINESTRING_NDR"), + std::make_tuple("LINESTRING (1 2,11 12)", wkbXDR, + "LINESTRING (2 4,12 14)", "LINESTRING_XDR"), + std::make_tuple("LINESTRING Z EMPTY", wkbNDR, "LINESTRING Z EMPTY", + "LINESTRING_Z_EMPTY"), + std::make_tuple("LINESTRING Z (1 2 3,11 12 13)", wkbNDR, + "LINESTRING Z (2 4 6,12 14 16)", + "LINESTRING_Z_NDR"), + std::make_tuple("LINESTRING Z (1 2 3,11 12 13)", wkbXDR, + "LINESTRING Z (2 4 6,12 14 16)", + "LINESTRING_Z_XDR"), + std::make_tuple("LINESTRING M EMPTY", wkbNDR, "LINESTRING M EMPTY", + "LINESTRING_M_EMPTY"), + std::make_tuple("LINESTRING M (1 2 -10,11 12 -20)", wkbNDR, + "LINESTRING M (2 4 -10,12 14 -20)", + "LINESTRING_M_NDR"), + std::make_tuple("LINESTRING M (1 2 -10,11 12 -20)", wkbXDR, + "LINESTRING M (2 4 -10,12 14 -20)", + "LINESTRING_M_XDR"), + std::make_tuple("LINESTRING ZM EMPTY", wkbNDR, + "LINESTRING ZM EMPTY", "LINESTRING_ZM_EMPTY"), + std::make_tuple("LINESTRING ZM (1 2 3 -10,11 12 13 -20)", wkbNDR, + "LINESTRING ZM (2 4 6 -10,12 14 16 -20)", + "LINESTRING_ZM_NDR"), + std::make_tuple("LINESTRING ZM (1 2 3 -10,11 12 13 -20)", wkbXDR, + "LINESTRING ZM (2 4 6 -10,12 14 16 -20)", + "LINESTRING_ZM_XDR"), + + // I know the polygon is invalid, but this is enough for our purposes + std::make_tuple("POLYGON EMPTY", wkbNDR, "POLYGON EMPTY", + "POLYGON_EMPTY"), + std::make_tuple("POLYGON ((1 2,11 12))", wkbNDR, + "POLYGON ((2 4,12 14))", "POLYGON_NDR"), + std::make_tuple("POLYGON ((1 2,11 12))", wkbXDR, + "POLYGON ((2 4,12 14))", "POLYGON_XDR"), + std::make_tuple("POLYGON ((1 2,11 12),(21 22,31 32))", wkbNDR, + "POLYGON ((2 4,12 14),(22 24,32 34))", + "POLYGON_TWO_RINGS"), + std::make_tuple("POLYGON Z EMPTY", wkbNDR, "POLYGON Z EMPTY", + "POLYGON_Z_EMPTY"), + std::make_tuple("POLYGON Z ((1 2 3,11 12 13))", wkbNDR, + "POLYGON Z ((2 4 6,12 14 16))", "POLYGON_Z_NDR"), + std::make_tuple("POLYGON Z ((1 2 3,11 12 13))", wkbXDR, + "POLYGON Z ((2 4 6,12 14 16))", "POLYGON_Z_XDR"), + std::make_tuple("POLYGON M EMPTY", wkbNDR, "POLYGON M EMPTY", + "POLYGON_M_EMPTY"), + std::make_tuple("POLYGON M ((1 2 -10,11 12 -20))", wkbNDR, + "POLYGON M ((2 4 -10,12 14 -20))", "POLYGON_M_NDR"), + std::make_tuple("POLYGON M ((1 2 -10,11 12 -20))", wkbXDR, + "POLYGON M ((2 4 -10,12 14 -20))", "POLYGON_M_XDR"), + std::make_tuple("POLYGON ZM EMPTY", wkbNDR, "POLYGON ZM EMPTY", + "POLYGON_ZM_EMPTY"), + std::make_tuple("POLYGON ZM ((1 2 3 -10,11 12 13 -20))", wkbNDR, + "POLYGON ZM ((2 4 6 -10,12 14 16 -20))", + "POLYGON_ZM_NDR"), + std::make_tuple("POLYGON ZM ((1 2 3 -10,11 12 13 -20))", wkbXDR, + "POLYGON ZM ((2 4 6 -10,12 14 16 -20))", + "POLYGON_ZM_XDR"), + + std::make_tuple("MULTIPOINT EMPTY", wkbNDR, "MULTIPOINT EMPTY", + "MULTIPOINT_EMPTY_NDR"), + std::make_tuple("MULTIPOINT ((1 2),(11 12))", wkbNDR, + "MULTIPOINT ((2 4),(12 14))", "MULTIPOINT_NDR"), + std::make_tuple("MULTIPOINT Z ((1 2 3),(11 12 13))", wkbXDR, + "MULTIPOINT Z ((2 4 6),(12 14 16))", + "MULTIPOINT_Z_XDR"), + + std::make_tuple("MULTILINESTRING ((1 2,11 12))", wkbNDR, + "MULTILINESTRING ((2 4,12 14))", + "MULTILINESTRING_NDR"), + + std::make_tuple("MULTIPOLYGON (((1 2,11 12)))", wkbNDR, + "MULTIPOLYGON (((2 4,12 14)))", "MULTIPOLYGON_NDR"), + + std::make_tuple("GEOMETRYCOLLECTION (POLYGON ((1 2,11 12)))", + wkbNDR, + "GEOMETRYCOLLECTION (POLYGON ((2 4,12 14)))", + "GEOMETRYCOLLECTION_NDR"), + + std::make_tuple("CIRCULARSTRING (1 2,11 12,21 22)", wkbNDR, + "CIRCULARSTRING (2 4,12 14,22 24)", + "CIRCULARSTRING_NDR"), + + std::make_tuple("COMPOUNDCURVE ((1 2,11 12))", wkbNDR, + "COMPOUNDCURVE ((2 4,12 14))", "COMPOUNDCURVE_NDR"), + + std::make_tuple("CURVEPOLYGON ((1 2,11 12,21 22,1 2))", wkbNDR, + "CURVEPOLYGON ((2 4,12 14,22 24,2 4))", + "CURVEPOLYGON_NDR"), + + std::make_tuple("MULTICURVE ((1 2,11 12))", wkbNDR, + "MULTICURVE ((2 4,12 14))", "MULTICURVE_NDR"), + + std::make_tuple("MULTISURFACE (((1 2,11 12)))", wkbNDR, + "MULTISURFACE (((2 4,12 14)))", "MULTISURFACE_NDR"), + + std::make_tuple("TRIANGLE ((1 2,11 12,21 22,1 2))", wkbNDR, + "TRIANGLE ((2 4,12 14,22 24,2 4))", "TRIANGLE_NDR"), + + std::make_tuple("POLYHEDRALSURFACE (((1 2,11 12,21 22,1 2)))", + wkbNDR, + "POLYHEDRALSURFACE (((2 4,12 14,22 24,2 4)))", + "POLYHEDRALSURFACE_NDR"), + + std::make_tuple("TIN (((1 2,11 12,21 22,1 2)))", wkbNDR, + "TIN (((2 4,12 14,22 24,2 4)))", "TIN_NDR"), + }; + } +}; + +struct MyCT final : public OGRCoordinateTransformation +{ + const bool m_bSuccess; + + explicit MyCT(bool bSuccess = true) : m_bSuccess(bSuccess) + { + } + + const OGRSpatialReference *GetSourceCS() const override + { + return nullptr; + } + + const OGRSpatialReference *GetTargetCS() const override + { + return nullptr; + } + + int Transform(size_t nCount, double *x, double *y, double *z, double *, + int *pabSuccess) override + { + for (size_t i = 0; i < nCount; ++i) + { + x[i] += 1; + y[i] += 2; + if (z) + z[i] += 3; + if (pabSuccess) + pabSuccess[i] = m_bSuccess; + } + return true; + } + + OGRCoordinateTransformation *Clone() const override + { + return new MyCT(); + } + + OGRCoordinateTransformation *GetInverse() const override + { + return nullptr; + } // unused +}; + +TEST_P(OGRWKBTransformFixture, test) +{ + const char *pszInput = std::get<0>(GetParam()); + OGRwkbByteOrder eByteOrder = std::get<1>(GetParam()); + const char *pszOutput = std::get<2>(GetParam()); + + MyCT oCT; + oCT.GetSourceCS(); // just for code coverage purpose + oCT.GetTargetCS(); // just for code coverage purpose + delete oCT.Clone(); // just for code coverage purpose + delete oCT.GetInverse(); // just for code coverage purpose + + OGRGeometry *poGeom = nullptr; + EXPECT_EQ(OGRGeometryFactory::createFromWkt(pszInput, nullptr, &poGeom), + OGRERR_NONE); + ASSERT_TRUE(poGeom != nullptr); + std::vector abyWkb(poGeom->WkbSize()); + poGeom->exportToWkb(eByteOrder, abyWkb.data(), wkbVariantIso); + delete poGeom; + + OGRWKBTransformCache oCache; + OGREnvelope3D sEnv; + EXPECT_TRUE( + OGRWKBTransform(abyWkb.data(), abyWkb.size(), &oCT, oCache, sEnv)); + const auto abyWkbOri = abyWkb; + + poGeom = nullptr; + OGRGeometryFactory::createFromWkb(abyWkb.data(), nullptr, &poGeom, + abyWkb.size()); + ASSERT_TRUE(poGeom != nullptr); + char *pszWKT = nullptr; + poGeom->exportToWkt(&pszWKT, wkbVariantIso); + delete poGeom; + EXPECT_STREQ(pszWKT, pszOutput); + CPLFree(pszWKT); + + { + CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler); + + // Truncated geometry + for (size_t i = 0; i < abyWkb.size(); ++i) + { + abyWkb = abyWkbOri; + EXPECT_FALSE(OGRWKBTransform(abyWkb.data(), i, &oCT, oCache, sEnv)); + } + + // Check altering all bytes + for (size_t i = 0; i < abyWkb.size(); ++i) + { + abyWkb = abyWkbOri; + abyWkb[i] = 0xff; + CPL_IGNORE_RET_VAL(OGRWKBTransform(abyWkb.data(), abyWkb.size(), + &oCT, oCache, sEnv)); + } + + if (abyWkb.size() > 9) + { + abyWkb = abyWkbOri; + if (!STARTS_WITH(pszInput, "POINT")) + { + // Corrupt number of sub-geometries + abyWkb[5] = 0xff; + abyWkb[6] = 0xff; + abyWkb[7] = 0xff; + abyWkb[8] = 0xff; + EXPECT_FALSE(OGRWKBTransform(abyWkb.data(), abyWkb.size(), &oCT, + oCache, sEnv)); + } + } + } +} + +INSTANTIATE_TEST_SUITE_P( + test_ogr_wkb, OGRWKBTransformFixture, + ::testing::ValuesIn(OGRWKBTransformFixture::GetTupleValues()), + [](const ::testing::TestParamInfo + &l_info) { return std::get<3>(l_info.param); }); + +TEST_F(test_ogr_wkb, OGRWKBTransformFixture_rec_collection) +{ + std::vector abyWkb; + constexpr int BEYOND_ALLOWED_RECURSION_LEVEL = 128; + for (int i = 0; i < BEYOND_ALLOWED_RECURSION_LEVEL; ++i) + { + abyWkb.push_back(wkbNDR); + abyWkb.push_back(wkbGeometryCollection); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(1); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + } + { + abyWkb.push_back(wkbNDR); + abyWkb.push_back(wkbGeometryCollection); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + abyWkb.push_back(0); + } + + MyCT oCT; + OGRWKBTransformCache oCache; + OGREnvelope3D sEnv; + EXPECT_FALSE( + OGRWKBTransform(abyWkb.data(), abyWkb.size(), &oCT, oCache, sEnv)); +} + +TEST_F(test_ogr_wkb, OGRWKBTransformFixture_ct_failure) +{ + MyCT oCT(/* bSuccess = */ false); + OGRWKBTransformCache oCache; + OGREnvelope3D sEnv; + { + OGRPoint p(1, 2); + std::vector abyWkb(p.WkbSize()); + static_cast(p).exportToWkb(wkbNDR, abyWkb.data(), + wkbVariantIso); + + EXPECT_FALSE( + OGRWKBTransform(abyWkb.data(), abyWkb.size(), &oCT, oCache, sEnv)); + } + { + OGRLineString ls; + ls.addPoint(1, 2); + std::vector abyWkb(ls.WkbSize()); + static_cast(ls).exportToWkb(wkbNDR, abyWkb.data(), + wkbVariantIso); + + EXPECT_FALSE( + OGRWKBTransform(abyWkb.data(), abyWkb.size(), &oCT, oCache, sEnv)); + } + { + OGRPolygon p; + auto poLR = std::make_unique(); + poLR->addPoint(1, 2); + p.addRing(std::move(poLR)); + std::vector abyWkb(p.WkbSize()); + static_cast(p).exportToWkb(wkbNDR, abyWkb.data(), + wkbVariantIso); + + EXPECT_FALSE( + OGRWKBTransform(abyWkb.data(), abyWkb.size(), &oCT, oCache, sEnv)); + } +} + } // namespace diff --git a/ogr/ogr_wkb.cpp b/ogr/ogr_wkb.cpp index 6b96f39e0144..c2b39d03b3ce 100644 --- a/ogr/ogr_wkb.cpp +++ b/ogr/ogr_wkb.cpp @@ -1072,6 +1072,357 @@ void OGRWKBFixupCounterClockWiseExternalRing(GByte *pabyWkb, size_t nWKBSize) pabyWkb, nWKBSize, iOffsetInOut, /* nRec = */ 0); } +/************************************************************************/ +/* OGRWKBPointUpdater() */ +/************************************************************************/ + +OGRWKBPointUpdater::OGRWKBPointUpdater() = default; + +/************************************************************************/ +/* OGRWKBIntersectsPointSequencePessimistic() */ +/************************************************************************/ + +static bool OGRWKBUpdatePointsSequence(uint8_t *data, const size_t size, + OGRWKBPointUpdater &oUpdater, + const OGRwkbByteOrder eByteOrder, + const int nDim, const bool bHasZ, + const bool bHasM, size_t &iOffsetInOut) +{ + const uint32_t nPoints = + OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut); + if (nPoints > (size - iOffsetInOut) / (nDim * sizeof(double))) + { + return false; + } + const bool bNeedSwap = OGR_SWAP(eByteOrder); + for (uint32_t j = 0; j < nPoints; j++) + { + void *pdfX = data + iOffsetInOut; + void *pdfY = data + iOffsetInOut + sizeof(double); + void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr; + void *pdfM = + bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double) + : nullptr; + if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM)) + return false; + + iOffsetInOut += nDim * sizeof(double); + } + + return true; +} + +/************************************************************************/ +/* OGRWKBVisitRingSequence() */ +/************************************************************************/ + +static bool OGRWKBVisitRingSequence(uint8_t *data, const size_t size, + OGRWKBPointUpdater &oUpdater, + const OGRwkbByteOrder eByteOrder, + const int nDim, const bool bHasZ, + const bool bHasM, size_t &iOffsetInOut) +{ + const uint32_t nRings = + OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut); + if (nRings > (size - iOffsetInOut) / sizeof(uint32_t)) + { + return false; + } + + for (uint32_t i = 0; i < nRings; ++i) + { + if (iOffsetInOut + sizeof(uint32_t) > size) + { + return false; + } + if (!OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder, nDim, + bHasZ, bHasM, iOffsetInOut)) + { + return false; + } + } + return true; +} + +/************************************************************************/ +/* OGRWKBUpdatePoints() */ +/************************************************************************/ + +static bool OGRWKBUpdatePoints(uint8_t *data, const size_t size, + OGRWKBPointUpdater &oUpdater, + size_t &iOffsetInOut, const int nRec) +{ + if (size - iOffsetInOut < MIN_WKB_SIZE) + { + return false; + } + const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(data[iOffsetInOut]); + if (!(nByteOrder == wkbXDR || nByteOrder == wkbNDR)) + { + return false; + } + const OGRwkbByteOrder eByteOrder = static_cast(nByteOrder); + + OGRwkbGeometryType eGeometryType = wkbUnknown; + OGRReadWKBGeometryType(data + iOffsetInOut, wkbVariantIso, &eGeometryType); + iOffsetInOut += 5; + const auto eFlatType = wkbFlatten(eGeometryType); + + if (eFlatType == wkbGeometryCollection || eFlatType == wkbCompoundCurve || + eFlatType == wkbCurvePolygon || eFlatType == wkbMultiPoint || + eFlatType == wkbMultiLineString || eFlatType == wkbMultiPolygon || + eFlatType == wkbMultiCurve || eFlatType == wkbMultiSurface || + eFlatType == wkbPolyhedralSurface || eFlatType == wkbTIN) + { + if (nRec == 128) + return false; + + const uint32_t nParts = + OGRWKBReadUInt32AtOffset(data, eByteOrder, iOffsetInOut); + if (nParts > (size - iOffsetInOut) / MIN_WKB_SIZE) + { + return false; + } + for (uint32_t k = 0; k < nParts; k++) + { + if (!OGRWKBUpdatePoints(data, size, oUpdater, iOffsetInOut, + nRec + 1)) + return false; + } + return true; + } + + const bool bHasZ = OGR_GT_HasZ(eGeometryType); + const bool bHasM = OGR_GT_HasM(eGeometryType); + const int nDim = 2 + (bHasZ ? 1 : 0) + (bHasM ? 1 : 0); + + if (eFlatType == wkbPoint) + { + if (size - iOffsetInOut < nDim * sizeof(double)) + return false; + void *pdfX = data + iOffsetInOut; + void *pdfY = data + iOffsetInOut + sizeof(double); + void *pdfZ = bHasZ ? data + iOffsetInOut + 2 * sizeof(double) : nullptr; + void *pdfM = + bHasM ? data + iOffsetInOut + (bHasZ ? 3 : 2) * sizeof(double) + : nullptr; + const bool bNeedSwap = OGR_SWAP(eByteOrder); + if (!oUpdater.update(bNeedSwap, pdfX, pdfY, pdfZ, pdfM)) + return false; + iOffsetInOut += nDim * sizeof(double); + return true; + } + + if (eFlatType == wkbLineString || eFlatType == wkbCircularString) + { + return OGRWKBUpdatePointsSequence(data, size, oUpdater, eByteOrder, + nDim, bHasZ, bHasM, iOffsetInOut); + } + + if (eFlatType == wkbPolygon || eFlatType == wkbTriangle) + { + return OGRWKBVisitRingSequence(data, size, oUpdater, eByteOrder, nDim, + bHasZ, bHasM, iOffsetInOut); + } + + CPLDebug("OGR", "Unknown WKB geometry type"); + return false; +} + +/** Visit all points of a WKB geometry to update them. + */ +bool OGRWKBUpdatePoints(GByte *pabyWkb, size_t nWKBSize, + OGRWKBPointUpdater &oUpdater) +{ + size_t iOffsetInOut = 0; + return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater, iOffsetInOut, + /* nRec = */ 0); +} + +/************************************************************************/ +/* OGRWKBTransformCache::clear() */ +/************************************************************************/ + +#ifdef OGR_WKB_TRANSFORM_ALL_AT_ONCE +void OGRWKBTransformCache::clear() +{ + abNeedSwap.clear(); + abIsEmpty.clear(); + apdfX.clear(); + apdfY.clear(); + apdfZ.clear(); + apdfM.clear(); + adfX.clear(); + adfY.clear(); + adfZ.clear(); + adfM.clear(); + anErrorCodes.clear(); +} +#endif + +/************************************************************************/ +/* OGRWKBTransform() */ +/************************************************************************/ + +/** Visit all points of a WKB geometry to transform them. + */ +bool OGRWKBTransform(GByte *pabyWkb, size_t nWKBSize, + OGRCoordinateTransformation *poCT, + [[maybe_unused]] OGRWKBTransformCache &oCache, + OGREnvelope3D &sEnvelope) +{ +#ifndef OGR_WKB_TRANSFORM_ALL_AT_ONCE + struct OGRWKBPointUpdaterReproj final : public OGRWKBPointUpdater + { + OGRCoordinateTransformation *m_poCT; + OGREnvelope3D &m_sEnvelope; + + explicit OGRWKBPointUpdaterReproj(OGRCoordinateTransformation *poCTIn, + OGREnvelope3D &sEnvelopeIn) + : m_poCT(poCTIn), m_sEnvelope(sEnvelopeIn) + { + } + + bool update(bool bNeedSwap, void *x, void *y, void *z, + void * /* m */) override + { + double dfX, dfY, dfZ; + memcpy(&dfX, x, sizeof(double)); + memcpy(&dfY, y, sizeof(double)); + if (bNeedSwap) + { + CPL_SWAP64PTR(&dfX); + CPL_SWAP64PTR(&dfY); + } + if (!(std::isnan(dfX) && std::isnan(dfY))) + { + if (z) + { + memcpy(&dfZ, z, sizeof(double)); + if (bNeedSwap) + { + CPL_SWAP64PTR(&dfZ); + } + } + else + dfZ = 0; + int nErrorCode = 0; + m_poCT->TransformWithErrorCodes(1, &dfX, &dfY, &dfZ, nullptr, + &nErrorCode); + if (nErrorCode) + return false; + m_sEnvelope.Merge(dfX, dfY, dfZ); + if (bNeedSwap) + { + CPL_SWAP64PTR(&dfX); + CPL_SWAP64PTR(&dfY); + CPL_SWAP64PTR(&dfZ); + } + memcpy(x, &dfX, sizeof(double)); + memcpy(y, &dfY, sizeof(double)); + if (z) + memcpy(z, &dfZ, sizeof(double)); + } + return true; + } + + private: + OGRWKBPointUpdaterReproj(const OGRWKBPointUpdaterReproj &) = delete; + OGRWKBPointUpdaterReproj & + operator=(const OGRWKBPointUpdaterReproj &) = delete; + }; + + sEnvelope = OGREnvelope3D(); + OGRWKBPointUpdaterReproj oUpdater(poCT, sEnvelope); + return OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater); + +#else + struct OGRWKBPointUpdaterReproj final : public OGRWKBPointUpdater + { + OGRWKBTransformCache &m_oCache; + + explicit OGRWKBPointUpdaterReproj(OGRWKBTransformCache &oCacheIn) + : m_oCache(oCacheIn) + { + } + + bool update(bool bNeedSwap, void *x, void *y, void *z, + void * /* m */) override + { + m_oCache.abNeedSwap.push_back(bNeedSwap); + m_oCache.apdfX.push_back(x); + m_oCache.apdfY.push_back(y); + m_oCache.apdfZ.push_back(z); + return true; + } + }; + + oCache.clear(); + OGRWKBPointUpdaterReproj oUpdater(oCache); + if (!OGRWKBUpdatePoints(pabyWkb, nWKBSize, oUpdater)) + return false; + + oCache.adfX.resize(oCache.apdfX.size()); + oCache.adfY.resize(oCache.apdfX.size()); + oCache.adfZ.resize(oCache.apdfX.size()); + + for (size_t i = 0; i < oCache.apdfX.size(); ++i) + { + memcpy(&oCache.adfX[i], oCache.apdfX[i], sizeof(double)); + memcpy(&oCache.adfY[i], oCache.apdfY[i], sizeof(double)); + if (oCache.apdfZ[i]) + memcpy(&oCache.adfZ[i], oCache.apdfZ[i], sizeof(double)); + if (oCache.abNeedSwap[i]) + { + CPL_SWAP64PTR(&oCache.adfX[i]); + CPL_SWAP64PTR(&oCache.adfY[i]); + CPL_SWAP64PTR(&oCache.adfZ[i]); + } + oCache.abIsEmpty.push_back(std::isnan(oCache.adfX[i]) && + std::isnan(oCache.adfY[i])); + } + + oCache.anErrorCodes.resize(oCache.apdfX.size()); + poCT->TransformWithErrorCodes(static_cast(oCache.apdfX.size()), + oCache.adfX.data(), oCache.adfY.data(), + oCache.adfZ.data(), nullptr, + oCache.anErrorCodes.data()); + + for (size_t i = 0; i < oCache.apdfX.size(); ++i) + { + if (!oCache.abIsEmpty[i] && oCache.anErrorCodes[i]) + return false; + } + + sEnvelope = OGREnvelope3D(); + for (size_t i = 0; i < oCache.apdfX.size(); ++i) + { + if (oCache.abIsEmpty[i]) + { + oCache.adfX[i] = std::numeric_limits::quiet_NaN(); + oCache.adfY[i] = std::numeric_limits::quiet_NaN(); + oCache.adfZ[i] = std::numeric_limits::quiet_NaN(); + } + else + { + sEnvelope.Merge(oCache.adfX[i], oCache.adfY[i], oCache.adfZ[i]); + } + if (oCache.abNeedSwap[i]) + { + CPL_SWAP64PTR(&oCache.adfX[i]); + CPL_SWAP64PTR(&oCache.adfY[i]); + CPL_SWAP64PTR(&oCache.adfZ[i]); + } + memcpy(oCache.apdfX[i], &oCache.adfX[i], sizeof(double)); + memcpy(oCache.apdfY[i], &oCache.adfY[i], sizeof(double)); + if (oCache.apdfZ[i]) + memcpy(oCache.apdfZ[i], &oCache.adfZ[i], sizeof(double)); + } + + return true; +#endif +} + /************************************************************************/ /* OGRAppendBuffer() */ /************************************************************************/ diff --git a/ogr/ogr_wkb.h b/ogr/ogr_wkb.h index 2c3a9770dc89..ef4bd409702f 100644 --- a/ogr/ogr_wkb.h +++ b/ogr/ogr_wkb.h @@ -32,6 +32,8 @@ #include "cpl_port.h" #include "ogr_core.h" +#include + bool CPL_DLL OGRWKBGetGeomType(const GByte *pabyWkb, size_t nWKBSize, bool &bNeedSwap, uint32_t &nType); bool OGRWKBPolygonGetArea(const GByte *&pabyWkb, size_t &nWKBSize, @@ -62,6 +64,46 @@ void CPL_DLL OGRWKBFixupCounterClockWiseExternalRing(GByte *pabyWkb, const GByte CPL_DLL *WKBFromEWKB(GByte *pabyEWKB, size_t nEWKBSize, size_t &nWKBSizeOut, int *pnSRIDOut); +/** Object to update point coordinates in a WKB geometry */ +class CPL_DLL OGRWKBPointUpdater +{ + public: + OGRWKBPointUpdater(); + virtual ~OGRWKBPointUpdater() = default; + + /** Update method */ + virtual bool update(bool bNeedSwap, void *x, void *y, void *z, void *m) = 0; +}; + +bool CPL_DLL OGRWKBUpdatePoints(GByte *pabyWkb, size_t nWKBSize, + OGRWKBPointUpdater &oUpdater); + +/** Transformation cache */ +struct CPL_DLL OGRWKBTransformCache +{ +#ifdef OGR_WKB_TRANSFORM_ALL_AT_ONCE + std::vector abNeedSwap{}; + std::vector abIsEmpty{}; + std::vector apdfX{}; + std::vector apdfY{}; + std::vector apdfZ{}; + std::vector apdfM{}; + std::vector adfX{}; + std::vector adfY{}; + std::vector adfZ{}; + std::vector adfM{}; + std::vector anErrorCodes{}; + + void clear(); +#endif +}; + +class OGRCoordinateTransformation; +bool CPL_DLL OGRWKBTransform(GByte *pabyWkb, size_t nWKBSize, + OGRCoordinateTransformation *poCT, + OGRWKBTransformCache &oCache, + OGREnvelope3D &sEnvelope); + /************************************************************************/ /* OGRAppendBuffer */ /************************************************************************/ From d4cafcd7d8e63ccfadad2a77c44d66c3e501a7d8 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 22:02:46 +0200 Subject: [PATCH 07/10] GPKG: ST_Transform(): use OGRWKBTransform() --- ogr/ogrsf_frmts/gpkg/ogr_geopackage.h | 3 + .../gpkg/ogrgeopackagedatasource.cpp | 35 ++++++-- ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.cpp | 81 +++++++++++++++++++ ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.h | 4 + 4 files changed, 115 insertions(+), 8 deletions(-) diff --git a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h index 46fc95e25509..720bdb56e9ed 100644 --- a/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h +++ b/ogr/ogrsf_frmts/gpkg/ogr_geopackage.h @@ -37,6 +37,7 @@ #include "cpl_threadsafe_queue.hpp" #include "ograrrowarrayhelper.h" #include "ogr_p.h" +#include "ogr_wkb.h" #include #include @@ -185,6 +186,8 @@ class GDALGeoPackageDataset final : public OGRSQLiteBaseDataSource, int m_nLastCachedCTSrcSRId = -1; int m_nLastCachedCTDstSRId = -1; std::unique_ptr m_poLastCachedCT{}; + OGRWKBTransformCache m_oWKBTransformCache{}; + std::vector m_abyWKBTransformCache{}; int m_nOverviewCount = 0; GDALGeoPackageDataset **m_papoOverviewDS = nullptr; diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp index 61ee94970372..3bc3888fbb9f 100644 --- a/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp +++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackagedatasource.cpp @@ -8905,21 +8905,40 @@ void OGRGeoPackageTransform(sqlite3_context *pContext, int argc, poCT = poDS->m_poLastCachedCT.get(); } - auto poGeom = std::unique_ptr( - GPkgGeometryToOGR(pabyBLOB, nBLOBLen, nullptr)); - if (poGeom == nullptr) + if (sHeader.nHeaderLen >= 8) { - // Try also spatialite geometry blobs - OGRGeometry *poGeomSpatialite = nullptr; - if (OGRSQLiteImportSpatiaLiteGeometry(pabyBLOB, nBLOBLen, - &poGeomSpatialite) != OGRERR_NONE) + std::vector &abyNewBLOB = poDS->m_abyWKBTransformCache; + abyNewBLOB.resize(nBLOBLen); + memcpy(abyNewBLOB.data(), pabyBLOB, nBLOBLen); + + OGREnvelope3D oEnv3d; + if (!OGRWKBTransform(abyNewBLOB.data() + sHeader.nHeaderLen, + nBLOBLen - sHeader.nHeaderLen, poCT, + poDS->m_oWKBTransformCache, oEnv3d) || + !GPkgUpdateHeader(abyNewBLOB.data(), nBLOBLen, nDestSRID, + oEnv3d.MinX, oEnv3d.MaxX, oEnv3d.MinY, + oEnv3d.MaxY, oEnv3d.MinZ, oEnv3d.MaxZ)) { CPLError(CE_Failure, CPLE_AppDefined, "Invalid geometry"); sqlite3_result_blob(pContext, nullptr, 0, nullptr); return; } - poGeom.reset(poGeomSpatialite); + + sqlite3_result_blob(pContext, abyNewBLOB.data(), nBLOBLen, + SQLITE_TRANSIENT); + return; + } + + // Try also spatialite geometry blobs + OGRGeometry *poGeomSpatialite = nullptr; + if (OGRSQLiteImportSpatiaLiteGeometry(pabyBLOB, nBLOBLen, + &poGeomSpatialite) != OGRERR_NONE) + { + CPLError(CE_Failure, CPLE_AppDefined, "Invalid geometry"); + sqlite3_result_blob(pContext, nullptr, 0, nullptr); + return; } + auto poGeom = std::unique_ptr(poGeomSpatialite); if (poGeom->transform(poCT) != OGRERR_NONE) { diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.cpp b/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.cpp index 70210463a029..4dd2b74f13a9 100644 --- a/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.cpp +++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.cpp @@ -487,6 +487,87 @@ OGRErr GPkgHeaderFromWKB(const GByte *pabyGpkg, size_t nGpkgLen, return OGRERR_NONE; } +bool GPkgUpdateHeader(GByte *pabyGpkg, size_t nGpkgLen, int nSrsId, double MinX, + double MaxX, double MinY, double MaxY, double MinZ, + double MaxZ) +{ + CPLAssert(nGpkgLen >= 8); + + /* Flags */ + const GByte byFlags = pabyGpkg[3]; + const auto eByteOrder = static_cast(byFlags & 0x01); + const OGRBoolean bSwap = OGR_SWAP(eByteOrder); + + /* SrsId */ + if (bSwap) + { + nSrsId = CPL_SWAP32(nSrsId); + } + memcpy(pabyGpkg + 4, &nSrsId, 4); + + /* Envelope */ + const int iEnvelope = (byFlags & (0x07 << 1)) >> 1; + int nEnvelopeDim = 0; + if (iEnvelope) + { + if (iEnvelope == 1) + { + nEnvelopeDim = 2; /* 2D envelope */ + } + else if (iEnvelope == 2) + { + nEnvelopeDim = 3; /* 2D+Z envelope */ + } + else if (iEnvelope == 3) + { + nEnvelopeDim = 3; /* 2D+M envelope */ + } + else if (iEnvelope == 4) + { + nEnvelopeDim = 4; /* 2D+ZM envelope */ + } + else + { + return false; + } + } + else + { + return true; + } + + if (nGpkgLen < static_cast(8 + 8 * 2 * nEnvelopeDim)) + { + // Not enough bytes + return false; + } + + /* Envelope */ + if (bSwap) + { + CPL_SWAPDOUBLE(&(MinX)); + CPL_SWAPDOUBLE(&(MaxX)); + CPL_SWAPDOUBLE(&(MinY)); + CPL_SWAPDOUBLE(&(MaxY)); + CPL_SWAPDOUBLE(&(MinZ)); + CPL_SWAPDOUBLE(&(MaxZ)); + } + + double *padPtr = reinterpret_cast(pabyGpkg + 8); + memcpy(&padPtr[0], &MinX, sizeof(double)); + memcpy(&padPtr[1], &MaxX, sizeof(double)); + memcpy(&padPtr[2], &MinY, sizeof(double)); + memcpy(&padPtr[3], &MaxY, sizeof(double)); + + if (iEnvelope == 2 || iEnvelope == 4) + { + memcpy(&padPtr[4], &MinZ, sizeof(double)); + memcpy(&padPtr[5], &MaxZ, sizeof(double)); + } + + return true; +} + OGRGeometry *GPkgGeometryToOGR(const GByte *pabyGpkg, size_t nGpkgLen, OGRSpatialReference *poSrs) { diff --git a/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.h b/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.h index 19b07248994e..57c252c0a3b5 100644 --- a/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.h +++ b/ogr/ogrsf_frmts/gpkg/ogrgeopackageutility.h @@ -72,4 +72,8 @@ bool OGRGeoPackageGetHeader(sqlite3_context *pContext, int /*argc*/, bool bNeedExtent, bool bNeedExtent3D, int iGeomIdx = 0); +bool GPkgUpdateHeader(GByte *pabyGpkg, size_t nGpkgLen, int nSrsId, double MinX, + double MaxX, double MinY, double MaxY, double MinZ, + double MaxZ); + #endif From 620cbba885e5225a0ab1811615f3d599b1e8150c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 22:03:12 +0200 Subject: [PATCH 08/10] Add OGRGeometryFactory::isTransformWithOptionsRegularTransform() --- ogr/ogr_geometry.h | 6 +++ ogr/ogrgeometryfactory.cpp | 89 ++++++++++++++++++++++++-------------- 2 files changed, 62 insertions(+), 33 deletions(-) diff --git a/ogr/ogr_geometry.h b/ogr/ogr_geometry.h index e1f9de8c3575..0c680130791d 100644 --- a/ogr/ogr_geometry.h +++ b/ogr/ogr_geometry.h @@ -4285,6 +4285,12 @@ class CPL_DLL OGRGeometryFactory ~TransformWithOptionsCache(); }; + //! @cond Doxygen_Suppress + static bool isTransformWithOptionsRegularTransform( + const OGRSpatialReference *poSourceCRS, + const OGRSpatialReference *poTargetCRS, CSLConstList papszOptions); + //! @endcond + static OGRGeometry *transformWithOptions( const OGRGeometry *poSrcGeom, OGRCoordinateTransformation *poCT, char **papszOptions, diff --git a/ogr/ogrgeometryfactory.cpp b/ogr/ogrgeometryfactory.cpp index 532028a25507..ff9bad541cb8 100644 --- a/ogr/ogrgeometryfactory.cpp +++ b/ogr/ogrgeometryfactory.cpp @@ -3855,6 +3855,50 @@ OGRGeometryFactory::TransformWithOptionsCache::~TransformWithOptionsCache() { } +/************************************************************************/ +/* isTransformWithOptionsRegularTransform() */ +/************************************************************************/ + +//! @cond Doxygen_Suppress +/*static */ +bool OGRGeometryFactory::isTransformWithOptionsRegularTransform( + [[maybe_unused]] const OGRSpatialReference *poSourceCRS, + [[maybe_unused]] const OGRSpatialReference *poTargetCRS, + CSLConstList papszOptions) +{ + if (papszOptions) + return false; + +#ifdef HAVE_GEOS + if (poSourceCRS->IsProjected() && poTargetCRS->IsGeographic() && + poTargetCRS->GetAxisMappingStrategy() == OAMS_TRADITIONAL_GIS_ORDER && + // check that angular units is degree + std::fabs(poTargetCRS->GetAngularUnits(nullptr) - + CPLAtof(SRS_UA_DEGREE_CONV)) <= + 1e-8 * CPLAtof(SRS_UA_DEGREE_CONV)) + { + double dfWestLong = 0.0; + double dfSouthLat = 0.0; + double dfEastLong = 0.0; + double dfNorthLat = 0.0; + if (poSourceCRS->GetAreaOfUse(&dfWestLong, &dfSouthLat, &dfEastLong, + &dfNorthLat, nullptr) && + !(dfSouthLat == -90.0 || dfNorthLat == 90.0 || + dfWestLong == -180.0 || dfEastLong == 180.0 || + dfWestLong > dfEastLong)) + { + // Not a global geographic CRS + return true; + } + return false; + } +#endif + + return true; +} + +//! @endcond + /************************************************************************/ /* transformWithOptions() */ /************************************************************************/ @@ -3889,41 +3933,20 @@ OGRGeometry *OGRGeometryFactory::transformWithOptions( cache.d->poSourceCRS = poSourceCRS; cache.d->poTargetCRS = poTargetCRS; cache.d->poCT = poCT; - if (poSourceCRS && poTargetCRS && poSourceCRS->IsProjected() && - poTargetCRS->IsGeographic() && - poTargetCRS->GetAxisMappingStrategy() == - OAMS_TRADITIONAL_GIS_ORDER && - // check that angular units is degree - std::fabs(poTargetCRS->GetAngularUnits(nullptr) - - CPLAtof(SRS_UA_DEGREE_CONV)) <= - 1e-8 * CPLAtof(SRS_UA_DEGREE_CONV)) + if (poSourceCRS && poTargetCRS && + !isTransformWithOptionsRegularTransform( + poSourceCRS, poTargetCRS, papszOptions)) { - double dfWestLong = 0.0; - double dfSouthLat = 0.0; - double dfEastLong = 0.0; - double dfNorthLat = 0.0; - if (poSourceCRS->GetAreaOfUse(&dfWestLong, &dfSouthLat, - &dfEastLong, &dfNorthLat, - nullptr) && - !(dfSouthLat == -90.0 || dfNorthLat == 90.0 || - dfWestLong == -180.0 || dfEastLong == 180.0 || - dfWestLong > dfEastLong)) + cache.d->poRevCT.reset(OGRCreateCoordinateTransformation( + poTargetCRS, poSourceCRS)); + cache.d->bIsNorthPolar = false; + cache.d->bIsPolar = false; + cache.d->poRevCT.reset(poCT->GetInverse()); + if (cache.d->poRevCT && + IsPolarToGeographic(poCT, cache.d->poRevCT.get(), + cache.d->bIsNorthPolar)) { - // Not a global geographic CRS - } - else - { - cache.d->poRevCT.reset(OGRCreateCoordinateTransformation( - poTargetCRS, poSourceCRS)); - cache.d->bIsNorthPolar = false; - cache.d->bIsPolar = false; - cache.d->poRevCT.reset(poCT->GetInverse()); - if (cache.d->poRevCT && - IsPolarToGeographic(poCT, cache.d->poRevCT.get(), - cache.d->bIsNorthPolar)) - { - cache.d->bIsPolar = true; - } + cache.d->bIsPolar = true; } } } From d4a125f92477686dad69b32fc872115cafad27c6 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 23:45:25 +0200 Subject: [PATCH 09/10] ogr2ogr: use OGRWKBTransform(), when possible, to speed-up -t_srs in Arrow code path Also use multi-threaded coordinate transformation by splitting the features within each batch in several sub-batches, each processed in its own thread. --- apps/ogr2ogr_lib.cpp | 329 +++++++++++++++++++++---- autotest/utilities/test_ogr2ogr_lib.py | 102 ++++++++ 2 files changed, 381 insertions(+), 50 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index 9d28764d9d98..5aafc7b22e3b 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -40,6 +40,7 @@ #include #include +#include #include #include #include @@ -47,6 +48,7 @@ #include #include #include +#include #include #include @@ -69,8 +71,10 @@ #include "ogr_p.h" #include "ogr_recordbatch.h" #include "ogr_spatialref.h" +#include "ogrlayerarrow.h" #include "ogrlayerdecorator.h" #include "ogrsf_frmts.h" +#include "ogr_wkb.h" typedef enum { @@ -520,38 +524,39 @@ class SetupTargetLayer bool &bError); public: - GDALDataset *m_poSrcDS; - GDALDataset *m_poDstDS; - char **m_papszLCO; - OGRSpatialReference *m_poOutputSRS; + GDALDataset *m_poSrcDS = nullptr; + GDALDataset *m_poDstDS = nullptr; + CSLConstList m_papszLCO = nullptr; + const OGRSpatialReference *m_poUserSourceSRS = nullptr; + const OGRSpatialReference *m_poOutputSRS = nullptr; bool m_bTransform = false; - bool m_bNullifyOutputSRS; + bool m_bNullifyOutputSRS = false; bool m_bSelFieldsSet = false; - char **m_papszSelFields; - bool m_bAppend; - bool m_bAddMissingFields; - int m_eGType; - GeomTypeConversion m_eGeomTypeConversion; - int m_nCoordDim; - bool m_bOverwrite; - char **m_papszFieldTypesToString; - char **m_papszMapFieldType; - bool m_bUnsetFieldWidth; - bool m_bExplodeCollections; - const char *m_pszZField; - char **m_papszFieldMap; - const char *m_pszWHERE; - bool m_bExactFieldNameMatch; - bool m_bQuiet; - bool m_bForceNullable; - bool m_bResolveDomains; - bool m_bUnsetDefault; - bool m_bUnsetFid; - bool m_bPreserveFID; - bool m_bCopyMD; - bool m_bNativeData; - bool m_bNewDataSource; - const char *m_pszCTPipeline; + CSLConstList m_papszSelFields = nullptr; + bool m_bAppend = false; + bool m_bAddMissingFields = false; + int m_eGType = 0; + GeomTypeConversion m_eGeomTypeConversion = GTC_DEFAULT; + int m_nCoordDim = 0; + bool m_bOverwrite = false; + CSLConstList m_papszFieldTypesToString = nullptr; + CSLConstList m_papszMapFieldType = nullptr; + bool m_bUnsetFieldWidth = false; + bool m_bExplodeCollections = false; + const char *m_pszZField = nullptr; + CSLConstList m_papszFieldMap = nullptr; + const char *m_pszWHERE = nullptr; + bool m_bExactFieldNameMatch = false; + bool m_bQuiet = false; + bool m_bForceNullable = false; + bool m_bResolveDomains = false; + bool m_bUnsetDefault = false; + bool m_bUnsetFid = false; + bool m_bPreserveFID = false; + bool m_bCopyMD = false; + bool m_bNativeData = false; + bool m_bNewDataSource = false; + const char *m_pszCTPipeline = nullptr; std::unique_ptr Setup(OGRLayer *poSrcLayer, const char *pszNewLayerName, @@ -560,11 +565,10 @@ class SetupTargetLayer class LayerTranslator { - static bool TranslateArrow(const TargetLayerInfo *psInfo, - GIntBig nCountLayerFeatures, - GIntBig *pnReadFeatureCount, - GDALProgressFunc pfnProgress, void *pProgressArg, - const GDALVectorTranslateOptions *psOptions); + bool TranslateArrow(TargetLayerInfo *psInfo, GIntBig nCountLayerFeatures, + GIntBig *pnReadFeatureCount, + GDALProgressFunc pfnProgress, void *pProgressArg, + const GDALVectorTranslateOptions *psOptions); public: GDALDataset *m_poSrcDS = nullptr; @@ -572,9 +576,9 @@ class LayerTranslator bool m_bTransform = false; bool m_bWrapDateline = false; CPLString m_osDateLineOffset{}; - OGRSpatialReference *m_poOutputSRS = nullptr; + const OGRSpatialReference *m_poOutputSRS = nullptr; bool m_bNullifyOutputSRS = false; - OGRSpatialReference *m_poUserSourceSRS = nullptr; + const OGRSpatialReference *m_poUserSourceSRS = nullptr; OGRCoordinateTransformation *m_poGCPCoordTrans = nullptr; int m_eGType = -1; GeomTypeConversion m_eGeomTypeConversion = GTC_DEFAULT; @@ -586,20 +590,20 @@ class LayerTranslator OGRGeometry *m_poClipSrcOri = nullptr; bool m_bWarnedClipSrcSRS = false; - std::unique_ptr m_poClipSrcReprojectedToSrcSRS; + std::unique_ptr m_poClipSrcReprojectedToSrcSRS{}; const OGRSpatialReference *m_poClipSrcReprojectedToSrcSRS_SRS = nullptr; OGREnvelope m_oClipSrcEnv{}; OGRGeometry *m_poClipDstOri = nullptr; bool m_bWarnedClipDstSRS = false; - std::unique_ptr m_poClipDstReprojectedToDstSRS; + std::unique_ptr m_poClipDstReprojectedToDstSRS{}; const OGRSpatialReference *m_poClipDstReprojectedToDstSRS_SRS = nullptr; OGREnvelope m_oClipDstEnv{}; bool m_bExplodeCollections = false; bool m_bNativeData = false; GIntBig m_nLimit = -1; - OGRGeometryFactory::TransformWithOptionsCache m_transformWithOptionsCache; + OGRGeometryFactory::TransformWithOptionsCache m_transformWithOptionsCache{}; bool Translate(OGRFeature *poFeatureIn, TargetLayerInfo *psInfo, GIntBig nCountLayerFeatures, GIntBig *pnReadFeatureCount, @@ -2910,6 +2914,7 @@ GDALDatasetH GDALVectorTranslate(const char *pszDest, GDALDatasetH hDstDS, oSetup.m_poOutputSRS = oOutputSRSHolder.get(); oSetup.m_bTransform = psOptions->bTransform; oSetup.m_bNullifyOutputSRS = psOptions->bNullifyOutputSRS; + oSetup.m_poUserSourceSRS = poSourceSRS; oSetup.m_bSelFieldsSet = psOptions->bSelFieldsSet; oSetup.m_papszSelFields = psOptions->aosSelFields.List(); oSetup.m_bAppend = bAppend; @@ -3807,8 +3812,8 @@ static OGRwkbGeometryType ConvertType(GeomTypeConversion eGeomTypeConversion, static void DoFieldTypeConversion(GDALDataset *poDstDS, OGRFieldDefn &oFieldDefn, - char **papszFieldTypesToString, - char **papszMapFieldType, + CSLConstList papszFieldTypesToString, + CSLConstList papszMapFieldType, bool bUnsetFieldWidth, bool bQuiet, bool bForceNullable, bool bUnsetDefault) { @@ -3953,6 +3958,47 @@ static void DoFieldTypeConversion(GDALDataset *poDstDS, } } +/************************************************************************/ +/* GetArrowGeomFieldIndex() */ +/************************************************************************/ + +static int GetArrowGeomFieldIndex(const struct ArrowSchema *psLayerSchema, + const char *pszFieldName) +{ + if (strcmp(psLayerSchema->format, "+s") == 0) // struct + { + for (int i = 0; i < psLayerSchema->n_children; ++i) + { + const auto psSchema = psLayerSchema->children[i]; + if (strcmp(psSchema->format, "z") == 0) // binary + { + if (strcmp(psSchema->name, pszFieldName) == 0) + { + return i; + } + else + { + // Check if ARROW:extension:name = ogc.wkb or geoarrow.wkb + const char *pabyMetadata = psSchema->metadata; + if (pabyMetadata) + { + const auto oMetadata = + OGRParseArrowMetadata(pabyMetadata); + auto oIter = oMetadata.find(ARROW_EXTENSION_NAME_KEY); + if (oIter != oMetadata.end() && + (oIter->second == EXTENSION_NAME_OGC_WKB || + oIter->second == EXTENSION_NAME_GEOARROW_WKB)) + { + return i; + } + } + } + } + } + } + return -1; +} + /************************************************************************/ /* SetupTargetLayer::CanUseWriteArrowBatch() */ /************************************************************************/ @@ -3980,11 +4026,11 @@ bool SetupTargetLayer::CanUseWriteArrowBatch( !psOptions->aosLCO.FetchNameValue("BATCH_SIZE") && CPLTestBool(CPLGetConfigOption("OGR2OGR_USE_ARROW_API", "YES"))) || CPLTestBool(CPLGetConfigOption("OGR2OGR_USE_ARROW_API", "NO"))) && - !psOptions->bSkipFailures && !psOptions->bTransform && - !psOptions->poClipSrc && !psOptions->poClipDst && - psOptions->oGCPs.nGCPCount == 0 && !psOptions->bWrapDateline && - !m_papszSelFields && !m_bAddMissingFields && - m_eGType == GEOMTYPE_UNCHANGED && psOptions->eGeomOp == GEOMOP_NONE && + !psOptions->bSkipFailures && !psOptions->poClipSrc && + !psOptions->poClipDst && psOptions->oGCPs.nGCPCount == 0 && + !psOptions->bWrapDateline && !m_papszSelFields && + !m_bAddMissingFields && m_eGType == GEOMTYPE_UNCHANGED && + psOptions->eGeomOp == GEOMOP_NONE && m_eGeomTypeConversion == GTC_DEFAULT && m_nCoordDim < 0 && !m_papszFieldTypesToString && !m_papszMapFieldType && !m_bUnsetFieldWidth && !m_bExplodeCollections && !m_pszZField && @@ -3993,6 +4039,26 @@ bool SetupTargetLayer::CanUseWriteArrowBatch( psOptions->dfXYRes == OGRGeomCoordinatePrecision::UNKNOWN && !psOptions->bMakeValid && !psOptions->bSkipInvalidGeom) { + if (psOptions->bTransform) + { + // To simplify implementation for now + if (poSrcLayer->GetLayerDefn()->GetGeomFieldCount() != 1 || + poDstLayer->GetLayerDefn()->GetGeomFieldCount() != 1) + { + return false; + } + auto poSrcSRS = m_poUserSourceSRS ? m_poUserSourceSRS + : poSrcLayer->GetLayerDefn() + ->GetGeomFieldDefn(0) + ->GetSpatialRef(); + if (!poSrcSRS || + !OGRGeometryFactory::isTransformWithOptionsRegularTransform( + poSrcSRS, m_poOutputSRS, nullptr)) + { + return false; + } + } + struct ArrowArrayStream streamSrc; const char *const apszOptions[] = {"SILENCE_GET_SCHEMA_ERROR=YES", nullptr}; @@ -4001,6 +4067,15 @@ bool SetupTargetLayer::CanUseWriteArrowBatch( struct ArrowSchema schemaSrc; if (streamSrc.get_schema(&streamSrc, &schemaSrc) == 0) { + if (psOptions->bTransform && + GetArrowGeomFieldIndex(&schemaSrc, + poSrcLayer->GetGeometryColumn()) < 0) + { + schemaSrc.release(&schemaSrc); + streamSrc.release(&streamSrc); + return false; + } + std::string osErrorMsg; if (poDstLayer->IsArrowSchemaSupported(&schemaSrc, nullptr, osErrorMsg)) @@ -5699,7 +5774,7 @@ SetupCT(TargetLayerInfo *psInfo, OGRLayer *poSrcLayer, bool bTransform, /************************************************************************/ bool LayerTranslator::TranslateArrow( - const TargetLayerInfo *psInfo, GIntBig nCountLayerFeatures, + TargetLayerInfo *psInfo, GIntBig nCountLayerFeatures, GIntBig *pnReadFeatureCount, GDALProgressFunc pfnProgress, void *pProgressArg, const GDALVectorTranslateOptions *psOptions) { @@ -5749,10 +5824,52 @@ bool LayerTranslator::TranslateArrow( return false; } + int iArrowGeomFieldIndex = -1; + if (m_bTransform) + { + iArrowGeomFieldIndex = GetArrowGeomFieldIndex( + &schema, psInfo->m_poSrcLayer->GetGeometryColumn()); + if (!SetupCT(psInfo, psInfo->m_poSrcLayer, m_bTransform, + m_bWrapDateline, m_osDateLineOffset, m_poUserSourceSRS, + nullptr, m_poOutputSRS, m_poGCPCoordTrans, false)) + { + return false; + } + } + bool bRet = true; GIntBig nCount = 0; bool bGoOn = true; + std::vector abyModifiedWKB; + const int nNumReprojectionThreads = []() + { + const int nNumCPUs = CPLGetNumCPUs(); + if (nNumCPUs <= 1) + { + return 1; + } + else + { + const char *pszNumThreads = + CPLGetConfigOption("GDAL_NUM_THREADS", nullptr); + if (pszNumThreads) + { + if (EQUAL(pszNumThreads, "ALL_CPUS")) + return CPLGetNumCPUs(); + return std::min(atoi(pszNumThreads), 1024); + } + else + { + return std::max(2, nNumCPUs / 2); + } + } + }(); + + // Somewhat arbitrary threshold (config option only/mostly for autotest purposes) + const int MIN_FEATURES_FOR_THREADED_REPROJ = atoi(CPLGetConfigOption( + "OGR2OGR_MIN_FEATURES_FOR_THREADED_REPROJ", "10000")); + while (bGoOn) { struct ArrowArray array; @@ -5789,9 +5906,121 @@ bool LayerTranslator::TranslateArrow( nCount += array.length; } + const auto nArrayLength = array.length; + + // Coordinate reprojection + const void *backupGeomArrayBuffers2 = nullptr; + if (m_bTransform) + { + auto *psGeomArray = array.children[iArrowGeomFieldIndex]; + GByte *pabyWKB = static_cast( + const_cast(psGeomArray->buffers[2])); + const uint32_t *panOffsets = + static_cast(psGeomArray->buffers[1]); + auto poCT = psInfo->m_aoReprojectionInfo[0].m_poCT.get(); + + try + { + abyModifiedWKB.resize(panOffsets[nArrayLength]); + } + catch (const std::exception &) + { + CPLError(CE_Failure, CPLE_OutOfMemory, "Out of memory"); + bRet = false; + if (array.release) + array.release(&array); + break; + } + memcpy(abyModifiedWKB.data(), pabyWKB, panOffsets[nArrayLength]); + backupGeomArrayBuffers2 = psGeomArray->buffers[2]; + psGeomArray->buffers[2] = abyModifiedWKB.data(); + + std::atomic atomicRet{true}; + const auto oReprojectionLambda = + [psGeomArray, nArrayLength, panOffsets, &atomicRet, + &abyModifiedWKB, &poCT](int iThread, int nThreads) + { + OGRWKBTransformCache oCache; + OGREnvelope3D sEnv3D; + auto poThisCT = + std::unique_ptr(poCT->Clone()); + if (!poThisCT) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Cannot clone OGRCoordinateTransformation"); + atomicRet = false; + return; + } + + const GByte *pabyValidity = + static_cast(psGeomArray->buffers[0]); + const size_t iStart = + static_cast(iThread * nArrayLength / nThreads); + const size_t iMax = static_cast( + (iThread + 1) * nArrayLength / nThreads); + for (size_t i = iStart; i < iMax; ++i) + { + const size_t iShifted = + static_cast(i + psGeomArray->offset); + if (!pabyValidity || (pabyValidity[iShifted >> 8] & + (1 << (iShifted % 8))) != 0) + { + const auto nWKBSize = + panOffsets[iShifted + 1] - panOffsets[iShifted]; + if (!OGRWKBTransform( + abyModifiedWKB.data() + panOffsets[iShifted], + nWKBSize, poThisCT.get(), oCache, sEnv3D)) + { + CPLError(CE_Failure, CPLE_AppDefined, + "Reprojection failed"); + atomicRet = false; + break; + } + } + } + }; + + if (nArrayLength >= MIN_FEATURES_FOR_THREADED_REPROJ && + nNumReprojectionThreads >= 2) + { + std::vector oThreads; + for (int iThread = 0; iThread < nNumReprojectionThreads; + ++iThread) + { + oThreads.emplace_back(oReprojectionLambda, iThread, + nNumReprojectionThreads); + } + for (auto &oThread : oThreads) + { + oThread.join(); + } + } + else + { + oReprojectionLambda(0, 1); + } + + bRet = atomicRet; + if (!bRet) + { + psGeomArray->buffers[2] = backupGeomArrayBuffers2; + if (array.release) + array.release(&array); + break; + } + } + // Write batch to target layer - if (!psInfo->m_poDstLayer->WriteArrowBatch( - &schema, &array, aosOptionsWriteArrowBatch.List())) + const bool bWriteOK = psInfo->m_poDstLayer->WriteArrowBatch( + &schema, &array, aosOptionsWriteArrowBatch.List()); + + if (backupGeomArrayBuffers2) + { + auto *psGeomArray = array.children[iArrowGeomFieldIndex]; + psGeomArray->buffers[2] = backupGeomArrayBuffers2; + } + + if (!bWriteOK) { CPLError(CE_Failure, CPLE_AppDefined, "WriteArrowBatch() failed"); if (array.release) diff --git a/autotest/utilities/test_ogr2ogr_lib.py b/autotest/utilities/test_ogr2ogr_lib.py index dbfec4cba68f..f51f7cd1ace3 100755 --- a/autotest/utilities/test_ogr2ogr_lib.py +++ b/autotest/utilities/test_ogr2ogr_lib.py @@ -2829,3 +2829,105 @@ def test_ogr2ogr_lib_skip_invalid(tmp_vsimem): lyr = ds.GetLayer(0) assert lyr.GetFeatureCount() == 1 ds = None + + +############################################################################### +# Test -t_srs in Arrow code path + + +@gdaltest.enable_exceptions() +@pytest.mark.parametrize("force_reproj_threading", [False, True]) +@pytest.mark.parametrize("source_driver", ["GPKG", "Parquet"]) +def test_ogr2ogr_lib_reproject_arrow(tmp_vsimem, source_driver, force_reproj_threading): + + src_driver = gdal.GetDriverByName(source_driver) + if src_driver is None: + pytest.skip(f"{source_driver} is not available") + src_filename = str(tmp_vsimem / ("in." + source_driver.lower())) + with src_driver.Create(src_filename, 0, 0, 0, gdal.GDT_Unknown) as srcDS: + srs = osr.SpatialReference() + srs.ImportFromEPSG(32631) + srcLayer = srcDS.CreateLayer("test", srs=srs) + f = ogr.Feature(srcLayer.GetLayerDefn()) + f.SetGeometry(ogr.CreateGeometryFromWkt("POINT(500000 4500000)")) + srcLayer.CreateFeature(f) + f = ogr.Feature(srcLayer.GetLayerDefn()) + srcLayer.CreateFeature(f) + f = ogr.Feature(srcLayer.GetLayerDefn()) + f.SetGeometry(ogr.CreateGeometryFromWkt("POINT(500000 4000000)")) + srcLayer.CreateFeature(f) + + config_options = {"CPL_DEBUG": "ON", "OGR2OGR_USE_ARROW_API": "YES"} + if force_reproj_threading: + config_options["OGR2OGR_MIN_FEATURES_FOR_THREADED_REPROJ"] = "0" + + with gdal.OpenEx(src_filename) as src_ds: + for i in range(2): + + got_msg = [] + + def my_handler(errorClass, errno, msg): + got_msg.append(msg) + return + + with gdaltest.error_handler(my_handler), gdaltest.config_options( + config_options + ): + ds = gdal.VectorTranslate( + "", src_ds, format="Memory", dstSRS="EPSG:4326" + ) + + assert "OGR2OGR: Using WriteArrowBatch()" in got_msg + + lyr = ds.GetLayer(0) + assert lyr.GetFeatureCount() == 3 + f = lyr.GetNextFeature() + ogrtest.check_feature_geometry(f, "POINT(3 40.65085651557158)") + f = lyr.GetNextFeature() + assert f.GetGeometryRef() is None + f = lyr.GetNextFeature() + ogrtest.check_feature_geometry(f, "POINT(3 36.14471809881776)") + + +############################################################################### +# Test -t_srs in Arrow code path in a situation where it cannot be triggered +# currently (source CRS is crossing anti-meridian) + + +@gdaltest.enable_exceptions() +@pytest.mark.require_geos +@pytest.mark.require_driver("GPKG") +def test_ogr2ogr_lib_reproject_arrow_optim_cannot_trigger(tmp_vsimem): + + src_filename = str(tmp_vsimem / "in.gpkg") + with gdal.GetDriverByName("GPKG").Create( + src_filename, 0, 0, 0, gdal.GDT_Unknown + ) as srcDS: + srs = osr.SpatialReference() + srs.ImportFromEPSG(32660) + srcLayer = srcDS.CreateLayer("test", srs=srs) + f = ogr.Feature(srcLayer.GetLayerDefn()) + f.SetGeometry( + ogr.CreateGeometryFromWkt( + "LINESTRING(657630.64 4984896.17,815261.43 4990738.26)" + ) + ) + srcLayer.CreateFeature(f) + + got_msg = [] + + def my_handler(errorClass, errno, msg): + got_msg.append(msg) + return + + config_options = {"CPL_DEBUG": "ON", "OGR2OGR_USE_ARROW_API": "YES"} + with gdaltest.error_handler(my_handler), gdaltest.config_options(config_options): + ds = gdal.VectorTranslate("", src_filename, format="Memory", dstSRS="EPSG:4326") + + assert "OGR2OGR: Using WriteArrowBatch()" not in got_msg + + lyr = ds.GetLayer(0) + assert lyr.GetFeatureCount() == 1 + f = lyr.GetNextFeature() + assert f.GetGeometryRef().GetGeometryType() == ogr.wkbMultiLineString + assert f.GetGeometryRef().GetGeometryCount() == 2 From d0131e4684e55ba9e1e9abfc4f9755d2f335c970 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 2 Sep 2024 23:44:45 +0200 Subject: [PATCH 10/10] ogr2ogr: using std::async rather than std::thread --- apps/ogr2ogr_lib.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/apps/ogr2ogr_lib.cpp b/apps/ogr2ogr_lib.cpp index 5aafc7b22e3b..377d27ec4baa 100644 --- a/apps/ogr2ogr_lib.cpp +++ b/apps/ogr2ogr_lib.cpp @@ -41,6 +41,7 @@ #include #include +#include #include #include #include @@ -48,7 +49,6 @@ #include #include #include -#include #include #include @@ -5983,16 +5983,17 @@ bool LayerTranslator::TranslateArrow( if (nArrayLength >= MIN_FEATURES_FOR_THREADED_REPROJ && nNumReprojectionThreads >= 2) { - std::vector oThreads; + std::vector> oTasks; for (int iThread = 0; iThread < nNumReprojectionThreads; ++iThread) { - oThreads.emplace_back(oReprojectionLambda, iThread, - nNumReprojectionThreads); + oTasks.emplace_back(std::async(std::launch::async, + oReprojectionLambda, iThread, + nNumReprojectionThreads)); } - for (auto &oThread : oThreads) + for (auto &oTask : oTasks) { - oThread.join(); + oTask.get(); } } else