Skip to content

Commit

Permalink
When the selected coordinate operation fails to project coordinates,
Browse files Browse the repository at this point in the history
try falling back to a default PROJ-provided operation to transform the
coordinates

This allows **some** transformation of coordinates to occur in the situation:
- user has selected a coordinate operation utilising a grid shift file
- one or more points which fall OUTSIDE the grid file are transformed

Not for merge -- likely needs some form of user-visible-warning when
this occurs

Refs qgis#33929
  • Loading branch information
nyalldawson committed Feb 17, 2020
1 parent 306262d commit ec78431
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 1 deletion.
52 changes: 52 additions & 0 deletions src/core/qgscoordinatetransform.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -637,6 +637,13 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
return;
}

double *xprev = new double[ numPoints ];
memcpy( xprev, x, sizeof( double ) * numPoints );
double *yprev = new double[ numPoints ];
memcpy( yprev, y, sizeof( double ) * numPoints );
double *zprev = new double[ numPoints ];
memcpy( zprev, z, sizeof( double ) * numPoints );

#ifdef COORDINATE_TRANSFORM_VERBOSE
double xorg = *x;
double yorg = *y;
Expand Down Expand Up @@ -670,9 +677,15 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
// bit more subtle than that, and I'm not completely sure the logic in
// pj_transform() was really sane & fully bullet proof
// So here just check proj_errno() for single point transform
int actualRes = 0;
if ( numPoints == 1 )
{
projResult = proj_errno( projData );
actualRes = projResult;
}
else
{
actualRes = proj_errno( projData );
}
#else
bool sourceIsLatLong = false;
Expand Down Expand Up @@ -707,6 +720,45 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
}
#endif

#if PROJ_VERSION_MAJOR>=6

if ( actualRes != 0 )
{
// fail #1 -- try with getting proj to auto-pick an appropriate coordinate operation for the points
if ( PJ *transform = d->threadLocalFallbackProjData() )
{
projResult = 0;
proj_errno_reset( transform );
proj_trans_generic( transform, direction == ForwardTransform ? PJ_FWD : PJ_INV,
xprev, sizeof( double ), numPoints,
yprev, sizeof( double ), numPoints,
zprev, sizeof( double ), numPoints,
nullptr, sizeof( double ), 0 );
// Try to - approximatively - emulate the behavior of pj_transform()...
// In the case of a single point transform, and a transformation error occurs,
// pj_transform() would return the errno. In cases of multiple point transform,
// it would continue (for non-transient errors, that is pipeline definition
// errors) and just set the resulting x,y to infinity. This is in fact a
// bit more subtle than that, and I'm not completely sure the logic in
// pj_transform() was really sane & fully bullet proof
// So here just check proj_errno() for single point transform
if ( numPoints == 1 )
{
projResult = proj_errno( transform );
}

if ( projResult == 0 )
{
memcpy( x, xprev, sizeof( double ) * numPoints );
memcpy( y, yprev, sizeof( double ) * numPoints );
memcpy( z, zprev, sizeof( double ) * numPoints );
}
}
}
delete [] xprev;
delete [] yprev;
delete [] zprev;
#endif
if ( projResult != 0 )
{
//something bad happened....
Expand Down
45 changes: 44 additions & 1 deletion src/core/qgscoordinatetransform_p.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -585,6 +585,33 @@ ProjData QgsCoordinateTransformPrivate::threadLocalProjData()
return res;
}

#if PROJ_VERSION_MAJOR>=6
ProjData QgsCoordinateTransformPrivate::threadLocalFallbackProjData()
{
QgsReadWriteLocker locker( mProjLock, QgsReadWriteLocker::Read );

PJ_CONTEXT *context = QgsProjContext::get();
QMap < uintptr_t, ProjData >::const_iterator it = mProjFallbackProjections.constFind( reinterpret_cast< uintptr_t>( context ) );

if ( it != mProjFallbackProjections.constEnd() )
{
ProjData res = it.value();
return res;
}

// proj projections don't exist yet, so we need to create
locker.changeMode( QgsReadWriteLocker::Write );

QgsProjUtils::proj_pj_unique_ptr transform( proj_create_crs_to_crs_from_pj( context, mSourceCRS.projObject(), mDestCRS.projObject(), nullptr, nullptr ) );
if ( transform )
transform.reset( proj_normalize_for_visualization( QgsProjContext::get(), transform.get() ) );

ProjData res = transform.release();
mProjFallbackProjections.insert( reinterpret_cast< uintptr_t>( context ), res );
return res;
}
#endif

void QgsCoordinateTransformPrivate::setCustomMissingRequiredGridHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QgsDatumTransform::GridDetails & )> &handler )
{
sMissingRequiredGridHandler = handler;
Expand Down Expand Up @@ -657,7 +684,7 @@ void QgsCoordinateTransformPrivate::addNullGridShifts( QString &srcProjString, Q
void QgsCoordinateTransformPrivate::freeProj()
{
QgsReadWriteLocker locker( mProjLock, QgsReadWriteLocker::Write );
if ( mProjProjections.isEmpty() )
if ( mProjProjections.isEmpty() && mProjFallbackProjections.isEmpty() )
return;
QMap < uintptr_t, ProjData >::const_iterator it = mProjProjections.constBegin();
#if PROJ_VERSION_MAJOR>=6
Expand All @@ -672,6 +699,14 @@ void QgsCoordinateTransformPrivate::freeProj()
proj_assign_context( it.value(), tmpContext );
proj_destroy( it.value() );
}

it = mProjFallbackProjections.constBegin();
for ( ; it != mProjFallbackProjections.constEnd(); ++it )
{
proj_assign_context( it.value(), tmpContext );
proj_destroy( it.value() );
}

proj_context_destroy( tmpContext );
#else
projCtx tmpContext = pj_ctx_alloc();
Expand All @@ -685,6 +720,7 @@ void QgsCoordinateTransformPrivate::freeProj()
pj_ctx_free( tmpContext );
#endif
mProjProjections.clear();
mProjFallbackProjections.clear();
}

#if PROJ_VERSION_MAJOR>=6
Expand All @@ -699,6 +735,13 @@ bool QgsCoordinateTransformPrivate::removeObjectsBelongingToCurrentThread( void
mProjProjections.erase( it );
}

it = mProjFallbackProjections.find( reinterpret_cast< uintptr_t>( pj_context ) );
if ( it != mProjFallbackProjections.end() )
{
proj_destroy( it.value() );
mProjFallbackProjections.erase( it );
}

return mProjProjections.isEmpty();
}
#endif
Expand Down
3 changes: 3 additions & 0 deletions src/core/qgscoordinatetransform_p.h
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,8 @@ class QgsCoordinateTransformPrivate : public QSharedData
ProjData threadLocalProjData();

#if PROJ_VERSION_MAJOR>=6
ProjData threadLocalFallbackProjData();

// Only meant to be called by QgsCoordinateTransform::removeFromCacheObjectsBelongingToCurrentThread()
bool removeObjectsBelongingToCurrentThread( void *pj_context );
#endif
Expand Down Expand Up @@ -152,6 +154,7 @@ class QgsCoordinateTransformPrivate : public QSharedData

QReadWriteLock mProjLock;
QMap < uintptr_t, ProjData > mProjProjections;
QMap < uintptr_t, ProjData > mProjFallbackProjections;

/**
* Sets a custom handler to use when a coordinate transform is created between \a sourceCrs and
Expand Down

0 comments on commit ec78431

Please sign in to comment.