QGIS API Documentation  3.9.0-Master (224899f119)
qgsellipsoidutils.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsellipsoidutils.cpp
3  ----------------------
4  Date : April 2017
5  Copyright : (C) 2017 by Nyall Dawson
6  email : nyall dot dawson at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include "qgsellipsoidutils.h"
17 #include "qgsapplication.h"
18 #include "qgslogger.h"
19 #include "qgsmessagelog.h"
20 #include <sqlite3.h>
21 #include <QCollator>
22 #include "qgsprojutils.h"
23 #include "qgsreadwritelocker.h"
24 
25 #if PROJ_VERSION_MAJOR>=6
26 #include <proj.h>
27 #include <mutex>
28 #endif
29 
30 QReadWriteLock QgsEllipsoidUtils::sEllipsoidCacheLock;
31 QHash< QString, QgsEllipsoidUtils::EllipsoidParameters > QgsEllipsoidUtils::sEllipsoidCache;
32 QReadWriteLock QgsEllipsoidUtils::sDefinitionCacheLock;
33 QList< QgsEllipsoidUtils::EllipsoidDefinition > QgsEllipsoidUtils::sDefinitionCache;
34 
35 // maps older QGIS ellipsoid acronyms to proj acronyms/names
36 const QMap< QString, QString > sProj6EllipsoidAcronymMap
37 {
38  { "clrk80", "clrk80ign" },
39  {"Adrastea2000", "ESRI:107909"},
40  {"Amalthea2000", "ESRI:107910"},
41  {"Ananke2000", "ESRI:107911"},
42  {"Ariel2000", "ESRI:107945"},
43  {"Atlas2000", "ESRI:107926"},
44  {"Belinda2000", "ESRI:107946"},
45  {"Bianca2000", "ESRI:107947"},
46  {"Callisto2000", "ESRI:107912"},
47  {"Calypso2000", "ESRI:107927"},
48  {"Carme2000", "ESRI:107913"},
49  {"Charon2000", "ESRI:107970"},
50  {"Cordelia2000", "ESRI:107948"},
51  {"Cressida2000", "ESRI:107949"},
52  {"Deimos2000", "ESRI:107906"},
53  {"Desdemona2000", "ESRI:107950"},
54  {"Despina2000", "ESRI:107961"},
55  {"Dione2000", "ESRI:107928"},
56  {"Elara2000", "ESRI:107914"},
57  {"Enceladus2000", "ESRI:107929"},
58  {"Epimetheus2000", "ESRI:107930"},
59  {"Europa2000", "ESRI:107915"},
60  {"Galatea2000", "ESRI:107962"},
61  {"Ganymede2000", "ESRI:107916"},
62  {"Helene2000", "ESRI:107931"},
63  {"Himalia2000", "ESRI:107917"},
64  {"Hyperion2000", "ESRI:107932"},
65  {"Iapetus2000", "ESRI:107933"},
66  {"Io2000", "ESRI:107918"},
67  {"Janus2000", "ESRI:107934"},
68  {"Juliet2000", "ESRI:107951"},
69  {"Jupiter2000", "ESRI:107908"},
70  {"Larissa2000", "ESRI:107963"},
71  {"Leda2000", "ESRI:107919"},
72  {"Lysithea2000", "ESRI:107920"},
73  {"Mars2000", "ESRI:107905"},
74  {"Mercury2000", "ESRI:107900"},
75  {"Metis2000", "ESRI:107921"},
76  {"Mimas2000", "ESRI:107935"},
77  {"Miranda2000", "ESRI:107952"},
78  {"Moon2000", "ESRI:107903"},
79  {"Naiad2000", "ESRI:107964"},
80  {"Neptune2000", "ESRI:107960"},
81  {"Nereid2000", "ESRI:107965"},
82  {"Oberon2000", "ESRI:107953"},
83  {"Ophelia2000", "ESRI:107954"},
84  {"Pan2000", "ESRI:107936"},
85  {"Pandora2000", "ESRI:107937"},
86  {"Pasiphae2000", "ESRI:107922"},
87  {"Phobos2000", "ESRI:107907"},
88  {"Phoebe2000", "ESRI:107938"},
89  {"Pluto2000", "ESRI:107969"},
90  {"Portia2000", "ESRI:107955"},
91  {"Prometheus2000", "ESRI:107939"},
92  {"Proteus2000", "ESRI:107966"},
93  {"Puck2000", "ESRI:107956"},
94  {"Rhea2000", "ESRI:107940"},
95  {"Rosalind2000", "ESRI:107957"},
96  {"Saturn2000", "ESRI:107925"},
97  {"Sinope2000", "ESRI:107923"},
98  {"Telesto2000", "ESRI:107941"},
99  {"Tethys2000", "ESRI:107942"},
100  {"Thalassa2000", "ESRI:107967"},
101  {"Thebe2000", "ESRI:107924"},
102  {"Titan2000", "ESRI:107943"},
103  {"Titania2000", "ESRI:107958"},
104  {"Triton2000", "ESRI:107968"},
105  {"Umbriel2000", "ESRI:107959"},
106  {"Uranus2000", "ESRI:107944"},
107  {"Venus2000", "ESRI:107902"},
108  {"IGNF:ELG053", "EPSG:7030"},
109  {"IGNF:ELG052", "EPSG:7043"},
110  {"IGNF:ELG102", "EPSG:7043"},
111  {"WGS66", "ESRI:107001"},
112  {"plessis", "EPSG:7027"},
113  {"IGNF:ELG017", "EPSG:7027"},
114  {"mod_airy", "EPSG:7002"},
115  {"IGNF:ELG037", "EPSG:7019"},
116  {"IGNF:ELG108", "EPSG:7036"},
117  {"cape", "EPSG:7034"},
118  {"IGNF:ELG010", "EPSG:7011"},
119  {"IGNF:ELG003", "EPSG:7012"},
120  {"IGNF:ELG004", "EPSG:7008"},
121  {"GSK2011", "EPSG:1025"},
122  {"airy", "EPSG:7001"},
123  {"aust_SA", "EPSG:7003"},
124  {"bessel", "EPSG:7004"},
125  {"clrk66", "EPSG:7008"},
126  {"clrk80ign", "EPSG:7011"},
127  {"evrst30", "EPSG:7015"},
128  {"evrstSS", "EPSG:7016"},
129  {"evrst48", "EPSG:7018"},
130  {"GRS80", "EPSG:7019"},
131  {"helmert", "EPSG:7020"},
132  {"intl", "EPSG:7022"},
133  {"krass", "EPSG:7024"},
134  {"NWL9D", "EPSG:7025"},
135  {"WGS84", "EPSG:7030"},
136  {"GRS67", "EPSG:7036"},
137  {"WGS72", "EPSG:7043"},
138  {"bess_nam", "EPSG:7046"},
139  {"IAU76", "EPSG:7049"},
140  {"sphere", "EPSG:7052"},
141  {"hough", "EPSG:7053"},
142  {"evrst69", "EPSG:7056"},
143  {"fschr60", "ESRI:107002"},
144  {"fschr68", "ESRI:107003"},
145  {"fschr60m", "ESRI:107004"},
146  {"walbeck", "ESRI:107007"},
147  {"IGNF:ELG001", "EPSG:7022"},
148  {"engelis", "EPSG:7054"},
149  {"evrst56", "EPSG:7044"},
150  {"SEasia", "ESRI:107004"},
151  {"SGS85", "EPSG:7054"},
152  {"andrae", "PROJ:ANDRAE"},
153  {"clrk80", "EPSG:7034"},
154  {"CPM", "PROJ:CPM"},
155  {"delmbr", "PROJ:DELMBR"},
156  {"Earth2000", "PROJ:EARTH2000"},
157  {"kaula", "PROJ:KAULA"},
158  {"lerch", "PROJ:LERCH"},
159  {"MERIT", "PROJ:MERIT"},
160  {"mprts", "PROJ:MPRTS"},
161  {"new_intl", "PROJ:NEW_INTL"},
162  {"WGS60", "PROJ:WGS60"}
163 };
164 
166 {
167  QString ellipsoid = e;
168 #if PROJ_VERSION_MAJOR >= 6
169  // ensure ellipsoid database is populated when first called
170  static std::once_flag initialized;
171  std::call_once( initialized, [ = ]
172  {
173  ( void )definitions();
174  } );
175 
176  ellipsoid = sProj6EllipsoidAcronymMap.value( ellipsoid, ellipsoid ); // silently upgrade older QGIS acronyms to proj acronyms
177 #endif
178 
179  // check cache
180  {
181  QgsReadWriteLocker locker( sEllipsoidCacheLock, QgsReadWriteLocker::Read );
182  QHash< QString, EllipsoidParameters >::const_iterator cacheIt = sEllipsoidCache.constFind( ellipsoid );
183  if ( cacheIt != sEllipsoidCache.constEnd() )
184  {
185  // found a match in the cache
186  QgsEllipsoidUtils::EllipsoidParameters params = cacheIt.value();
187  return params;
188  }
189  }
190 
191  EllipsoidParameters params;
192 
193  // Check if we have a custom projection, and set from text string.
194  // Format is "PARAMETER:<semi-major axis>:<semi minor axis>
195  // Numbers must be with (optional) decimal point and no other separators (C locale)
196  // Distances in meters. Flattening is calculated.
197  if ( ellipsoid.startsWith( QLatin1String( "PARAMETER" ) ) )
198  {
199  QStringList paramList = ellipsoid.split( ':' );
200  bool semiMajorOk, semiMinorOk;
201  double semiMajor = paramList[1].toDouble( & semiMajorOk );
202  double semiMinor = paramList[2].toDouble( & semiMinorOk );
203  if ( semiMajorOk && semiMinorOk )
204  {
205  params.semiMajor = semiMajor;
206  params.semiMinor = semiMinor;
207  params.inverseFlattening = semiMajor / ( semiMajor - semiMinor );
208  params.useCustomParameters = true;
209  }
210  else
211  {
212  params.valid = false;
213  }
214 
215  QgsReadWriteLocker locker( sEllipsoidCacheLock, QgsReadWriteLocker::Write );
216  sEllipsoidCache.insert( ellipsoid, params );
217  return params;
218  }
219 
220 #if PROJ_VERSION_MAJOR< 6
221  // cache miss - get from database
222  // NOT REQUIRED FOR PROJ >= 6 -- we populate known types once by calling definitions() above
223 
224  QString radius, parameter2;
225  //
226  // SQLITE3 stuff - get parameters for selected ellipsoid
227  //
230  // Continue with PROJ list of ellipsoids.
231 
232  //check the db is available
233  int result = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
234  if ( result )
235  {
236  QgsMessageLog::logMessage( QObject::tr( "Can not open srs database (%1): %2" ).arg( QgsApplication::srsDatabaseFilePath(), database.errorMessage() ) );
237  // XXX This will likely never happen since on open, sqlite creates the
238  // database if it does not exist.
239  return params;
240  }
241  // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
242  QString sql = "select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid + '\'';
243  statement = database.prepare( sql, result );
244  // XXX Need to free memory from the error msg if one is set
245  if ( result == SQLITE_OK )
246  {
247  if ( statement.step() == SQLITE_ROW )
248  {
249  radius = statement.columnAsText( 0 );
250  parameter2 = statement.columnAsText( 1 );
251  }
252  }
253  // row for this ellipsoid wasn't found?
254  if ( radius.isEmpty() || parameter2.isEmpty() )
255  {
256  QgsDebugMsg( QStringLiteral( "setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
257  params.valid = false;
258  sEllipsoidCacheLock.lockForWrite();
259  sEllipsoidCache.insert( ellipsoid, params );
260  sEllipsoidCacheLock.unlock();
261  return params;
262  }
263 
264  // get major semiaxis
265  if ( radius.left( 2 ) == QLatin1String( "a=" ) )
266  params.semiMajor = radius.midRef( 2 ).toDouble();
267  else
268  {
269  QgsDebugMsg( QStringLiteral( "setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
270  params.valid = false;
271  sEllipsoidCacheLock.lockForWrite();
272  sEllipsoidCache.insert( ellipsoid, params );
273  sEllipsoidCacheLock.unlock();
274  return params;
275  }
276 
277  // get second parameter
278  // one of values 'b' or 'f' is in field parameter2
279  // second one must be computed using formula: invf = a/(a-b)
280  if ( parameter2.left( 2 ) == QLatin1String( "b=" ) )
281  {
282  params.semiMinor = parameter2.midRef( 2 ).toDouble();
283  params.inverseFlattening = params.semiMajor / ( params.semiMajor - params.semiMinor );
284  }
285  else if ( parameter2.left( 3 ) == QLatin1String( "rf=" ) )
286  {
287  params.inverseFlattening = parameter2.midRef( 3 ).toDouble();
288  params.semiMinor = params.semiMajor - ( params.semiMajor / params.inverseFlattening );
289  }
290  else
291  {
292  QgsDebugMsg( QStringLiteral( "setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
293  params.valid = false;
294  sEllipsoidCacheLock.lockForWrite();
295  sEllipsoidCache.insert( ellipsoid, params );
296  sEllipsoidCacheLock.unlock();
297  return params;
298  }
299 
300  QgsDebugMsgLevel( QStringLiteral( "setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( params.semiMajor ).arg( params.semiMinor ).arg( params.inverseFlattening ), 4 );
301 
302 
303  // get spatial ref system for ellipsoid
304  QString proj4 = "+proj=longlat +ellps=" + ellipsoid + " +no_defs";
306  //TODO: createFromProj4 used to save to the user database any new CRS
307  // this behavior was changed in order to separate creation and saving.
308  // Not sure if it necessary to save it here, should be checked by someone
309  // familiar with the code (should also give a more descriptive name to the generated CRS)
310  if ( destCRS.srsid() == 0 )
311  {
312  QString name = QStringLiteral( " * %1 (%2)" )
313  .arg( QObject::tr( "Generated CRS", "A CRS automatically generated from layer info get this prefix for description" ),
314  destCRS.toProj4() );
315  destCRS.saveAsUserCrs( name );
316  }
317  //
318 
319  // set transformation from project CRS to ellipsoid coordinates
320  params.crs = destCRS;
321 
322  sEllipsoidCacheLock.lockForWrite();
323  sEllipsoidCache.insert( ellipsoid, params );
324  sEllipsoidCacheLock.unlock();
325  return params;
326 #else
327  params.valid = false;
328 
329  QgsReadWriteLocker l( sEllipsoidCacheLock, QgsReadWriteLocker::Write );
330  sEllipsoidCache.insert( ellipsoid, params );
331 
332  return params;
333 #endif
334 }
335 
336 QList<QgsEllipsoidUtils::EllipsoidDefinition> QgsEllipsoidUtils::definitions()
337 {
338  QgsReadWriteLocker defLocker( sDefinitionCacheLock, QgsReadWriteLocker::Read );
339  if ( !sDefinitionCache.isEmpty() )
340  {
341  QList<QgsEllipsoidUtils::EllipsoidDefinition> defs = sDefinitionCache;
342  return defs;
343  }
345 
346  QList<QgsEllipsoidUtils::EllipsoidDefinition> defs;
347 
348 #if PROJ_VERSION_MAJOR>=6
349  QgsReadWriteLocker locker( sEllipsoidCacheLock, QgsReadWriteLocker::Write );
350 
351  PJ_CONTEXT *context = QgsProjContext::get();
352  if ( PROJ_STRING_LIST authorities = proj_get_authorities_from_database( context ) )
353  {
354  PROJ_STRING_LIST authoritiesIt = authorities;
355  while ( char *authority = *authoritiesIt )
356  {
357  if ( PROJ_STRING_LIST codes = proj_get_codes_from_database( context, authority, PJ_TYPE_ELLIPSOID, 0 ) )
358  {
359  PROJ_STRING_LIST codesIt = codes;
360  while ( char *code = *codesIt )
361  {
362  QgsProjUtils::proj_pj_unique_ptr ellipsoid( proj_create_from_database( context, authority, code, PJ_CATEGORY_ELLIPSOID, 0, nullptr ) );
363  if ( ellipsoid.get() )
364  {
366  QString name = QString( proj_get_name( ellipsoid.get() ) );
367  def.acronym = QStringLiteral( "%1:%2" ).arg( authority, code );
368  name.replace( '_', ' ' );
369  def.description = QStringLiteral( "%1 (%2:%3)" ).arg( name, authority, code );
370 
371  double semiMajor, semiMinor, invFlattening;
372  int semiMinorComputed = 0;
373  if ( proj_ellipsoid_get_parameters( context, ellipsoid.get(), &semiMajor, &semiMinor, &semiMinorComputed, &invFlattening ) )
374  {
375  def.parameters.semiMajor = semiMajor;
376  def.parameters.semiMinor = semiMinor;
377  def.parameters.inverseFlattening = invFlattening;
378  if ( !semiMinorComputed )
379  def.parameters.crs = QgsCoordinateReferenceSystem::fromProj4( QStringLiteral( "+proj=longlat +a=%1 +b=%2 +no_defs +type=crs" ).arg( def.parameters.semiMajor, 0, 'g', 17 ).arg( def.parameters.semiMinor, 0, 'g', 17 ) );
380  else if ( !qgsDoubleNear( def.parameters.inverseFlattening, 0.0 ) )
381  def.parameters.crs = QgsCoordinateReferenceSystem::fromProj4( QStringLiteral( "+proj=longlat +a=%1 +rf=%2 +no_defs +type=crs" ).arg( def.parameters.semiMajor, 0, 'g', 17 ).arg( def.parameters.inverseFlattening, 0, 'g', 17 ) );
382  else
383  def.parameters.crs = QgsCoordinateReferenceSystem::fromProj4( QStringLiteral( "+proj=longlat +a=%1 +no_defs +type=crs" ).arg( def.parameters.semiMajor, 0, 'g', 17 ) );
384  }
385  else
386  {
387  def.parameters.valid = false;
388  }
389 
390  defs << def;
391  sEllipsoidCache.insert( def.acronym, def.parameters );
392  }
393 
394  codesIt++;
395  }
396  proj_string_list_destroy( codes );
397  }
398 
399  authoritiesIt++;
400  }
401  proj_string_list_destroy( authorities );
402  }
403  locker.unlock();
404 
405 #else
408  int result;
409 
410  //check the db is available
411  result = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
412  if ( result )
413  {
414  QgsDebugMsg( QStringLiteral( "Can't open database: %1" ).arg( database.errorMessage() ) );
415  // XXX This will likely never happen since on open, sqlite creates the
416  // database if it does not exist.
417  Q_ASSERT( result == 0 );
418  }
419 
420  // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
421  QString sql = QStringLiteral( "select acronym, name from tbl_ellipsoid order by name" );
422  statement = database.prepare( sql, result );
423 
424  if ( result == SQLITE_OK )
425  {
426  while ( statement.step() == SQLITE_ROW )
427  {
429  def.acronym = statement.columnAsText( 0 );
430  def.description = statement.columnAsText( 1 );
431 
432  // use ellipsoidParameters so that result is cached
434 
435  defs << def;
436  }
437  }
438 
439 #endif
440 
441  QCollator collator;
442  collator.setCaseSensitivity( Qt::CaseInsensitive );
443  std::sort( defs.begin(), defs.end(), [&collator]( const EllipsoidDefinition & a, const EllipsoidDefinition & b )
444  {
445  return collator.compare( a.description, b.description ) < 0;
446  } );
447  sDefinitionCache = defs;
448 
449  return defs;
450 }
451 
453 {
454  QStringList result;
455  const auto constDefinitions = definitions();
456  for ( const QgsEllipsoidUtils::EllipsoidDefinition &def : constDefinitions )
457  {
458  result << def.acronym;
459  }
460  return result;
461 }
static QgsCoordinateReferenceSystem fromProj4(const QString &proj4)
Creates a CRS from a proj4 style formatted string.
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)
Contains definition of an ellipsoid.
Unique pointer for sqlite3 prepared statements, which automatically finalizes the statement when the ...
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:265
QString toProj4() const
Returns a Proj4 string representation of this CRS.
Contains parameters for an ellipsoid.
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
long saveAsUserCrs(const QString &name)
Save the proj4-string as a custom CRS.
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
QString errorMessage() const
Returns the most recent error message encountered by the database.
The QgsReadWriteLocker class is a convenience class that simplifies locking and unlocking QReadWriteL...
bool valid
Whether ellipsoid parameters are valid.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
int step()
Steps to the next record in the statement, returning the sqlite3 result code.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
QString columnAsText(int column) const
Returns the column value from the current statement row as a string.
sqlite3_statement_unique_ptr prepare(const QString &sql, int &resultCode) const
Prepares a sql statement, returning the result.
const QMap< QString, QString > sProj6EllipsoidAcronymMap
QString acronym
authority:code for QGIS builds with proj version 6 or greater, or custom acronym for ellipsoid for ea...
Unique pointer for sqlite3 databases, which automatically closes the database when the pointer goes o...
void unlock()
Unlocks the lock.
int open_v2(const QString &path, int flags, const char *zVfs)
Opens the database at the specified file path.
void PJ_CONTEXT
Definition: qgsprojutils.h:135
This class represents a coordinate reference system (CRS).
double inverseFlattening
Inverse flattening.
void changeMode(Mode mode)
Change the mode of the lock to mode.
static QString srsDatabaseFilePath()
Returns the path to the srs.db file.
QgsEllipsoidUtils::EllipsoidParameters parameters
Ellipsoid parameters.
static QList< QgsEllipsoidUtils::EllipsoidDefinition > definitions()
Returns a list of the definitions for all known ellipsoids from the internal ellipsoid database...
static PJ_CONTEXT * get()
Returns a thread local instance of a proj context, safe for use in the current thread.
static QStringList acronyms()
Returns a list of all known ellipsoid acronyms from the internal ellipsoid database.
long srsid() const
Returns the internal CRS ID, if available.
QString description
Description of ellipsoid.