qgsprojutils.cpp

Go to the documentation of this file.

/***************************************************************************
                             qgsprojutils.h
                             -------------------
    begin                : March 2019
    copyright            : (C) 2019 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 "qgsprojutils.h"
 
#include <proj.h>
#include <proj_experimental.h>
 
#include "qgis.h"
#include "qgscoordinatetransform.h"
#include "qgsexception.h"
#include "qgslogger.h"
 
#include <QDate>
#include <QRegularExpression>
#include <QSet>
#include <QString>
 
using namespace Qt::StringLiterals;
 
#if defined(USE_THREAD_LOCAL) && !defined(Q_OS_WIN)
thread_local QgsProjContext QgsProjContext::sProjContext;
#else
QThreadStorage< QgsProjContext * > QgsProjContext::sProjContext;
#endif
 
QgsProjContext::QgsProjContext()
{
  mContext = proj_context_create();
}
 
QgsProjContext::~QgsProjContext()
{
  // Call removeFromCacheObjectsBelongingToCurrentThread() before
  // destroying the context
  QgsCoordinateTransform::removeFromCacheObjectsBelongingToCurrentThread( mContext );
  QgsCoordinateReferenceSystem::removeFromCacheObjectsBelongingToCurrentThread( mContext );
  proj_context_destroy( mContext );
}
 
PJ_CONTEXT *QgsProjContext::get()
{
#if defined(USE_THREAD_LOCAL) && !defined(Q_OS_WIN)
  return sProjContext.mContext;
#else
  PJ_CONTEXT *pContext = nullptr;
  if ( sProjContext.hasLocalData() )
  {
    pContext = sProjContext.localData()->mContext;
  }
  else
  {
    sProjContext.setLocalData( new QgsProjContext() );
    pContext = sProjContext.localData()->mContext;
  }
  return pContext;
#endif
}
 
void QgsProjUtils::ProjPJDeleter::operator()( PJ *object ) const
{
  proj_destroy( object );
}
 
bool QgsProjUtils::usesAngularUnit( const QString &projDef )
{
  const QString crsDef = u"%1 +type=crs"_s.arg( projDef );
  PJ_CONTEXT *context = QgsProjContext::get();
  QgsProjUtils::proj_pj_unique_ptr projSingleOperation( proj_create( context, crsDef.toUtf8().constData() ) );
  if ( !projSingleOperation )
    return false;
 
  QgsProjUtils::proj_pj_unique_ptr coordinateSystem( proj_crs_get_coordinate_system( context, projSingleOperation.get() ) );
  if ( !coordinateSystem )
    return false;
 
  const int axisCount = proj_cs_get_axis_count( context, coordinateSystem.get() );
  if ( axisCount > 0 )
  {
    const char *outUnitAuthName = nullptr;
    const char *outUnitAuthCode = nullptr;
    // Read only first axis
    proj_cs_get_axis_info( context, coordinateSystem.get(), 0,
                           nullptr,
                           nullptr,
                           nullptr,
                           nullptr,
                           nullptr,
                           &outUnitAuthName,
                           &outUnitAuthCode );
 
    if ( outUnitAuthName && outUnitAuthCode )
    {
      const char *unitCategory = nullptr;
      if ( proj_uom_get_info_from_database( context, outUnitAuthName, outUnitAuthCode, nullptr, nullptr, &unitCategory ) )
      {
        return QString( unitCategory ).compare( "angular"_L1, Qt::CaseInsensitive ) == 0;
      }
    }
  }
  return false;
}
 
bool QgsProjUtils::axisOrderIsSwapped( const PJ *crs )
{
  //ported from https://github.com/pramsey/postgis/blob/7ecf6839c57a838e2c8540001a3cd35b78a730db/liblwgeom/lwgeom_transform.c#L299
  //and GDAL OGRSpatialReference::isNorthEastAxisOrder https://github.com/OSGeo/gdal/blob/release/3.10/ogr/ogrspatialreference.cpp#L419
  if ( !crs )
    return false;
 
  PJ_CONTEXT *context = QgsProjContext::get();
 
  QgsProjUtils::proj_pj_unique_ptr horizCrs = crsToHorizontalCrs( crs );
  if ( !horizCrs )
    return false;
 
  QgsProjUtils::proj_pj_unique_ptr pjCs( proj_crs_get_coordinate_system( context, horizCrs.get() ) );
  if ( !pjCs )
    return false;
 
  const int axisCount = proj_cs_get_axis_count( context, pjCs.get() );
  if ( axisCount > 0 )
  {
    const char *outDirection0 = nullptr;
    const char *outDirection1 = nullptr;
    const char *outName0 = nullptr;
    const char *outName1 = nullptr;
 
    proj_cs_get_axis_info( context, pjCs.get(), 0,
                           &outName0,
                           nullptr,
                           &outDirection0,
                           nullptr,
                           nullptr,
                           nullptr,
                           nullptr
                         );
 
    proj_cs_get_axis_info( context, pjCs.get(), 1,
                           &outName1,
                           nullptr,
                           &outDirection1,
                           nullptr,
                           nullptr,
                           nullptr,
                           nullptr
                         );
 
    if ( QString( outDirection0 ).compare( "north"_L1, Qt::CaseInsensitive ) == 0 &&
         QString( outDirection1 ).compare( "east"_L1, Qt::CaseInsensitive ) == 0 )
    {
      return true;
    }
 
    // Handle polar projections with NE-order
    if ( ( QString( outDirection0 ).compare( "north"_L1, Qt::CaseInsensitive ) == 0 &&
           QString( outDirection1 ).compare( "north"_L1, Qt::CaseInsensitive ) == 0 ) ||
         ( QString( outDirection0 ).compare( "south"_L1, Qt::CaseInsensitive ) == 0 &&
           QString( outDirection1 ).compare( "south"_L1, Qt::CaseInsensitive ) == 0 ) )
    {
      return QString( outName0 ).startsWith( "northing"_L1, Qt::CaseInsensitive ) &&
             QString( outName1 ).startsWith( "easting"_L1, Qt::CaseInsensitive ) ;
    }
  }
  return false;
}
 
bool QgsProjUtils::isDynamic( const PJ *crs )
{
  // ported from GDAL OGRSpatialReference::IsDynamic()
  bool isDynamic = false;
  PJ_CONTEXT *context = QgsProjContext::get();
 
  // prefer horizontal crs if possible
  proj_pj_unique_ptr candidate = crsToHorizontalCrs( crs );
  if ( !crs )
    candidate = unboundCrs( crs );
 
  proj_pj_unique_ptr datum( candidate ? proj_crs_get_datum( context, candidate.get() ) : nullptr );
  if ( datum )
  {
    const PJ_TYPE type = proj_get_type( datum.get() );
    isDynamic = type == PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME ||
                type == PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME;
    if ( !isDynamic )
    {
      const QString authName( proj_get_id_auth_name( datum.get(), 0 ) );
      const QString code( proj_get_id_code( datum.get(), 0 ) );
      if ( authName == "EPSG"_L1 && code == "6326"_L1 )
      {
        isDynamic = true;
      }
    }
  }
  else
  {
    proj_pj_unique_ptr ensemble( candidate ? proj_crs_get_datum_ensemble( context, candidate.get() ) : nullptr );
    if ( ensemble )
    {
      proj_pj_unique_ptr member( proj_datum_ensemble_get_member( context, ensemble.get(), 0 ) );
      if ( member )
      {
        const PJ_TYPE type = proj_get_type( member.get() );
        isDynamic = type == PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME ||
                    type == PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME;
      }
    }
  }
  return isDynamic;
}
 
QgsProjUtils::proj_pj_unique_ptr QgsProjUtils::crsToHorizontalCrs( const PJ *crs )
{
  if ( !crs )
    return nullptr;
 
  PJ_CONTEXT *context = QgsProjContext::get();
  switch ( proj_get_type( crs ) )
  {
    case PJ_TYPE_COMPOUND_CRS:
    {
      int i = 0;
      QgsProjUtils::proj_pj_unique_ptr res( proj_crs_get_sub_crs( context, crs, i ) );
      while ( res && ( proj_get_type( res.get() ) == PJ_TYPE_VERTICAL_CRS || proj_get_type( res.get() ) == PJ_TYPE_TEMPORAL_CRS ) )
      {
        i++;
        res.reset( proj_crs_get_sub_crs( context, crs, i ) );
      }
      return res;
    }
 
    case PJ_TYPE_VERTICAL_CRS:
      return nullptr;
 
    // maybe other types to handle??
 
    default:
      return unboundCrs( crs );
  }
 
#ifndef _MSC_VER  // unreachable
  return nullptr;
#endif
}
 
QgsProjUtils::proj_pj_unique_ptr QgsProjUtils::crsToVerticalCrs( const PJ *crs )
{
  if ( !crs )
    return nullptr;
 
  PJ_CONTEXT *context = QgsProjContext::get();
  switch ( proj_get_type( crs ) )
  {
    case PJ_TYPE_COMPOUND_CRS:
    {
      int i = 0;
      QgsProjUtils::proj_pj_unique_ptr res( proj_crs_get_sub_crs( context, crs, i ) );
      while ( res && ( proj_get_type( res.get() ) != PJ_TYPE_VERTICAL_CRS ) )
      {
        i++;
        res.reset( proj_crs_get_sub_crs( context, crs, i ) );
      }
      return res;
    }
 
    case PJ_TYPE_VERTICAL_CRS:
      return QgsProjUtils::proj_pj_unique_ptr( proj_clone( context, crs ) );
 
    // maybe other types to handle??
 
    default:
      return nullptr;
  }
 
  BUILTIN_UNREACHABLE
}
 
bool QgsProjUtils::hasVerticalAxis( const PJ *crs )
{
  if ( !crs )
    return false;
 
  PJ_CONTEXT *context = QgsProjContext::get();
 
  switch ( proj_get_type( crs ) )
  {
    case PJ_TYPE_COMPOUND_CRS:
    {
      int i = 0;
      QgsProjUtils::proj_pj_unique_ptr res( proj_crs_get_sub_crs( context, crs, i ) );
      while ( res )
      {
        if ( hasVerticalAxis( res.get() ) )
          return true;
        i++;
        res.reset( proj_crs_get_sub_crs( context, crs, i ) );
      }
      return false;
    }
 
    case PJ_TYPE_BOUND_CRS:
    {
      return hasVerticalAxis( proj_get_source_crs( context, crs ) );
    }
 
    // maybe other types to handle like this??
 
    default:
      break;
  }
 
  QgsProjUtils::proj_pj_unique_ptr pjCs( proj_crs_get_coordinate_system( context, crs ) );
  if ( !pjCs )
    return false;
 
  const int axisCount = proj_cs_get_axis_count( context, pjCs.get() );
  for ( int axisIndex = 0; axisIndex < axisCount; ++axisIndex )
  {
    const char *outDirection = nullptr;
    proj_cs_get_axis_info( context, pjCs.get(), axisIndex,
                           nullptr,
                           nullptr,
                           &outDirection,
                           nullptr,
                           nullptr,
                           nullptr,
                           nullptr
                         );
    const QString outDirectionString = QString( outDirection );
    if ( outDirectionString.compare( "geocentricZ"_L1, Qt::CaseInsensitive ) == 0
         || outDirectionString.compare( "up"_L1, Qt::CaseInsensitive ) == 0
         || outDirectionString.compare( "down"_L1, Qt::CaseInsensitive ) == 0 )
    {
      return true;
    }
  }
  return false;
}
 
QgsProjUtils::proj_pj_unique_ptr QgsProjUtils::unboundCrs( const PJ *crs )
{
  if ( !crs )
    return nullptr;
 
  PJ_CONTEXT *context = QgsProjContext::get();
  switch ( proj_get_type( crs ) )
  {
    case PJ_TYPE_BOUND_CRS:
      return QgsProjUtils::proj_pj_unique_ptr( proj_get_source_crs( context, crs ) );
 
    // maybe other types to handle??
 
    default:
      return QgsProjUtils::proj_pj_unique_ptr( proj_clone( context, crs ) );
  }
 
#ifndef _MSC_VER  // unreachable
  return nullptr;
#endif
}
 
QgsProjUtils::proj_pj_unique_ptr QgsProjUtils::crsToDatumEnsemble( const PJ *crs )
{
  if ( !crs )
    return nullptr;
 
  PJ_CONTEXT *context = QgsProjContext::get();
  QgsProjUtils::proj_pj_unique_ptr candidate = crsToHorizontalCrs( crs );
  if ( !candidate ) // purely vertical CRS
    candidate = unboundCrs( crs );
 
  if ( !candidate )
    return nullptr;
 
  return QgsProjUtils::proj_pj_unique_ptr( proj_crs_get_datum_ensemble( context, candidate.get() ) );
}
 
void QgsProjUtils::proj_collecting_logger( void *user_data, int /*level*/, const char *message )
{
  QStringList *dest = reinterpret_cast< QStringList * >( user_data );
  QString messageString( message );
  messageString.replace( "internal_proj_create: "_L1, QString() );
  dest->append( messageString );
}
 
