QGIS API Documentation  3.4.15-Madeira (e83d02e274)
qgscoordinatetransform_p.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgscoordinatetransform_p.cpp
3  ----------------------------
4  begin : May 2017
5  copyright : (C) 2017 Nyall Dawson
6  email : nyall dot dawson at gmail dot com
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
19 #include "qgslogger.h"
20 #include "qgsapplication.h"
21 
22 extern "C"
23 {
24 #ifndef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
25 #define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
26 #endif
27 #include <proj_api.h>
28 }
29 #include <sqlite3.h>
30 
31 #include <QStringList>
32 
34 
35 #ifdef USE_THREAD_LOCAL
36 thread_local QgsProjContextStore QgsCoordinateTransformPrivate::mProjContext;
37 #else
38 QThreadStorage< QgsProjContextStore * > QgsCoordinateTransformPrivate::mProjContext;
39 #endif
40 
41 QgsProjContextStore::QgsProjContextStore()
42 {
43  context = pj_ctx_alloc();
44 }
45 
46 QgsProjContextStore::~QgsProjContextStore()
47 {
48  pj_ctx_free( context );
49 }
50 
51 QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate()
52 {
53  setFinder();
54 }
55 
56 QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinateReferenceSystem &source,
57  const QgsCoordinateReferenceSystem &destination,
58  const QgsCoordinateTransformContext &context )
59  : mSourceCRS( source )
60  , mDestCRS( destination )
61 {
62  setFinder();
63  calculateTransforms( context );
64 }
65 
66 QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, int sourceDatumTransform, int destDatumTransform )
67  : mSourceCRS( source )
68  , mDestCRS( destination )
69  , mSourceDatumTransform( sourceDatumTransform )
70  , mDestinationDatumTransform( destDatumTransform )
71 {
72  setFinder();
73 }
74 
75 QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinateTransformPrivate &other )
76  : QSharedData( other )
77  , mIsValid( other.mIsValid )
78  , mShortCircuit( other.mShortCircuit )
79  , mSourceCRS( other.mSourceCRS )
80  , mDestCRS( other.mDestCRS )
81  , mSourceDatumTransform( other.mSourceDatumTransform )
82  , mDestinationDatumTransform( other.mDestinationDatumTransform )
83 {
84  //must reinitialize to setup mSourceProjection and mDestinationProjection
85  initialize();
86 }
87 
88 QgsCoordinateTransformPrivate::~QgsCoordinateTransformPrivate()
89 {
90  // free the proj objects
91  freeProj();
92 }
93 
94 bool QgsCoordinateTransformPrivate::checkValidity()
95 {
96  if ( !mSourceCRS.isValid() || !mDestCRS.isValid() )
97  {
98  invalidate();
99  return false;
100  }
101  return true;
102 }
103 
104 void QgsCoordinateTransformPrivate::invalidate()
105 {
106  mShortCircuit = true;
107  mIsValid = false;
108 }
109 
110 bool QgsCoordinateTransformPrivate::initialize()
111 {
112  invalidate();
113  if ( !mSourceCRS.isValid() )
114  {
115  // Pass through with no projection since we have no idea what the layer
116  // coordinates are and projecting them may not be appropriate
117  QgsDebugMsgLevel( QStringLiteral( "Source CRS is invalid!" ), 4 );
118  return false;
119  }
120 
121  if ( !mDestCRS.isValid() )
122  {
123  //No destination projection is set so we set the default output projection to
124  //be the same as input proj.
125  mDestCRS = mSourceCRS;
126  QgsDebugMsgLevel( QStringLiteral( "Destination CRS is invalid!" ), 4 );
127  return false;
128  }
129 
130  mIsValid = true;
131 
132  int sourceDatumTransform = mSourceDatumTransform;
133  int destDatumTransform = mDestinationDatumTransform;
134  bool useDefaultDatumTransform = ( sourceDatumTransform == - 1 && destDatumTransform == -1 );
135 
136  // init the projections (destination and source)
137  freeProj();
138 
139  mSourceProjString = mSourceCRS.toProj4();
140  if ( !useDefaultDatumTransform )
141  {
142  mSourceProjString = stripDatumTransform( mSourceProjString );
143  }
144  if ( sourceDatumTransform != -1 )
145  {
146  mSourceProjString += ( ' ' + QgsDatumTransform::datumTransformToProj( sourceDatumTransform ) );
147  }
148 
149  mDestProjString = mDestCRS.toProj4();
150  if ( !useDefaultDatumTransform )
151  {
152  mDestProjString = stripDatumTransform( mDestProjString );
153  }
154  if ( destDatumTransform != -1 )
155  {
156  mDestProjString += ( ' ' + QgsDatumTransform::datumTransformToProj( destDatumTransform ) );
157  }
158 
159  if ( !useDefaultDatumTransform )
160  {
161  addNullGridShifts( mSourceProjString, mDestProjString, sourceDatumTransform, destDatumTransform );
162  }
163 
164  // create proj projections for current thread
165  QPair<projPJ, projPJ> res = threadLocalProjData();
166 
167 #ifdef COORDINATE_TRANSFORM_VERBOSE
168  QgsDebugMsg( "From proj : " + mSourceCRS.toProj4() );
169  QgsDebugMsg( "To proj : " + mDestCRS.toProj4() );
170 #endif
171 
172  if ( !res.first || !res.second )
173  {
174  mIsValid = false;
175  }
176 
177 #ifdef COORDINATE_TRANSFORM_VERBOSE
178  if ( mIsValid )
179  {
180  QgsDebugMsg( QStringLiteral( "------------------------------------------------------------" ) );
181  QgsDebugMsg( QStringLiteral( "The OGR Coordinate transformation for this layer was set to" ) );
182  QgsLogger::debug<QgsCoordinateReferenceSystem>( "Input", mSourceCRS, __FILE__, __FUNCTION__, __LINE__ );
183  QgsLogger::debug<QgsCoordinateReferenceSystem>( "Output", mDestCRS, __FILE__, __FUNCTION__, __LINE__ );
184  QgsDebugMsg( QStringLiteral( "------------------------------------------------------------" ) );
185  }
186  else
187  {
188  QgsDebugMsg( QStringLiteral( "------------------------------------------------------------" ) );
189  QgsDebugMsg( QStringLiteral( "The OGR Coordinate transformation FAILED TO INITIALIZE!" ) );
190  QgsDebugMsg( QStringLiteral( "------------------------------------------------------------" ) );
191  }
192 #else
193  if ( !mIsValid )
194  {
195  QgsDebugMsg( QStringLiteral( "Coordinate transformation failed to initialize!" ) );
196  }
197 #endif
198 
199  //XXX todo overload == operator for QgsCoordinateReferenceSystem
200  //at the moment srs.parameters contains the whole proj def...soon it won't...
201  //if (mSourceCRS->toProj4() == mDestCRS->toProj4())
202  if ( mSourceCRS == mDestCRS )
203  {
204  // If the source and destination projection are the same, set the short
205  // circuit flag (no transform takes place)
206  mShortCircuit = true;
207  QgsDebugMsgLevel( QStringLiteral( "Source/Dest CRS equal, shortcircuit is set." ), 3 );
208  }
209  else
210  {
211  // Transform must take place
212  mShortCircuit = false;
213  QgsDebugMsgLevel( QStringLiteral( "Source/Dest CRS not equal, shortcircuit is not set." ), 3 );
214  }
215  return mIsValid;
216 }
217 
218 void QgsCoordinateTransformPrivate::calculateTransforms( const QgsCoordinateTransformContext &context )
219 {
220  // recalculate datum transforms from context
221  QgsDatumTransform::TransformPair transforms = context.calculateDatumTransforms( mSourceCRS, mDestCRS );
222  mSourceDatumTransform = transforms.sourceTransformId;
223  mDestinationDatumTransform = transforms.destinationTransformId;
224 }
225 
226 QPair<projPJ, projPJ> QgsCoordinateTransformPrivate::threadLocalProjData()
227 {
228  mProjLock.lockForRead();
229 
230 #ifdef USE_THREAD_LOCAL
231  QMap < uintptr_t, QPair< projPJ, projPJ > >::const_iterator it = mProjProjections.constFind( reinterpret_cast< uintptr_t>( mProjContext.get() ) );
232 #else
233  projCtx pContext = nullptr;
234  if ( mProjContext.hasLocalData() )
235  {
236  pContext = mProjContext.localData()->get();
237  }
238  else
239  {
240  mProjContext.setLocalData( new QgsProjContextStore() );
241  pContext = mProjContext.localData()->get();
242  }
243  QMap < uintptr_t, QPair< projPJ, projPJ > >::const_iterator it = mProjProjections.constFind( reinterpret_cast< uintptr_t>( pContext ) );
244 #endif
245 
246  if ( it != mProjProjections.constEnd() )
247  {
248  QPair<projPJ, projPJ> res = it.value();
249  mProjLock.unlock();
250  return res;
251  }
252 
253  // proj projections don't exist yet, so we need to create
254  mProjLock.unlock();
255  mProjLock.lockForWrite();
256 
257 #ifdef USE_THREAD_LOCAL
258  QPair<projPJ, projPJ> res = qMakePair( pj_init_plus_ctx( mProjContext.get(), mSourceProjString.toUtf8() ),
259  pj_init_plus_ctx( mProjContext.get(), mDestProjString.toUtf8() ) );
260  mProjProjections.insert( reinterpret_cast< uintptr_t>( mProjContext.get() ), res );
261 #else
262  QPair<projPJ, projPJ> res = qMakePair( pj_init_plus_ctx( pContext, mSourceProjString.toUtf8() ),
263  pj_init_plus_ctx( pContext, mDestProjString.toUtf8() ) );
264  mProjProjections.insert( reinterpret_cast< uintptr_t>( pContext ), res );
265 #endif
266  mProjLock.unlock();
267  return res;
268 }
269 
270 QString QgsCoordinateTransformPrivate::stripDatumTransform( const QString &proj4 ) const
271 {
272  QStringList parameterSplit = proj4.split( '+', QString::SkipEmptyParts );
273  QString currentParameter;
274  QString newProjString;
275 
276  for ( int i = 0; i < parameterSplit.size(); ++i )
277  {
278  currentParameter = parameterSplit.at( i );
279  if ( !currentParameter.startsWith( QLatin1String( "towgs84" ), Qt::CaseInsensitive )
280  && !currentParameter.startsWith( QLatin1String( "nadgrids" ), Qt::CaseInsensitive ) )
281  {
282  newProjString.append( '+' );
283  newProjString.append( currentParameter );
284  newProjString.append( ' ' );
285  }
286  }
287  return newProjString;
288 }
289 
290 void QgsCoordinateTransformPrivate::addNullGridShifts( QString &srcProjString, QString &destProjString,
291  int sourceDatumTransform, int destinationDatumTransform ) const
292 {
293  //if one transformation uses ntv2, the other one needs to be null grid shift
294  if ( destinationDatumTransform == -1 && srcProjString.contains( QLatin1String( "+nadgrids" ) ) ) //add null grid if source transformation is ntv2
295  {
296  destProjString += QLatin1String( " [email protected]" );
297  return;
298  }
299  if ( sourceDatumTransform == -1 && destProjString.contains( QLatin1String( "+nadgrids" ) ) )
300  {
301  srcProjString += QLatin1String( " [email protected]" );
302  return;
303  }
304 
305  //add null shift grid for google mercator
306  //(see e.g. http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator)
307  if ( mSourceCRS.authid().compare( QLatin1String( "EPSG:3857" ), Qt::CaseInsensitive ) == 0 && sourceDatumTransform == -1 )
308  {
309  srcProjString += QLatin1String( " [email protected]" );
310  }
311  if ( mDestCRS.authid().compare( QLatin1String( "EPSG:3857" ), Qt::CaseInsensitive ) == 0 && destinationDatumTransform == -1 )
312  {
313  destProjString += QLatin1String( " [email protected]" );
314  }
315 }
316 
317 void QgsCoordinateTransformPrivate::setFinder()
318 {
319 #if 0
320  // Attention! It should be possible to set PROJ_LIB
321  // but it can happen that it was previously set by installer
322  // (version 0.7) and the old installation was deleted
323 
324  // Another problem: PROJ checks if pj_finder was set before
325  // PROJ_LIB environment variable. pj_finder is probably set in
326  // GRASS gproj library when plugin is loaded, consequently
327  // PROJ_LIB is ignored
328 
329  pj_set_finder( finder );
330 #endif
331 }
332 
333 void QgsCoordinateTransformPrivate::freeProj()
334 {
335  mProjLock.lockForWrite();
336  if ( !mProjProjections.isEmpty() )
337  {
338  projCtx tmpContext = pj_ctx_alloc();
339  QMap < uintptr_t, QPair< projPJ, projPJ > >::const_iterator it = mProjProjections.constBegin();
340  for ( ; it != mProjProjections.constEnd(); ++it )
341  {
342  pj_set_ctx( it.value().first, tmpContext );
343  pj_free( it.value().first );
344  pj_set_ctx( it.value().second, tmpContext );
345  pj_free( it.value().second );
346  }
347  pj_ctx_free( tmpContext );
348  mProjProjections.clear();
349  }
350  mProjLock.unlock();
351 }
352 
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
QgsDatumTransform::TransformPair calculateDatumTransforms(const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination) const
Returns the pair of source and destination datum transforms to use for a transform from the specified...
Contains information about the context in which a coordinate transform is executed.
int destinationTransformId
ID for the datum transform to use when projecting to the destination CRS.
Contains datum transform information.
This class represents a coordinate reference system (CRS).
static QString datumTransformToProj(int datumTransformId)
Returns a proj string representing the specified datumTransformId datum transform ID...
int sourceTransformId
ID for the datum transform to use when projecting from the source CRS.
const char * finder(const char *name)