Skip to content

Commit

Permalink
Profile vertical transform fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@nyalldawson
nyalldawson committed Jul 11, 2024
1 parent c6d1be2 commit a95704b
Show file tree
Hide file tree
Showing 3 changed files with 197 additions and 2 deletions.
2 changes: 1 addition & 1 deletion src/core/qgsfeatureiterator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 & )
Expand Down
3 changes: 2 additions & 1 deletion src/core/vector/qgsvectorlayerprofilegenerator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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() );
Expand Down
194 changes: 194 additions & 0 deletions tests/src/python/test_qgsvectorlayerprofilegenerator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@
QgsSymbolLayer,
QgsVectorLayer,
QgsWkbTypes,
QgsPoint,
QgsCoordinateTransform
)
import unittest
from qgis.testing import start_app, QgisTestCase
Expand Down Expand Up @@ -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()

0 comments on commit a95704b

Please sign in to comment.