void QgsProjUtils::proj_silent_logger( void * /*user_data*/, int /*level*/, const char * /*message*/ )
{
}
 
void QgsProjUtils::proj_logger( void *, int level, const char *message )
{
#ifdef QGISDEBUG
  if ( level == PJ_LOG_ERROR )
  {
    const QString messageString( message );
    if ( messageString == "push: Invalid latitude"_L1 )
    {
      // these messages tend to spam the console as they can be repeated 1000s of times
      QgsDebugMsgLevel( messageString, 3 );
    }
    else
    {
      QgsDebugError( messageString );
    }
  }
  else if ( level == PJ_LOG_DEBUG )
  {
    QgsDebugMsgLevel( QString( message ), 3 );
  }
#else
  ( void )level;
  ( void )message;
#endif
}
 
QgsProjUtils::proj_pj_unique_ptr QgsProjUtils::createCompoundCrs( const PJ *horizontalCrs, const PJ *verticalCrs, QStringList *errors )
{
  if ( !horizontalCrs || !verticalCrs )
    return nullptr;
 
  PJ_CONTEXT *context = QgsProjContext::get();
  // collect errors instead of dumping them to terminal
 
  QgsScopedProjCollectingLogger projLogger;
 
  // const cast here is for compatibility with proj < 9.5
  QgsProjUtils::proj_pj_unique_ptr compoundCrs( proj_create_compound_crs( context,
      nullptr,
      const_cast< PJ *>( horizontalCrs ),
      const_cast< PJ * >( verticalCrs ) ) );
 
  if ( errors )
    *errors = projLogger.errors();
 
  return compoundCrs;
}
 
