diff --git a/src/core/qgsfeatureiterator.cpp b/src/core/qgsfeatureiterator.cpp index e931b3140a67..952d92537be7 100644 --- a/src/core/qgsfeatureiterator.cpp +++ b/src/core/qgsfeatureiterator.cpp @@ -102,7 +102,7 @@ void QgsAbstractFeatureIterator::geometryToDestinationCrs( QgsFeature &feature, try { QgsGeometry g = feature.geometry(); - g.transform( transform ); + g.transform( transform, Qgis::TransformDirection::Forward, true ); feature.setGeometry( g ); } catch ( QgsCsException & ) diff --git a/src/core/vector/qgsvectorlayerprofilegenerator.cpp b/src/core/vector/qgsvectorlayerprofilegenerator.cpp index dc342f760da0..d27ab5023fc3 100644 --- a/src/core/vector/qgsvectorlayerprofilegenerator.cpp +++ b/src/core/vector/qgsvectorlayerprofilegenerator.cpp @@ -806,7 +806,8 @@ bool QgsVectorLayerProfileGenerator::generateProfileForPoints() { // get features from layer QgsFeatureRequest request; - request.setDestinationCrs( mTargetCrs, mTransformContext ); + request.setCoordinateTransform( QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext ) ); +// request.setDestinationCrs(mTargetCrs, mTransformContext ); request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance ); request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields ); request.setFeedback( mFeedback.get() ); diff --git a/tests/src/python/test_qgsvectorlayerprofilegenerator.py b/tests/src/python/test_qgsvectorlayerprofilegenerator.py index 2a68466da5e4..30df090d7ad2 100644 --- a/tests/src/python/test_qgsvectorlayerprofilegenerator.py +++ b/tests/src/python/test_qgsvectorlayerprofilegenerator.py @@ -39,6 +39,8 @@ QgsSymbolLayer, QgsVectorLayer, QgsWkbTypes, + QgsPoint, + QgsCoordinateTransform ) import unittest from qgis.testing import start_app, QgisTestCase @@ -2309,6 +2311,198 @@ def testPolygonGenerationFeature(self): self.doCheckLine(req, 10, vl, [168, 172, 206, 210, 231, 267, 275, 282, 283, 284, 306, 307, 319, 321], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], Qgis.GeometryType.Line) self.doCheckLine(req, 11, vl, [168, 172, 206, 210, 231, 237, 255, 267, 275, 282, 283, 284, 306, 307, 319, 321], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], Qgis.GeometryType.Line) + def testProfileTransformGDA2020toAVWS(self): + """ + Test a profile requiring a vertical datum transform from GDA2020 to AVWS + """ + # GDA2020 vertical CRS + vl = QgsVectorLayer('PointZ?crs=EPSG:7843', 'gda2020points', 'memory') + self.assertTrue(vl.isValid()) + self.assertEqual(vl.crs().authid(), 'EPSG:7843') + + self.assertEqual(vl.crs3D().horizontalCrs().authid(), 'EPSG:7843') + + f = QgsFeature() + f.setGeometry(QgsPoint(134.445567853, + -23.445567853, + 5543.325)) + self.assertTrue(vl.dataProvider().addFeature(f)) + + profile_crs, msg = QgsCoordinateReferenceSystem.createCompoundCrs( + QgsCoordinateReferenceSystem('EPSG:7844'), + QgsCoordinateReferenceSystem('EPSG:9458')) + self.assertFalse(msg) + self.assertTrue(profile_crs.isValid()) + self.assertEqual(profile_crs.horizontalCrs().authid(), 'EPSG:7844') + self.assertEqual(profile_crs.verticalCrs().authid(), 'EPSG:9458') + + # for debugging only + # transform = QgsCoordinateTransform( + # vl.crs3D(), + # profile_crs, + # QgsCoordinateTransformContext() + #) + # self.assertEqual(transform.instantiatedCoordinateOperationDetails().proj, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AGQG_20201120.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg') + + curve = QgsLineString() + curve.fromWkt( + 'LineString (134.445567853 -23.445567853, 135.445567853 -23.445567853)') + request = QgsProfileRequest(curve) + request.setCrs(profile_crs) + request.setTolerance(1) + + generator = vl.createProfileGenerator(request) + self.assertTrue(generator.generateProfile()) + + r = generator.takeResults() + res = r.asGeometries()[0] + self.assertAlmostEqual(res.constGet().x(), 134.445567853, 6) + self.assertAlmostEqual(res.constGet().y(), -23.445567853, 6) + # comparing against results from https://geodesyapps.ga.gov.au/avws + self.assertAlmostEqual(res.constGet().z(), 5524.13969, 3) + + def testProfileTransformAVWStoGDA2020(self): + """ + Test a profile requiring a AVWS vertical datum transform to GDA2020 + """ + # AVWS vertical CRS + vl = QgsVectorLayer('PointZ?crs=EPSG:7844', 'gda2020points', 'memory') + self.assertTrue(vl.isValid()) + self.assertEqual(vl.crs().authid(), 'EPSG:7844') + vl.setVerticalCrs(QgsCoordinateReferenceSystem('EPSG:9458')) + + self.assertEqual(vl.crs3D().horizontalCrs().authid(), 'EPSG:7844') + self.assertEqual(vl.crs3D().verticalCrs().authid(), 'EPSG:9458') + + f = QgsFeature() + f.setGeometry(QgsPoint(134.445567853, + -23.445567853, + 5524.13969)) + self.assertTrue(vl.dataProvider().addFeature(f)) + + profile_crs = QgsCoordinateReferenceSystem('EPSG:7843') + + # for debugging only + # transform = QgsCoordinateTransform( + # vl.crs3D(), + # profile_crs, + # QgsCoordinateTransformContext() + #) + # self.assertEqual(transform.instantiatedCoordinateOperationDetails().proj, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AGQG_20201120.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg') + + curve = QgsLineString() + curve.fromWkt( + 'LineString (134.445567853 -23.445567853, 135.445567853 -23.445567853)') + request = QgsProfileRequest(curve) + request.setCrs(profile_crs) + request.setTolerance(1) + + generator = vl.createProfileGenerator(request) + self.assertTrue(generator.generateProfile()) + + r = generator.takeResults() + res = r.asGeometries()[0] + self.assertAlmostEqual(res.constGet().x(), 134.445567853, 6) + self.assertAlmostEqual(res.constGet().y(), -23.445567853, 6) + # comparing against results from https://geodesyapps.ga.gov.au/avws + self.assertAlmostEqual(res.constGet().z(), 5543.325, 3) + + def testProfileTransformGDA2020toAHD(self): + """ + Test a profile requiring a vertical datum transform from GDA2020 to AHD + """ + # GDA2020 vertical CRS + vl = QgsVectorLayer('PointZ?crs=EPSG:7843', 'gda2020points', 'memory') + self.assertTrue(vl.isValid()) + self.assertEqual(vl.crs().authid(), 'EPSG:7843') + + self.assertEqual(vl.crs3D().horizontalCrs().authid(), 'EPSG:7843') + + f = QgsFeature() + f.setGeometry(QgsPoint(134.445567853, + -23.445567853, + 5543.325)) + self.assertTrue(vl.dataProvider().addFeature(f)) + + profile_crs, msg = QgsCoordinateReferenceSystem.createCompoundCrs( + QgsCoordinateReferenceSystem('EPSG:7844'), + QgsCoordinateReferenceSystem('EPSG:5711')) + self.assertFalse(msg) + self.assertTrue(profile_crs.isValid()) + self.assertEqual(profile_crs.horizontalCrs().authid(), 'EPSG:7844') + self.assertEqual(profile_crs.verticalCrs().authid(), 'EPSG:5711') + + # for debugging only + # transform = QgsCoordinateTransform( + # vl.crs3D(), + # profile_crs, + # QgsCoordinateTransformContext() + #) + # self.assertEqual(transform.instantiatedCoordinateOperationDetails().proj, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AGQG_20201120.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg') + + curve = QgsLineString() + curve.fromWkt( + 'LineString (134.445567853 -23.445567853, 135.445567853 -23.445567853)') + request = QgsProfileRequest(curve) + request.setCrs(profile_crs) + request.setTolerance(1) + + generator = vl.createProfileGenerator(request) + self.assertTrue(generator.generateProfile()) + + r = generator.takeResults() + res = r.asGeometries()[0] + self.assertAlmostEqual(res.constGet().x(), 134.445567853, 6) + self.assertAlmostEqual(res.constGet().y(), -23.445567853, 6) + # comparing against results from https://geodesyapps.ga.gov.au/ausgeoid2020 + self.assertAlmostEqual(res.constGet().z(), 5523.598, 3) + + def testProfileTransformAHDtoGDA2020(self): + """ + Test a profile requiring a AHD vertical datum transform to GDA2020 + """ + # AHD vertical CRS + vl = QgsVectorLayer('PointZ?crs=EPSG:7844', 'gda2020points', 'memory') + self.assertTrue(vl.isValid()) + self.assertEqual(vl.crs().authid(), 'EPSG:7844') + vl.setVerticalCrs(QgsCoordinateReferenceSystem('EPSG:5711')) + + self.assertEqual(vl.crs3D().horizontalCrs().authid(), 'EPSG:7844') + self.assertEqual(vl.crs3D().verticalCrs().authid(), 'EPSG:5711') + + f = QgsFeature() + f.setGeometry(QgsPoint(134.445567853, + -23.445567853, + 5523.598)) + self.assertTrue(vl.dataProvider().addFeature(f)) + + profile_crs = QgsCoordinateReferenceSystem('EPSG:7843') + + # for debugging only + # transform = QgsCoordinateTransform( + # vl.crs3D(), + # profile_crs, + # QgsCoordinateTransformContext() + #) + # self.assertEqual(transform.instantiatedCoordinateOperationDetails().proj, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AGQG_20201120.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg') + + curve = QgsLineString() + curve.fromWkt( + 'LineString (134.445567853 -23.445567853, 135.445567853 -23.445567853)') + request = QgsProfileRequest(curve) + request.setCrs(profile_crs) + request.setTolerance(1) + + generator = vl.createProfileGenerator(request) + self.assertTrue(generator.generateProfile()) + + r = generator.takeResults() + res = r.asGeometries()[0] + self.assertAlmostEqual(res.constGet().x(), 134.445567853, 6) + self.assertAlmostEqual(res.constGet().y(), -23.445567853, 6) + # comparing against results from https://geodesyapps.ga.gov.au/ausgeoid2020 + self.assertAlmostEqual(res.constGet().z(), 5543.325, 3) + if __name__ == '__main__': unittest.main()