From b979610d2020b54804c71f2d74927af2c71befdf Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 7 Jan 2025 18:03:16 +0100 Subject: [PATCH] gml:CircleByCenterPoint: correctly take into account radius.uom for projected CRS Fixes #11597 --- autotest/ogr/ogr_gml_geom.py | 30 ++++++++++++++++ ogr/gml2ogrgeometry.cpp | 70 +++++++++++++++++++++--------------- 2 files changed, 71 insertions(+), 29 deletions(-) diff --git a/autotest/ogr/ogr_gml_geom.py b/autotest/ogr/ogr_gml_geom.py index 7a2a1ef55ada..086a1b635726 100755 --- a/autotest/ogr/ogr_gml_geom.py +++ b/autotest/ogr/ogr_gml_geom.py @@ -2295,6 +2295,36 @@ def test_gml_ArcByCenterPoint(): geom, ogr.CreateGeometryFromWkt("CIRCULARSTRING (1 4,-1 2,1 0)") ) + gml = "1 20.00290270" + geom = ogr.CreateGeometryFromGML(gml) + + ogrtest.check_feature_geometry( + geom, ogr.CreateGeometryFromWkt("CIRCULARSTRING (1 4,-1 2,1 0)") + ) + + gml = "1 20.00290270" + geom = ogr.CreateGeometryFromGML(gml) + + ogrtest.check_feature_geometry( + geom, ogr.CreateGeometryFromWkt("CIRCULARSTRING (1 4,-1 2,1 0)") + ) + + # EPSG:2222 : unit is ft + gml = "1 20.609690270" + geom = ogr.CreateGeometryFromGML(gml) + + ogrtest.check_feature_geometry( + geom, ogr.CreateGeometryFromWkt("CIRCULARSTRING (1 4,-1 2,1 0)") + ) + + # EPSG:2222 : unit is ft + gml = "1 2290270" + geom = ogr.CreateGeometryFromGML(gml) + + ogrtest.check_feature_geometry( + geom, ogr.CreateGeometryFromWkt("CIRCULARSTRING (1 4,-1 2,1 0)") + ) + ############################################################################### # Test compound curve of ArcByCenterPoint whose ends don't exactly match diff --git a/ogr/gml2ogrgeometry.cpp b/ogr/gml2ogrgeometry.cpp index 92c9716a490d..bd4a59adf386 100644 --- a/ogr/gml2ogrgeometry.cpp +++ b/ogr/gml2ogrgeometry.cpp @@ -811,26 +811,26 @@ static bool GML2OGRGeometry_AddToMultiSurface( } /************************************************************************/ -/* GetDistanceInMetre() */ +/* GetUOMInMetre() */ /************************************************************************/ -static double GetDistanceInMetre(double dfDistance, const char *pszUnits, - const char *pszAttribute, const char *pszId) +static double GetUOMInMetre(const char *pszUnits, const char *pszAttribute, + const char *pszId) { - if (EQUAL(pszUnits, "m")) - return dfDistance; + if (!pszUnits || EQUAL(pszUnits, "m")) + return 1.0; if (EQUAL(pszUnits, "km")) - return dfDistance * 1000; + return 1000.0; if (EQUAL(pszUnits, "nm") || EQUAL(pszUnits, "[nmi_i]")) - return dfDistance * CPLAtof(SRS_UL_INTL_NAUT_MILE_CONV); + return CPLAtof(SRS_UL_INTL_NAUT_MILE_CONV); if (EQUAL(pszUnits, "mi")) - return dfDistance * CPLAtof(SRS_UL_INTL_STAT_MILE_CONV); + return CPLAtof(SRS_UL_INTL_STAT_MILE_CONV); if (EQUAL(pszUnits, "ft")) - return dfDistance * CPLAtof(SRS_UL_INTL_FOOT_CONV); + return CPLAtof(SRS_UL_INTL_FOOT_CONV); if (pszId) { @@ -1211,8 +1211,12 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, const CPLXMLNode *psRadius = FindBareXMLChild(psChild, "radius"); if (psRadius && psRadius->eType == CXT_Element) { - double dfRadius = CPLAtof(CPLGetXMLValue(psRadius, nullptr, "0")); const char *pszUnits = CPLGetXMLValue(psRadius, "uom", nullptr); + const double dfUOMConv = GetUOMInMetre(pszUnits, "radius", pszId); + const double dfRadiusRaw = + CPLAtof(CPLGetXMLValue(psRadius, nullptr, "0")); + const double dfRadius = + dfUOMConv > 0 ? dfRadiusRaw * dfUOMConv : dfRadiusRaw; bool bSRSUnitIsDegree = false; bool bInvertedAxisOrder = false; if (l_pszSRSName != nullptr) @@ -1231,9 +1235,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, } } } - if (bSRSUnitIsDegree && pszUnits != nullptr && - (dfRadius = GetDistanceInMetre(dfRadius, pszUnits, "radius", - pszId)) > 0) + if (bSRSUnitIsDegree && dfUOMConv > 0) { bIsApproximateArc = true; dfLastCurveApproximateArcRadius = dfRadius; @@ -1808,8 +1810,11 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, ReportFailure("Missing radius element."); return nullptr; } - const double dfRadius = CPLAtof(CPLGetXMLValue(psChild, nullptr, "0")); const char *pszUnits = CPLGetXMLValue(psChild, "uom", nullptr); + const double dfUOMConv = GetUOMInMetre(pszUnits, "radius", pszId); + const double dfRadiusRaw = + CPLAtof(CPLGetXMLValue(psChild, nullptr, "0")); + double dfRadius = dfUOMConv > 0 ? dfRadiusRaw * dfUOMConv : dfRadiusRaw; psChild = FindBareXMLChild(psNode, "startAngle"); if (psChild == nullptr || psChild->eType != CXT_Element) @@ -1856,6 +1861,11 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, dfSemiMajor = GetSemiMajor(&oSRS); bInvertedAxisOrder = CPL_TO_BOOL(oSRS.EPSGTreatsAsNorthingEasting()); + + const double dfSRSUnitsToMetre = + oSRS.GetLinearUnits(nullptr); + if (dfSRSUnitsToMetre > 0) + dfRadius /= dfSRSUnitsToMetre; } } } @@ -1863,10 +1873,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, double dfCenterX = p.getX(); double dfCenterY = p.getY(); - double dfDistance; - if (bSRSUnitIsDegree && pszUnits != nullptr && - (dfDistance = - GetDistanceInMetre(dfRadius, pszUnits, "radius", pszId)) > 0) + if (bSRSUnitIsDegree && dfUOMConv > 0) { auto poLS = std::make_unique(); // coverity[tainted_data] @@ -1882,7 +1889,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, if (bInvertedAxisOrder) { OGR_GreatCircle_ExtendPosition( - dfCenterX, dfCenterY, dfDistance, + dfCenterX, dfCenterY, dfRadius, // See // https://ext.eurocontrol.int/aixm_confluence/display/ACG/ArcByCenterPoint+Interpretation+Summary dfAngle, dfSemiMajor, &dfLat, &dfLong); @@ -1892,7 +1899,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, else { OGR_GreatCircle_ExtendPosition( - dfCenterY, dfCenterX, dfDistance, 90 - dfAngle, + dfCenterY, dfCenterX, dfRadius, 90 - dfAngle, dfSemiMajor, &dfLat, &dfLong); p.setX(dfLong); p.setY(dfLat); @@ -1904,7 +1911,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, double dfLat = 0.0; if (bInvertedAxisOrder) { - OGR_GreatCircle_ExtendPosition(dfCenterX, dfCenterY, dfDistance, + OGR_GreatCircle_ExtendPosition(dfCenterX, dfCenterY, dfRadius, dfEndAngle, dfSemiMajor, &dfLat, &dfLong); p.setX(dfLat); // yes, external code will do the swap later @@ -1912,7 +1919,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, } else { - OGR_GreatCircle_ExtendPosition(dfCenterY, dfCenterX, dfDistance, + OGR_GreatCircle_ExtendPosition(dfCenterY, dfCenterX, dfRadius, 90 - dfEndAngle, dfSemiMajor, &dfLat, &dfLong); p.setX(dfLong); @@ -1955,8 +1962,11 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, ReportFailure("Missing radius element."); return nullptr; } - const double dfRadius = CPLAtof(CPLGetXMLValue(psChild, nullptr, "0")); const char *pszUnits = CPLGetXMLValue(psChild, "uom", nullptr); + const double dfUOMConv = GetUOMInMetre(pszUnits, "radius", pszId); + const double dfRadiusRaw = + CPLAtof(CPLGetXMLValue(psChild, nullptr, "0")); + double dfRadius = dfUOMConv > 0 ? dfRadiusRaw * dfUOMConv : dfRadiusRaw; OGRPoint p; if (!ParseGMLCoordinates(psNode, &p, nSRSDimension)) @@ -1985,6 +1995,11 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, dfSemiMajor = GetSemiMajor(&oSRS); bInvertedAxisOrder = CPL_TO_BOOL(oSRS.EPSGTreatsAsNorthingEasting()); + + const double dfSRSUnitsToMetre = + oSRS.GetLinearUnits(nullptr); + if (dfSRSUnitsToMetre > 0) + dfRadius /= dfSRSUnitsToMetre; } } } @@ -1992,10 +2007,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, double dfCenterX = p.getX(); double dfCenterY = p.getY(); - double dfDistance; - if (bSRSUnitIsDegree && pszUnits != nullptr && - (dfDistance = - GetDistanceInMetre(dfRadius, pszUnits, "radius", pszId)) > 0) + if (bSRSUnitIsDegree && dfUOMConv > 0) { auto poLS = std::make_unique(); const double dfStep = @@ -2007,7 +2019,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, if (bInvertedAxisOrder) { OGR_GreatCircle_ExtendPosition( - dfCenterX, dfCenterY, dfDistance, dfAngle, dfSemiMajor, + dfCenterX, dfCenterY, dfRadius, dfAngle, dfSemiMajor, &dfLat, &dfLong); p.setX(dfLat); // yes, external code will do the swap later p.setY(dfLong); @@ -2015,7 +2027,7 @@ GML2OGRGeometry_XMLNode_Internal(const CPLXMLNode *psNode, const char *pszId, else { OGR_GreatCircle_ExtendPosition( - dfCenterY, dfCenterX, dfDistance, dfAngle, dfSemiMajor, + dfCenterY, dfCenterX, dfRadius, dfAngle, dfSemiMajor, &dfLat, &dfLong); p.setX(dfLong); p.setY(dfLat);