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