bool QgsProjUtils::identifyCrs( const PJ *crs, QString &authName, QString &authCode, IdentifyFlags flags )
{
  authName.clear();
  authCode.clear();
 
  if ( !crs )
    return false;
 
  int *confidence = nullptr;
  if ( PJ_OBJ_LIST *crsList = proj_identify( QgsProjContext::get(), crs, nullptr, nullptr, &confidence ) )
  {
    const int count = proj_list_get_count( crsList );
    int bestConfidence = 0;
    QgsProjUtils::proj_pj_unique_ptr matchedCrs;
    for ( int i = 0; i < count; ++i )
    {
      if ( confidence[i] >= bestConfidence )
      {
        QgsProjUtils::proj_pj_unique_ptr candidateCrs( proj_list_get( QgsProjContext::get(), crsList, i ) );
        switch ( proj_get_type( candidateCrs.get() ) )
        {
          case PJ_TYPE_BOUND_CRS:
            // proj_identify also matches bound CRSes to the source CRS. But they are not the same as the source CRS, so we don't
            // consider them a candidate for a match here (depending on the identify flags, that is!)
            if ( flags & FlagMatchBoundCrsToUnderlyingSourceCrs )
              break;
            else
              continue;
 
          default:
            break;
        }
 
        candidateCrs = QgsProjUtils::unboundCrs( candidateCrs.get() );
        const QString authName( proj_get_id_auth_name( candidateCrs.get(), 0 ) );
        // if a match is identical confidence, we prefer EPSG codes for compatibility with earlier qgis conversions
        if ( confidence[i] > bestConfidence || ( confidence[i] == bestConfidence && authName == "EPSG"_L1 ) )
        {
          bestConfidence = confidence[i];
          matchedCrs = std::move( candidateCrs );
        }
      }
    }
    proj_list_destroy( crsList );
    proj_int_list_destroy( confidence );
    if ( matchedCrs && bestConfidence >= 70 )
    {
      authName = QString( proj_get_id_auth_name( matchedCrs.get(), 0 ) );
      authCode = QString( proj_get_id_code( matchedCrs.get(), 0 ) );
    }
  }
  return !authName.isEmpty() && !authCode.isEmpty();
}
 
