Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ensure vector layer elevation profiles respect layer/map vert datums #58148

Merged
merged 1 commit into from
Jul 23, 2024
Merged

Ensure vector layer elevation profiles respect layer/map vert datums #58148

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions src/app/elevation/qgselevationprofilewidget.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -479,6 +479,7 @@ QgsElevationProfileWidget::QgsElevationProfileWidget( const QString &name )
mLayerTreeView->populateInitialLayers( QgsProject::instance() );

connect( QgsProject::instance()->elevationProperties(), &QgsProjectElevationProperties::changed, this, &QgsElevationProfileWidget::onProjectElevationPropertiesChanged );
connect( QgsProject::instance(), &QgsProject::crs3DChanged, this, &QgsElevationProfileWidget::onProjectElevationPropertiesChanged );

updateCanvasLayers();
}
Expand Down Expand Up @@ -691,7 +692,7 @@ void QgsElevationProfileWidget::onCanvasPointHovered( const QgsPointXY &, const
void QgsElevationProfileWidget::updatePlot()
{
mCanvas->setTolerance( mSettingsAction->toleranceSpinBox()->value() );
mCanvas->setCrs( mMainCanvas->mapSettings().destinationCrs() );
mCanvas->setCrs( QgsProject::instance()->crs3D() );

if ( !mProfileCurve.isEmpty() )
{
Expand Down Expand Up @@ -880,7 +881,7 @@ void QgsElevationProfileWidget::exportResults( Qgis::ProfileExportType type )
file = QgsFileUtils::ensureFileNameHasExtension( file, QgsFileUtils::extensionsFromFilter( selectedFilter ) );

QgsProfileRequest request( profileCurve.release() );
request.setCrs( mMainCanvas->mapSettings().destinationCrs() );
request.setCrs( QgsProject::instance()->crs3D() );
request.setTolerance( mSettingsAction->toleranceSpinBox()->value() );
request.setTransformContext( QgsProject::instance()->transformContext() );
request.setTerrainProvider( QgsProject::instance()->elevationProperties()->terrainProvider() ? QgsProject::instance()->elevationProperties()->terrainProvider()->clone() : nullptr );
Expand Down
4 changes: 2 additions & 2 deletions src/core/vector/qgsvectorlayerprofilegenerator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -673,7 +673,7 @@ QgsVectorLayerProfileGenerator::QgsVectorLayerProfileGenerator( QgsVectorLayer *
, mProfileCurve( request.profileCurve() ? request.profileCurve()->clone() : nullptr )
, mTerrainProvider( request.terrainProvider() ? request.terrainProvider()->clone() : nullptr )
, mTolerance( request.tolerance() )
, mSourceCrs( layer->crs() )
, mSourceCrs( layer->crs3D() )
, mTargetCrs( request.crs() )
, mTransformContext( request.transformContext() )
, mExtent( layer->extent() )
Expand Down Expand Up @@ -806,7 +806,7 @@ bool QgsVectorLayerProfileGenerator::generateProfileForPoints()
{
// get features from layer
QgsFeatureRequest request;
request.setDestinationCrs( mTargetCrs, mTransformContext );
request.setCoordinateTransform( QgsCoordinateTransform( mSourceCrs, mTargetCrs, mTransformContext ) );
request.setDistanceWithin( QgsGeometry( mProfileCurve->clone() ), mTolerance );
request.setSubsetOfAttributes( mDataDefinedProperties.referencedFields( mExpressionContext ), mFields );
request.setFeedback( mFeedback.get() );
Expand Down
232 changes: 232 additions & 0 deletions tests/src/python/test_qgsvectorlayerprofilegenerator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,9 @@
QgsSymbolLayer,
QgsVectorLayer,
QgsWkbTypes,
QgsPoint,
QgsCoordinateTransform,
QgsDatumTransform
)
import unittest
from qgis.testing import start_app, QgisTestCase
Expand Down Expand Up @@ -2309,6 +2312,235 @@ 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 test_vertical_transformation_4978_to_4985(self):
"""
Test vertical transformations are correctly handled during profile generation

EPSG:4979 to EPSG:4985
"""

vl = QgsVectorLayer('PointZ?crs=EPSG:4979', '4979_points', 'memory')
self.assertTrue(vl.isValid())
self.assertEqual(vl.crs().authid(), 'EPSG:4979')

self.assertEqual(vl.crs3D().horizontalCrs().authid(), 'EPSG:4979')

f = QgsFeature()
f.setGeometry(QgsPoint(134.445567853,
-23.445567853,
5543.325))
self.assertTrue(vl.dataProvider().addFeature(f))

profile_crs = QgsCoordinateReferenceSystem('EPSG:4985')
self.assertTrue(profile_crs.isValid())
self.assertEqual(profile_crs.horizontalCrs().authid(), 'EPSG:4985')

curve = QgsLineString()
curve.fromWkt(
'LineString (134.405567853 -23.435567853, 134.485567853 -23.455567853)')
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.4454139641, 6)
self.assertAlmostEqual(res.constGet().y(), -23.4456037763, 6)
self.assertAlmostEqual(res.constGet().z(), 5545.6857, 3)

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')

available_operations = QgsDatumTransform.operations(vl.crs3D(),
profile_crs)
self.assertEqual(len(available_operations[0].grids), 1)
self.assertEqual(available_operations[0].grids[0].shortName,
'au_ga_AGQG_20201120.tif')
if not available_operations[0].isAvailable:
self.skipTest(
f'Required grid {available_operations[0].grids[0].shortName} not available on system')

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')

available_operations = QgsDatumTransform.operations(vl.crs3D(),
profile_crs)
self.assertEqual(len(available_operations[0].grids), 1)
self.assertEqual(available_operations[0].grids[0].shortName, 'au_ga_AGQG_20201120.tif')
if not available_operations[0].isAvailable:
self.skipTest(f'Required grid {available_operations[0].grids[0].shortName} not available on system')

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')

available_operations = QgsDatumTransform.operations(vl.crs3D(),
profile_crs)
self.assertEqual(len(available_operations[0].grids), 1)
self.assertEqual(available_operations[0].grids[0].shortName, 'au_ga_AUSGeoid2020_20180201.tif')
if not available_operations[0].isAvailable:
self.skipTest(f'Required grid {available_operations[0].grids[0].shortName} not available on system')

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')

available_operations = QgsDatumTransform.operations(vl.crs3D(),
profile_crs)
self.assertEqual(len(available_operations[0].grids), 1)
self.assertEqual(available_operations[0].grids[0].shortName, 'au_ga_AUSGeoid2020_20180201.tif')
if not available_operations[0].isAvailable:
self.skipTest(f'Required grid {available_operations[0].grids[0].shortName} not available on system')

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()
Loading