qgsellipsoidutils.cpp

Go to the documentation of this file.

/***************************************************************************
  qgsellipsoidutils.cpp
 ----------------------
  Date                 : April 2017
  Copyright            : (C) 2017 by Nyall Dawson
  email                : nyall dot dawson at gmail dot com
 ***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgsellipsoidutils.h"
 
#include <mutex>
#include <proj.h>
#include <sqlite3.h>
 
#include "qgsapplication.h"
#include "qgscelestialbody.h"
#include "qgscoordinatereferencesystemregistry.h"
#include "qgslogger.h"
#include "qgsmessagelog.h"
#include "qgsprojutils.h"
#include "qgsreadwritelocker.h"
#include "qgsruntimeprofiler.h"
 
#include <QCollator>
#include <QMatrix3x3>
#include <QString>
 
using namespace Qt::StringLiterals;
 
Q_GLOBAL_STATIC( QReadWriteLock, sEllipsoidCacheLock )
typedef QHash< QString, QgsEllipsoidUtils::EllipsoidParameters > EllipsoidParamCache;
Q_GLOBAL_STATIC( EllipsoidParamCache, sEllipsoidCache )
 
Q_GLOBAL_STATIC( QReadWriteLock, sDefinitionCacheLock );
typedef QList< QgsEllipsoidUtils::EllipsoidDefinition > EllipsoidDefinitionCache;
Q_GLOBAL_STATIC( EllipsoidDefinitionCache, sDefinitionCache )
 
static bool sDisableCache = false;
 
QgsEllipsoidUtils::EllipsoidParameters QgsEllipsoidUtils::ellipsoidParameters( const QString &e )
{
  // maps older QGIS ellipsoid acronyms to proj acronyms/names
  static const QMap< QString, QString > sProj6EllipsoidAcronymMap {
    { "clrk80", "clrk80ign" },
    { "Adrastea2000", "ESRI:107909" },
    { "Amalthea2000", "ESRI:107910" },
    { "Ananke2000", "ESRI:107911" },
    { "Ariel2000", "ESRI:107945" },
    { "Atlas2000", "ESRI:107926" },
    { "Belinda2000", "ESRI:107946" },
    { "Bianca2000", "ESRI:107947" },
    { "Callisto2000", "ESRI:107912" },
    { "Calypso2000", "ESRI:107927" },
    { "Carme2000", "ESRI:107913" },
    { "Charon2000", "ESRI:107970" },
    { "Cordelia2000", "ESRI:107948" },
    { "Cressida2000", "ESRI:107949" },
    { "Deimos2000", "ESRI:107906" },
    { "Desdemona2000", "ESRI:107950" },
    { "Despina2000", "ESRI:107961" },
    { "Dione2000", "ESRI:107928" },
    { "Elara2000", "ESRI:107914" },
    { "Enceladus2000", "ESRI:107929" },
    { "Epimetheus2000", "ESRI:107930" },
    { "Europa2000", "ESRI:107915" },
    { "Galatea2000", "ESRI:107962" },
    { "Ganymede2000", "ESRI:107916" },
    { "Helene2000", "ESRI:107931" },
    { "Himalia2000", "ESRI:107917" },
    { "Hyperion2000", "ESRI:107932" },
    { "Iapetus2000", "ESRI:107933" },
    { "Io2000", "ESRI:107918" },
    { "Janus2000", "ESRI:107934" },
    { "Juliet2000", "ESRI:107951" },
    { "Jupiter2000", "ESRI:107908" },
    { "Larissa2000", "ESRI:107963" },
    { "Leda2000", "ESRI:107919" },
    { "Lysithea2000", "ESRI:107920" },
    { "Mars2000", "ESRI:107905" },
    { "Mercury2000", "ESRI:107900" },
    { "Metis2000", "ESRI:107921" },
    { "Mimas2000", "ESRI:107935" },
    { "Miranda2000", "ESRI:107952" },
    { "Moon2000", "ESRI:107903" },
    { "Naiad2000", "ESRI:107964" },
    { "Neptune2000", "ESRI:107960" },
    { "Nereid2000", "ESRI:107965" },
    { "Oberon2000", "ESRI:107953" },
    { "Ophelia2000", "ESRI:107954" },
    { "Pan2000", "ESRI:107936" },
    { "Pandora2000", "ESRI:107937" },
    { "Pasiphae2000", "ESRI:107922" },
    { "Phobos2000", "ESRI:107907" },
    { "Phoebe2000", "ESRI:107938" },
    { "Pluto2000", "ESRI:107969" },
    { "Portia2000", "ESRI:107955" },
    { "Prometheus2000", "ESRI:107939" },
    { "Proteus2000", "ESRI:107966" },
    { "Puck2000", "ESRI:107956" },
    { "Rhea2000", "ESRI:107940" },
    { "Rosalind2000", "ESRI:107957" },
    { "Saturn2000", "ESRI:107925" },
    { "Sinope2000", "ESRI:107923" },
    { "Telesto2000", "ESRI:107941" },
    { "Tethys2000", "ESRI:107942" },
    { "Thalassa2000", "ESRI:107967" },
    { "Thebe2000", "ESRI:107924" },
    { "Titan2000", "ESRI:107943" },
    { "Titania2000", "ESRI:107958" },
    { "Triton2000", "ESRI:107968" },
    { "Umbriel2000", "ESRI:107959" },
    { "Uranus2000", "ESRI:107944" },
    { "Venus2000", "ESRI:107902" },
    { "IGNF:ELG053", "EPSG:7030" },
    { "IGNF:ELG052", "EPSG:7043" },
    { "IGNF:ELG102", "EPSG:7043" },
    { "WGS66", "ESRI:107001" },
    { "plessis", "EPSG:7027" },
    { "IGNF:ELG017", "EPSG:7027" },
    { "mod_airy", "EPSG:7002" },
    { "IGNF:ELG037", "EPSG:7019" },
    { "IGNF:ELG108", "EPSG:7036" },
    { "cape", "EPSG:7034" },
    { "IGNF:ELG010", "EPSG:7011" },
    { "IGNF:ELG003", "EPSG:7012" },
    { "IGNF:ELG004", "EPSG:7008" },
    { "GSK2011", "EPSG:1025" },
    { "airy", "EPSG:7001" },
    { "aust_SA", "EPSG:7003" },
    { "bessel", "EPSG:7004" },
    { "clrk66", "EPSG:7008" },
    { "clrk80ign", "EPSG:7011" },
    { "evrst30", "EPSG:7015" },
    { "evrstSS", "EPSG:7016" },
    { "evrst48", "EPSG:7018" },
    { "GRS80", "EPSG:7019" },
    { "helmert", "EPSG:7020" },
    { "intl", "EPSG:7022" },
    { "krass", "EPSG:7024" },
    { "NWL9D", "EPSG:7025" },
    { "WGS84", "EPSG:7030" },
    { "GRS67", "EPSG:7036" },
    { "WGS72", "EPSG:7043" },
    { "bess_nam", "EPSG:7046" },
    { "IAU76", "EPSG:7049" },
    { "sphere", "EPSG:7052" },
    { "hough", "EPSG:7053" },
    { "evrst69", "EPSG:7056" },
    { "fschr60", "ESRI:107002" },
    { "fschr68", "ESRI:107003" },
    { "fschr60m", "ESRI:107004" },
    { "walbeck", "ESRI:107007" },
    { "IGNF:ELG001", "EPSG:7022" },
    { "engelis", "EPSG:7054" },
    { "evrst56", "EPSG:7044" },
    { "SEasia", "ESRI:107004" },
    { "SGS85", "EPSG:7054" },
    { "andrae", "PROJ:ANDRAE" },
    { "clrk80", "EPSG:7034" },
    { "CPM", "PROJ:CPM" },
    { "delmbr", "PROJ:DELMBR" },
    { "Earth2000", "PROJ:EARTH2000" },
    { "kaula", "PROJ:KAULA" },
    { "lerch", "PROJ:LERCH" },
    { "MERIT", "PROJ:MERIT" },
    { "mprts", "PROJ:MPRTS" },
    { "new_intl", "PROJ:NEW_INTL" },
    { "WGS60", "PROJ:WGS60" }
  };
 
  QString ellipsoid = e;
  // ensure ellipsoid database is populated when first called
  static std::once_flag initialized;
  std::call_once( initialized, [] {
    const QgsScopedRuntimeProfile profile( QObject::tr( "Initialize ellipsoids" ) );
    ( void ) definitions();
  } );
 
  ellipsoid = sProj6EllipsoidAcronymMap.value( ellipsoid, ellipsoid ); // silently upgrade older QGIS acronyms to proj acronyms
 
  // check cache
  {
    const QgsReadWriteLocker locker( *sEllipsoidCacheLock(), QgsReadWriteLocker::Read );
    if ( !sDisableCache )
    {
      const QHash< QString, EllipsoidParameters >::const_iterator cacheIt = sEllipsoidCache()->constFind( ellipsoid );
      if ( cacheIt != sEllipsoidCache()->constEnd() )
      {
        // found a match in the cache
        QgsEllipsoidUtils::EllipsoidParameters params = cacheIt.value();
        return params;
      }
    }
  }
 
  EllipsoidParameters params;
 
  // Check if we have a custom projection, and set from text string.
  // Format is "PARAMETER:<semi-major axis>:<semi minor axis>
  // Numbers must be with (optional) decimal point and no other separators (C locale)
  // Distances in meters.  Flattening is calculated.
  if ( ellipsoid.startsWith( "PARAMETER"_L1 ) )
  {
    QStringList paramList = ellipsoid.split( ':' );
    bool semiMajorOk, semiMinorOk;
    const double semiMajor = paramList[1].toDouble( &semiMajorOk );
    const double semiMinor = paramList[2].toDouble( &semiMinorOk );
    if ( semiMajorOk && semiMinorOk )
    {
      params.semiMajor = semiMajor;
      params.semiMinor = semiMinor;
      params.inverseFlattening = semiMajor / ( semiMajor - semiMinor );
      params.useCustomParameters = true;
      params.crs.createFromProj( u"+proj=longlat +a=%1 +b=%2 +no_defs +type=crs"_s.arg( params.semiMajor, 0, 'g', 17 ).arg( params.semiMinor, 0, 'g', 17 ) );
    }
    else
    {
      params.valid = false;
    }
 
    const QgsReadWriteLocker locker( *sEllipsoidCacheLock(), QgsReadWriteLocker::Write );
    if ( !sDisableCache )
    {
      sEllipsoidCache()->insert( ellipsoid, params );
    }
    return params;
  }
  params.valid = false;
 
  const QgsReadWriteLocker l( *sEllipsoidCacheLock(), QgsReadWriteLocker::Write );
  if ( !sDisableCache )
  {
    sEllipsoidCache()->insert( ellipsoid, params );
  }
 
  return params;
}
 
QList<QgsEllipsoidUtils::EllipsoidDefinition> QgsEllipsoidUtils::definitions()
{
  QgsReadWriteLocker defLocker( *sDefinitionCacheLock(), QgsReadWriteLocker::Read );
  if ( !sDefinitionCache()->isEmpty() )
  {
    return *sDefinitionCache();
  }
  defLocker.changeMode( QgsReadWriteLocker::Write );
 
  QList<QgsEllipsoidUtils::EllipsoidDefinition> defs;
 
  QgsReadWriteLocker locker( *sEllipsoidCacheLock(), QgsReadWriteLocker::Write );
 
  PJ_CONTEXT *context = QgsProjContext::get();
  if ( PROJ_STRING_LIST authorities = proj_get_authorities_from_database( context ) )
  {
    PROJ_STRING_LIST authoritiesIt = authorities;
    while ( char *authority = *authoritiesIt )
    {
      if ( PROJ_STRING_LIST codes = proj_get_codes_from_database( context, authority, PJ_TYPE_ELLIPSOID, 0 ) )
      {
        PROJ_STRING_LIST codesIt = codes;
        while ( char *code = *codesIt )
        {
          const QgsProjUtils::proj_pj_unique_ptr ellipsoid( proj_create_from_database( context, authority, code, PJ_CATEGORY_ELLIPSOID, 0, nullptr ) );
          if ( ellipsoid )
          {
            EllipsoidDefinition def;
            QString name = QString( proj_get_name( ellipsoid.get() ) );
            def.acronym = u"%1:%2"_s.arg( authority, code );
            name.replace( '_', ' ' );
            def.description = u"%1 (%2:%3)"_s.arg( name, authority, code );
 
            def.celestialBodyName = proj_get_celestial_body_name( context, ellipsoid.get() );
 
            double semiMajor, semiMinor, invFlattening;
            int semiMinorComputed = 0;
            if ( proj_ellipsoid_get_parameters( context, ellipsoid.get(), &semiMajor, &semiMinor, &semiMinorComputed, &invFlattening ) )
            {
              def.parameters.semiMajor = semiMajor;
              def.parameters.semiMinor = semiMinor;
              def.parameters.inverseFlattening = invFlattening;
              if ( !semiMinorComputed )
                def.parameters.crs.createFromProj( u"+proj=longlat +a=%1 +b=%2 +no_defs +type=crs"_s.arg( def.parameters.semiMajor, 0, 'g', 17 ).arg( def.parameters.semiMinor, 0, 'g', 17 ), false );
              else if ( !qgsDoubleNear( def.parameters.inverseFlattening, 0.0 ) )
                def.parameters.crs.createFromProj( u"+proj=longlat +a=%1 +rf=%2 +no_defs +type=crs"_s.arg( def.parameters.semiMajor, 0, 'g', 17 ).arg( def.parameters.inverseFlattening, 0, 'g', 17 ), false );
              else
                def.parameters.crs.createFromProj( u"+proj=longlat +a=%1 +no_defs +type=crs"_s.arg( def.parameters.semiMajor, 0, 'g', 17 ), false );
            }
            else
            {
              def.parameters.valid = false;
            }
 
            defs << def;
            if ( !sDisableCache )
            {
              sEllipsoidCache()->insert( def.acronym, def.parameters );
            }
          }
 
          codesIt++;
        }
        proj_string_list_destroy( codes );
      }
 
      authoritiesIt++;
    }
    proj_string_list_destroy( authorities );
  }
  locker.unlock();
 
  QCollator collator;
  collator.setCaseSensitivity( Qt::CaseInsensitive );
  std::sort( defs.begin(), defs.end(), [&collator]( const EllipsoidDefinition &a, const EllipsoidDefinition &b ) { return collator.compare( a.description, b.description ) < 0; } );
  if ( !sDisableCache )
  {
    *sDefinitionCache() = defs;
  }
 
  return defs;
}
 
QStringList QgsEllipsoidUtils::acronyms()
{
  QStringList result;
  const QList<QgsEllipsoidUtils::EllipsoidDefinition> defs = definitions();
  result.reserve( defs.size() );
  for ( const QgsEllipsoidUtils::EllipsoidDefinition &def : defs )
  {
    result << def.acronym;
  }
  return result;
}
 
QList<QgsCelestialBody> QgsEllipsoidUtils::celestialBodies()
{
  return QgsApplication::coordinateReferenceSystemRegistry()->celestialBodies();
}
 
void QgsEllipsoidUtils::invalidateCache( bool disableCache )
{
  const QgsReadWriteLocker locker1( *sEllipsoidCacheLock(), QgsReadWriteLocker::Write );
  const QgsReadWriteLocker locker2( *sDefinitionCacheLock(), QgsReadWriteLocker::Write );
 
  if ( !sDisableCache )
  {
    if ( disableCache )
      sDisableCache = true;
    sEllipsoidCache()->clear();
    sDefinitionCache()->clear();
  }
}
 
QQuaternion QgsEllipsoidUtils::quaternionFromNormalUpRight( const QVector3D &normalUp, const QVector3D &normalRight )
{
  const QVector3D right = normalRight.normalized();
  const QVector3D up = normalUp.normalized();
  // In a right-handed coordinate system with X=right, Z=up:
  // Y = cross(up, right) = forward (pointing North for ENU)
  const QVector3D forward = QVector3D::crossProduct( up, right ).normalized();
 
  // Build rotation matrix with columns [right, forward, up]
  // QGenericMatrix constructor takes row-major input
  const float matData[9] = { right.x(), forward.x(), up.x(), right.y(), forward.y(), up.y(), right.z(), forward.z(), up.z() };
  const QMatrix3x3 rotMatrix( matData );
  return QQuaternion::fromRotationMatrix( rotMatrix );
}
 
QQuaternion QgsEllipsoidUtils::ellipsoidEastNorthUpRotation( const QgsVector3D &position, double semiMajorAxis, double semiMinorAxis )
{
  QgsVector3D up( position.x() / ( semiMajorAxis * semiMajorAxis ), position.y() / ( semiMajorAxis * semiMajorAxis ), position.z() / ( semiMinorAxis * semiMinorAxis ) );
  up.normalize();
 
  QgsVector3D east( -position.y(), position.x(), 0.0 );
  if ( east.length() < 1e-6 )
  {
    east = QgsVector3D( 1.0, 0.0, 0.0 );
  }
  else
  {
    east.normalize();
  }
 
  return quaternionFromNormalUpRight( up.toVector3D(), east.toVector3D() );
}