bool QgsProjUtils::coordinateOperationIsAvailable( const QString &projDef )
{
  if ( projDef.isEmpty() )
    return true;
 
  PJ_CONTEXT *context = QgsProjContext::get();
  QgsProjUtils::proj_pj_unique_ptr coordinateOperation( proj_create( context, projDef.toUtf8().constData() ) );
  if ( !coordinateOperation )
    return false;
 
  return static_cast< bool >( proj_coordoperation_is_instantiable( context, coordinateOperation.get() ) );
}
 
QList<QgsDatumTransform::GridDetails> QgsProjUtils::gridsUsed( const QString &proj )
{
  const thread_local QRegularExpression regex( u"\\+(?:nad)?grids=(.*?)\\s"_s );
 
  QList< QgsDatumTransform::GridDetails > grids;
  QRegularExpressionMatchIterator matches = regex.globalMatch( proj );
  while ( matches.hasNext() )
  {
    const QRegularExpressionMatch match = matches.next();
    const QString gridName = match.captured( 1 );
    QgsDatumTransform::GridDetails grid;
    grid.shortName = gridName;
    const char *fullName = nullptr;
    const char *packageName = nullptr;
    const char *url = nullptr;
    int directDownload = 0;
    int openLicense = 0;
    int available = 0;
    proj_grid_get_info_from_database( QgsProjContext::get(), gridName.toUtf8().constData(), &fullName, &packageName, &url, &directDownload, &openLicense, &available );
    grid.fullName = QString( fullName );
    grid.packageName = QString( packageName );
    grid.url = QString( url );
    grid.directDownload = directDownload;
    grid.openLicense = openLicense;
    grid.isAvailable = available;
    grids.append( grid );
  }
  return grids;
}
 
