Custom geoid support

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

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

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

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

app/position/providers/androidpositionprovider.cpp

@@ -56,16 +56,21 @@ void jniOnPositionUpdated( JNIEnv *env, jclass clazz, jint instanceId, jobject l
     const jdouble ellipsoidHeight = location.callMethod<jdouble>( "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;
+      }
     }
   }
 

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 );
   }
 }

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

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<double>::quiet_NaN();
+      position.elevation_diff = std::numeric_limits<double>::quiet_NaN();
+    }
   }
   else
   {
-    position.elevation = qQNaN();
-    position.elevation_diff = qQNaN();
+    position.elevation = std::numeric_limits<double>::quiet_NaN();
+    position.elevation_diff = std::numeric_limits<double>::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<double>::quiet_NaN();
+    position.elevation_diff = std::numeric_limits<double>::quiet_NaN();
+  }
   position.utcDateTime = QDateTime::currentDateTime();
   position.direction = 360 - static_cast<int>( mAngle ) % 360;
   position.hacc = ( *mGenerator )() % 20;

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 {