diff --git a/app/inpututils.cpp b/app/inpututils.cpp index bed368ff3..146f9871e 100644 --- a/app/inpututils.cpp +++ b/app/inpututils.cpp @@ -905,6 +905,42 @@ QgsPoint InputUtils::transformPoint( const QgsCoordinateReferenceSystem &srcCrs, return {}; } +QgsPoint InputUtils::transformPoint( const QgsCoordinateReferenceSystem &srcCrs, + const QgsCoordinateReferenceSystem &destCrs, + const QgsCoordinateTransformContext &context, + const QgsPoint &srcPoint, bool &fallbackOperationOccurred ) +{ + // we do not want to transform empty points, + // QGIS would convert them to a valid (0, 0) points + if ( srcPoint.isEmpty() ) + { + return {}; + } + + try + { + const QgsCoordinateTransform ct( srcCrs, destCrs, context ); + if ( ct.isValid() ) + { + if ( !ct.isShortCircuited() ) + { + const QgsVector3D transformed = ct.transform( QgsVector3D( srcPoint.x(), srcPoint.y(), srcPoint.z() ) ); + fallbackOperationOccurred = ct.fallbackOperationOccurred(); + const QgsPoint pt( transformed.x(), transformed.y(), transformed.z(), srcPoint.m() ); + return pt; + } + + return srcPoint; + } + } + catch ( QgsCsException &cse ) + { + Q_UNUSED( cse ) + } + + return {}; +} + QPointF InputUtils::transformPointToScreenCoordinates( const QgsCoordinateReferenceSystem &srcCrs, InputMapSettings *mapSettings, const QgsPoint &srcPoint ) { if ( !mapSettings || srcPoint.isEmpty() ) diff --git a/app/inpututils.h b/app/inpututils.h index a31148d3a..33db41997 100644 --- a/app/inpututils.h +++ b/app/inpututils.h @@ -305,6 +305,14 @@ class InputUtils: public QObject const QgsCoordinateTransformContext &context, const QgsPoint &srcPoint ); + /** + * Overload of transformPoint function, which also notifies the caller if PROJ fallback operation occurred + */ + static QgsPoint transformPoint( const QgsCoordinateReferenceSystem &srcCrs, + const QgsCoordinateReferenceSystem &destCrs, + const QgsCoordinateTransformContext &context, + const QgsPoint &srcPoint, bool &fallbackOperationOccurred ); + /** * Transforms point between CRS and screen pixels * Return empty QgsPoint if the transformation could not be applied or srcPoint is empty diff --git a/app/position/positionkit.cpp b/app/position/positionkit.cpp index f44af4140..8c98590f0 100644 --- a/app/position/positionkit.cpp +++ b/app/position/positionkit.cpp @@ -35,9 +35,29 @@ PositionKit::PositionKit( QObject *parent ) QgsCoordinateReferenceSystem PositionKit::positionCrs3D() { + bool crsExists = false; + const QString crsWktDef = QgsProject::instance()->readEntry( QStringLiteral( "Mergin" ), QStringLiteral( "TargetVerticalCRS" ), QString(), &crsExists ); + if ( crsExists ) + { + const QgsCoordinateReferenceSystem verticalCrs = QgsCoordinateReferenceSystem::fromWkt( crsWktDef ); + QString compoundCrsError{}; + const QgsCoordinateReferenceSystem compoundCrs = QgsCoordinateReferenceSystem::createCompoundCrs( positionCrs2D(), verticalCrs, compoundCrsError ); + if ( compoundCrs.isValid() && compoundCrsError.isEmpty() ) + { + return compoundCrs; + } + CoreUtils::log( QStringLiteral( "PositionKit" ), QStringLiteral( "Failed to create custom compound crs: %1" ).arg( compoundCrsError ) ); + } + return QgsCoordinateReferenceSystem::fromEpsgId( 9707 ); } +QString PositionKit::positionCrs3DGeoidModelName() +{ + const QgsCoordinateReferenceSystem crs = positionCrs3D().verticalCrs(); + return crs.description(); +} + QgsCoordinateReferenceSystem PositionKit::positionCrs2D() { return QgsCoordinateReferenceSystem::fromEpsgId( 4326 ); diff --git a/app/position/positionkit.h b/app/position/positionkit.h index d31596bda..0d21e339b 100644 --- a/app/position/positionkit.h +++ b/app/position/positionkit.h @@ -124,8 +124,13 @@ class PositionKit : public QObject static QgsCoordinateReferenceSystem positionCrs2D(); // Coordinate reference system - WGS84 + ellipsoid height (EPSG:4979) static QgsCoordinateReferenceSystem positionCrs3DEllipsoidHeight(); - // Coordinate reference system of position - WGS84 + geoid height - egm96_15 (EPSG:9707) + /** + * Coordinate reference system of position (WGS84 + geoid height) - can use custom geoid model + * \note by default we use egm96_15 model (EPSG:9707) + */ static QgsCoordinateReferenceSystem positionCrs3D(); + // Returns the model name used for elevation transformations + Q_INVOKABLE static QString positionCrs3DGeoidModelName(); Q_INVOKABLE static AbstractPositionProvider *constructProvider( const QString &type, const QString &id, const QString &name = QString() ); Q_INVOKABLE static AbstractPositionProvider *constructActiveProvider( AppSettings *appsettings ); diff --git a/app/position/providers/androidpositionprovider.cpp b/app/position/providers/androidpositionprovider.cpp index 6581052b5..9de51ae2e 100644 --- a/app/position/providers/androidpositionprovider.cpp +++ b/app/position/providers/androidpositionprovider.cpp @@ -56,16 +56,21 @@ void jniOnPositionUpdated( JNIEnv *env, jclass clazz, jint instanceId, jobject l const jdouble ellipsoidHeight = location.callMethod( "getAltitude" ); if ( !qFuzzyIsNull( ellipsoidHeight ) ) { + bool positionOutsideGeoidModelArea = false; // transform the altitude from EPSG:4979 (WGS84 (EPSG:4326) + ellipsoidal height) to specified geoid model const QgsPoint geoidPosition = InputUtils::transformPoint( PositionKit::positionCrs3DEllipsoidHeight(), PositionKit::positionCrs3D(), QgsProject::instance()->transformContext(), - {longitude, latitude, ellipsoidHeight} ); - pos.elevation = geoidPosition.z(); - - const double geoidSeparation = ellipsoidHeight - geoidPosition.z(); - pos.elevation_diff = geoidSeparation; + {longitude, latitude, ellipsoidHeight}, + positionOutsideGeoidModelArea ); + if ( !positionOutsideGeoidModelArea ) + { + pos.elevation = geoidPosition.z(); + + const double geoidSeparation = ellipsoidHeight - geoidPosition.z(); + pos.elevation_diff = geoidSeparation; + } } } diff --git a/app/position/providers/bluetoothpositionprovider.cpp b/app/position/providers/bluetoothpositionprovider.cpp index 566347d10..3fdcc7b47 100644 --- a/app/position/providers/bluetoothpositionprovider.cpp +++ b/app/position/providers/bluetoothpositionprovider.cpp @@ -216,14 +216,18 @@ void BluetoothPositionProvider::positionUpdateReceived() // The geoid models used in GNSS devices can be often times unreliable, thus we apply the transformations ourselves // GNSS supplied orthometric elevation -> ellipsoid elevation -> orthometric elevation based on our model const double ellipsoidElevation = positionData.elevation + positionData.elevation_diff; + bool positionOutsideGeoidModelArea = false; const QgsPoint geoidPosition = InputUtils::transformPoint( PositionKit::positionCrs3DEllipsoidHeight(), PositionKit::positionCrs3D(), QgsProject::instance()->transformContext(), - {positionData.longitude, positionData.latitude, ellipsoidElevation} ); - positionData.elevation = geoidPosition.z(); - positionData.elevation_diff = ellipsoidElevation - geoidPosition.z(); - + {positionData.longitude, positionData.latitude, ellipsoidElevation}, + positionOutsideGeoidModelArea ); + if ( !positionOutsideGeoidModelArea ) + { + positionData.elevation = geoidPosition.z(); + positionData.elevation_diff = ellipsoidElevation - geoidPosition.z(); + } emit positionChanged( positionData ); } } diff --git a/app/position/providers/internalpositionprovider.cpp b/app/position/providers/internalpositionprovider.cpp index d7758aba1..b036c32c6 100644 --- a/app/position/providers/internalpositionprovider.cpp +++ b/app/position/providers/internalpositionprovider.cpp @@ -143,29 +143,34 @@ void InternalPositionProvider::parsePositionUpdate( const QGeoPositionInfo &posi } // transform the altitude from EPSG:4979 (WGS84 (EPSG:4326) + ellipsoidal height) to specified geoid model + bool positionOutsideGeoidModelArea = false; const QgsPoint geoidPosition = InputUtils::transformPoint( PositionKit::positionCrs3DEllipsoidHeight(), PositionKit::positionCrs3D(), QgsProject::instance()->transformContext(), - {position.coordinate().longitude(), position.coordinate().latitude(), position.coordinate().altitude()} ); - if ( !qgsDoubleNear( geoidPosition.z(), mLastPosition.elevation ) ) + {position.coordinate().longitude(), position.coordinate().latitude(), position.coordinate().altitude()}, + positionOutsideGeoidModelArea ); + if ( !positionOutsideGeoidModelArea ) { - mLastPosition.elevation = geoidPosition.z(); - positionDataHasChanged = true; - } - - // QGeoCoordinate::altitude() docs claim that it is above the sea level (i.e. geoid) altitude, - // but that's not really true in our case: - // - on Android - it is MSL altitude only if "useMslAltitude" parameter is passed to the Android - // Qt positioning plugin, which we don't do - see https://doc.qt.io/qt-6/position-plugin-android.html - // - on iOS - it would return MSL altitude, but we have a custom patch in vcpkg to return - // ellipsoid altitude (so we do not rely on geoid model of unknown quality/resolution) - const double ellipsoidAltitude = position.coordinate().altitude(); - const double geoidSeparation = ellipsoidAltitude - geoidPosition.z(); - if ( !qgsDoubleNear( geoidSeparation, mLastPosition.elevation_diff ) ) - { - mLastPosition.elevation_diff = geoidSeparation; - positionDataHasChanged = true; + if ( !qgsDoubleNear( geoidPosition.z(), mLastPosition.elevation ) ) + { + mLastPosition.elevation = geoidPosition.z(); + positionDataHasChanged = true; + } + + // QGeoCoordinate::altitude() docs claim that it is above the sea level (i.e. geoid) altitude, + // but that's not really true in our case: + // - on Android - it is MSL altitude only if "useMslAltitude" parameter is passed to the Android + // Qt positioning plugin, which we don't do - see https://doc.qt.io/qt-6/position-plugin-android.html + // - on iOS - it would return MSL altitude, but we have a custom patch in vcpkg to return + // ellipsoid altitude (so we do not rely on geoid model of unknown quality/resolution) + const double ellipsoidAltitude = position.coordinate().altitude(); + const double geoidSeparation = ellipsoidAltitude - geoidPosition.z(); + if ( !qgsDoubleNear( geoidSeparation, mLastPosition.elevation_diff ) ) + { + mLastPosition.elevation_diff = geoidSeparation; + positionDataHasChanged = true; + } } bool hasSpeedInfo = position.hasAttribute( QGeoPositionInfo::GroundSpeed ); diff --git a/app/position/providers/simulatedpositionprovider.cpp b/app/position/providers/simulatedpositionprovider.cpp index 762bbc416..994fe0e89 100644 --- a/app/position/providers/simulatedpositionprovider.cpp +++ b/app/position/providers/simulatedpositionprovider.cpp @@ -92,18 +92,28 @@ void SimulatedPositionProvider::generateRadiusPosition() double ellipsoidAltitude = ( *mGenerator )() % 40 + 80; // rand altitude <80,115>m and lost (NaN) if ( ellipsoidAltitude <= 115 ) { + bool positionOutsideGeoidModelArea = false; const QgsPoint geoidPosition = InputUtils::transformPoint( PositionKit::positionCrs3DEllipsoidHeight(), PositionKit::positionCrs3D(), QgsCoordinateTransformContext(), - {longitude, latitude, ellipsoidAltitude} ); - position.elevation = geoidPosition.z(); - position.elevation_diff = ellipsoidAltitude - position.elevation; + {longitude, latitude, ellipsoidAltitude}, + positionOutsideGeoidModelArea ); + if ( !positionOutsideGeoidModelArea ) + { + position.elevation = geoidPosition.z(); + position.elevation_diff = ellipsoidAltitude - position.elevation; + } + else + { + position.elevation = std::numeric_limits::quiet_NaN(); + position.elevation_diff = std::numeric_limits::quiet_NaN(); + } } else { - position.elevation = qQNaN(); - position.elevation_diff = qQNaN(); + position.elevation = std::numeric_limits::quiet_NaN(); + position.elevation_diff = std::numeric_limits::quiet_NaN(); } const QDateTime timestamp = QDateTime::currentDateTime(); @@ -132,13 +142,23 @@ void SimulatedPositionProvider::generateConstantPosition() position.latitude = mLatitude; position.longitude = mLongitude; // we take 100 as elevation returned by WGS84 ellipsoid and recalculate it to geoid + bool positionOutsideGeoidModelArea = false; const QgsPoint geoidPosition = InputUtils::transformPoint( PositionKit::positionCrs3DEllipsoidHeight(), PositionKit::positionCrs3D(), QgsCoordinateTransformContext(), - {mLongitude, mLatitude, 100} ); - position.elevation = geoidPosition.z(); - position.elevation_diff = 100 - position.elevation; + {mLongitude, mLatitude, 100}, + positionOutsideGeoidModelArea ); + if ( !positionOutsideGeoidModelArea ) + { + position.elevation = geoidPosition.z(); + position.elevation_diff = 100 - position.elevation; + } + else + { + position.elevation = std::numeric_limits::quiet_NaN(); + position.elevation_diff = std::numeric_limits::quiet_NaN(); + } position.utcDateTime = QDateTime::currentDateTime(); position.direction = 360 - static_cast( mAngle ) % 360; position.hacc = ( *mGenerator )() % 20; diff --git a/app/qml/gps/MMGpsDataDrawer.qml b/app/qml/gps/MMGpsDataDrawer.qml index 35a36fc3e..efe5ee7ad 100644 --- a/app/qml/gps/MMGpsDataDrawer.qml +++ b/app/qml/gps/MMGpsDataDrawer.qml @@ -208,7 +208,7 @@ MMComponents.MMDrawer { } alignmentRight: Positioner.index % 2 === 1 - desc: qsTr("Orthometric height, using EGM96 geoid") + desc: qsTr(("Orthometric height, using %1 geoid").arg(PositionKit.positionCrs3DGeoidModelName())) } MMGpsComponents.MMGpsDataText {