#if 0
QStringList QgsProjUtils::nonAvailableGrids( const QString &projDef )
{
  if ( projDef.isEmpty() )
    return QStringList();
 
  PJ_CONTEXT *context = QgsProjContext::get();
  QgsProjUtils::proj_pj_unique_ptr op( proj_create( context, projDef.toUtf8().constData() ) ); < ---- - this always fails if grids are missing
  if ( !op )
      return QStringList();
 
  QStringList res;
  for ( int j = 0; j < proj_coordoperation_get_grid_used_count( context, op.get() ); ++j )
  {
    const char *shortName = nullptr;
    int isAvailable = 0;
    proj_coordoperation_get_grid_used( context, op.get(), j, &shortName, nullptr, nullptr, nullptr, nullptr, nullptr, &isAvailable );
    if ( !isAvailable )
      res << QString( shortName );
  }
  return res;
}
#endif
 
int QgsProjUtils::projVersionMajor()
{
  return PROJ_VERSION_MAJOR;
}
 
int QgsProjUtils::projVersionMinor()
{
  return PROJ_VERSION_MINOR;
}
 
QString QgsProjUtils::epsgRegistryVersion()
{
  PJ_CONTEXT *context = QgsProjContext::get();
  const char *version = proj_context_get_database_metadata( context, "EPSG.VERSION" );
  return QString( version );
}
 
QDate QgsProjUtils::epsgRegistryDate()
{
  PJ_CONTEXT *context = QgsProjContext::get();
  const char *date = proj_context_get_database_metadata( context, "EPSG.DATE" );
  return QDate::fromString( date, Qt::DateFormat::ISODate );
}
 
QString QgsProjUtils::esriDatabaseVersion()
{
  PJ_CONTEXT *context = QgsProjContext::get();
  const char *version = proj_context_get_database_metadata( context, "ESRI.VERSION" );
  return QString( version );
}
 
QDate QgsProjUtils::esriDatabaseDate()
{
  PJ_CONTEXT *context = QgsProjContext::get();
  const char *date = proj_context_get_database_metadata( context, "ESRI.DATE" );
  return QDate::fromString( date, Qt::DateFormat::ISODate );
}
 
QString QgsProjUtils::ignfDatabaseVersion()
{
  PJ_CONTEXT *context = QgsProjContext::get();
  const char *version = proj_context_get_database_metadata( context, "IGNF.VERSION" );
  return QString( version );
}
 
QDate QgsProjUtils::ignfDatabaseDate()
{
  PJ_CONTEXT *context = QgsProjContext::get();
  const char *date = proj_context_get_database_metadata( context, "IGNF.DATE" );
  return QDate::fromString( date, Qt::DateFormat::ISODate );
}
 
QStringList QgsProjUtils::searchPaths()
{
  const QString path( proj_info().searchpath );
  QStringList paths;
#ifdef Q_OS_WIN
  paths = path.split( ';' );
#else
  paths = path.split( ':' );
#endif
 
  QSet<QString> existing;
  // thin out duplicates from paths -- see https://github.com/OSGeo/proj.4/pull/1498
  QStringList res;
  res.reserve( paths.count() );
  for ( const QString &p : std::as_const( paths ) )
  {
    if ( existing.contains( p ) )
      continue;
 
    existing.insert( p );
    res << p;
  }
  return res;
}
 
//
// QgsScopedProjCollectingLogger
//
 
QgsScopedProjCollectingLogger::QgsScopedProjCollectingLogger()
{
  proj_log_func( QgsProjContext::get(), &mProjErrors, QgsProjUtils::proj_collecting_logger );
}
 
QgsScopedProjCollectingLogger::~QgsScopedProjCollectingLogger()
{
  // reset logger back to terminal output
  proj_log_func( QgsProjContext::get(), nullptr, QgsProjUtils::proj_logger );
}
 
//
// QgsScopedProjSilentLogger
//
 
QgsScopedProjSilentLogger::QgsScopedProjSilentLogger()
{
  proj_log_func( QgsProjContext::get(), nullptr, QgsProjUtils::proj_silent_logger );
}
 
QgsScopedProjSilentLogger::~QgsScopedProjSilentLogger()
{
  // reset logger back to terminal output
  proj_log_func( QgsProjContext::get(), nullptr, QgsProjUtils::proj_logger );
}