Skip to content
Permalink
Browse files
Add basic support for coordinate epoch of dynamic (not plate fixed) crs
to QgsCoordinateTransform

QgsCoordinateTransform class can perform time-dependent transformations
between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
however dynamic CRS to dynamic CRS transformations are not currently supported.

Using the same approach as the GDAL changeset in OSGeo/gdal#3810
  • Loading branch information
nyalldawson committed May 13, 2021
1 parent 420a6c1 commit bd26ae1d2923b22447ed210e3c7e37e0043c1591
@@ -850,6 +850,13 @@ CRS (see :py:func:`~QgsCoordinateReferenceSystem.isDynamic`).

:param epoch: Coordinate epoch as decimal year (e.g. 2021.3)

.. warning::

The :py:class:`QgsCoordinateTransform` class can perform time-dependent transformations
between a static and dynamic CRS based on either the source or destination CRS coordinate epoch,
however dynamic CRS to dynamic CRS transformations are not currently supported.


.. seealso:: :py:func:`coordinateEpoch`


@@ -870,8 +877,15 @@ observation, and not to the CRS, however it is often more practical to
bind it to the CRS. The coordinate epoch should be specified for dynamic
CRS (see :py:func:`~QgsCoordinateReferenceSystem.isDynamic`).

.. warning::

The :py:class:`QgsCoordinateTransform` class can perform time-dependent transformations
between a static and dynamic CRS based on either the source or destination CRS coordinate epoch,
however dynamic CRS to dynamic CRS transformations are not currently supported.

:return: Coordinate epoch as decimal year (e.g. 2021.3), or NaN if not set, or relevant.


.. seealso:: :py:func:`setCoordinateEpoch`


@@ -29,6 +29,12 @@ transforms coordinates from the layer's coordinate system to the map canvas.

Since QGIS 3.0 :py:class:`QgsCoordinateReferenceSystem` objects are implicitly shared.

.. warning::

Since QGIS 3.20 The :py:class:`QgsCoordinateTransform` class can perform time-dependent transformations
between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
however dynamic CRS to dynamic CRS transformations are not currently supported.

.. seealso:: :py:class:`QgsDatumTransform`

.. seealso:: :py:class:`QgsCoordinateTransformContext`
@@ -64,6 +70,12 @@ to utilize.
Python scripts should generally use the constructor variant which accepts
a :py:class:`QgsProject` instance instead of this constructor.

.. warning::

Since QGIS 3.20 The QgsCoordinateTransform class can perform time-dependent transformations
between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
however dynamic CRS to dynamic CRS transformations are not currently supported.

.. warning::

Do NOT use an empty/default constructed QgsCoordinateTransformContext()
@@ -73,6 +85,7 @@ a :py:class:`QgsProject` instance instead of this constructor.
based on the current code context, or use the constructor variant which
accepts a :py:class:`QgsProject` argument instead.


.. versionadded:: 3.0
%End

@@ -97,6 +110,12 @@ correctly respected during coordinate transforms. E.g.
transform = QgsCoordinateTransform(QgsCoordinateReferenceSystem("EPSG:3111"),
QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance())

.. warning::

Since QGIS 3.20 The QgsCoordinateTransform class can perform time-dependent transformations
between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
however dynamic CRS to dynamic CRS transformations are not currently supported.

.. versionadded:: 3.0
%End

@@ -783,6 +783,10 @@ class CORE_EXPORT QgsCoordinateReferenceSystem
*
* \param epoch Coordinate epoch as decimal year (e.g. 2021.3)
*
* \warning The QgsCoordinateTransform class can perform time-dependent transformations
* between a static and dynamic CRS based on either the source or destination CRS coordinate epoch,
* however dynamic CRS to dynamic CRS transformations are not currently supported.
*
* \see coordinateEpoch()
*
* \since QGIS 3.20
@@ -802,6 +806,10 @@ class CORE_EXPORT QgsCoordinateReferenceSystem
* bind it to the CRS. The coordinate epoch should be specified for dynamic
* CRS (see isDynamic()).
*
* \warning The QgsCoordinateTransform class can perform time-dependent transformations
* between a static and dynamic CRS based on either the source or destination CRS coordinate epoch,
* however dynamic CRS to dynamic CRS transformations are not currently supported.
*
* \returns Coordinate epoch as decimal year (e.g. 2021.3), or NaN if not set, or relevant.
*
* \see setCoordinateEpoch()
@@ -654,6 +654,9 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
std::vector< double > zprev( numPoints );
memcpy( zprev.data(), z, sizeof( double ) * numPoints );

const bool useTime = !std::isnan( d->mDefaultTime );
std::vector< double > t( useTime ? numPoints : 0, d->mDefaultTime );

#ifdef COORDINATE_TRANSFORM_VERBOSE
double xorg = *x;
double yorg = *y;
@@ -678,7 +681,7 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
x, sizeof( double ), numPoints,
y, sizeof( double ), numPoints,
z, sizeof( double ), numPoints,
nullptr, sizeof( double ), 0 );
useTime ? t.data() : nullptr, sizeof( double ), useTime ? numPoints : 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,
@@ -723,7 +726,7 @@ void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *
xprev.data(), sizeof( double ), numPoints,
yprev.data(), sizeof( double ), numPoints,
zprev.data(), sizeof( double ), numPoints,
nullptr, sizeof( double ), 0 );
useTime ? t.data() : nullptr, sizeof( double ), useTime ? numPoints : 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,
@@ -896,7 +899,16 @@ bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &s
const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( sourceKey, destKey ) );
for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
{
if ( ( *valIt ).coordinateOperation() == coordinateOperationProj && ( *valIt ).allowFallbackTransforms() == allowFallback )
if ( ( *valIt ).coordinateOperation() == coordinateOperationProj
&& ( *valIt ).allowFallbackTransforms() == allowFallback

// careful here, nan != nan, so we need to explicitly handle the case when both crses have nan coordinateEpoch
&& ( std::isnan( ( *valIt ).sourceCrs().coordinateEpoch() ) == std::isnan( src.coordinateEpoch() )
&& ( std::isnan( src.coordinateEpoch() ) || src.coordinateEpoch() == ( *valIt ).sourceCrs().coordinateEpoch() ) )

&& ( std::isnan( ( *valIt ).destinationCrs().coordinateEpoch() ) == std::isnan( dest.coordinateEpoch() )
&& ( std::isnan( dest.coordinateEpoch() ) || dest.coordinateEpoch() == ( *valIt ).destinationCrs().coordinateEpoch() ) )
)
{
// need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
QgsCoordinateTransformContext context = mContext;
@@ -47,6 +47,10 @@ class QgsVector3D;
* transforms coordinates from the layer's coordinate system to the map canvas.
* \note Since QGIS 3.0 QgsCoordinateReferenceSystem objects are implicitly shared.
*
* \warning Since QGIS 3.20 The QgsCoordinateTransform class can perform time-dependent transformations
* between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
* however dynamic CRS to dynamic CRS transformations are not currently supported.
*
* \see QgsDatumTransform
* \see QgsCoordinateTransformContext
*/
@@ -76,6 +80,10 @@ class CORE_EXPORT QgsCoordinateTransform
* Python scripts should generally use the constructor variant which accepts
* a QgsProject instance instead of this constructor.
*
* \warning Since QGIS 3.20 The QgsCoordinateTransform class can perform time-dependent transformations
* between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
* however dynamic CRS to dynamic CRS transformations are not currently supported.
*
* \warning Do NOT use an empty/default constructed QgsCoordinateTransformContext()
* object when creating QgsCoordinateTransform objects. This prevents correct
* datum transform handling and may result in inaccurate transformations. Always
@@ -107,6 +115,9 @@ class CORE_EXPORT QgsCoordinateTransform
* QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance())
* \endcode
*
* \warning Since QGIS 3.20 The QgsCoordinateTransform class can perform time-dependent transformations
* between a static and dynamic CRS based on either the source OR destination CRS coordinate epoch,
* however dynamic CRS to dynamic CRS transformations are not currently supported.
* \since QGIS 3.0
*/
explicit QgsCoordinateTransform( const QgsCoordinateReferenceSystem &source,
@@ -701,7 +712,15 @@ class CORE_EXPORT QgsCoordinateTransform

// cache
static QReadWriteLock sCacheLock;
static QMultiHash< QPair< QString, QString >, QgsCoordinateTransform > sTransforms; //same auth_id pairs might have different datum transformations

/**
* Map of cached transforms. The keys are formed by pairs of strings uniquely identifying the source and
* destination CRS, using the auth:id were available or a full WKT2 definition where an auth:id is not available.
*
* The same auth_id pairs might have multiple transformations, as they can be based on different coordinate
* operations, allowance of ballpark transforms, and the source or destination coordinate epoch.
*/
static QMultiHash< QPair< QString, QString >, QgsCoordinateTransform > sTransforms;
static bool sDisableCache;


@@ -87,6 +87,11 @@ QgsCoordinateTransformPrivate::QgsCoordinateTransformPrivate( const QgsCoordinat
, mProjCoordinateOperation( other.mProjCoordinateOperation )
, mShouldReverseCoordinateOperation( other.mShouldReverseCoordinateOperation )
, mAllowFallbackTransforms( other.mAllowFallbackTransforms )
, mSourceIsDynamic( other.mSourceIsDynamic )
, mDestIsDynamic( other.mDestIsDynamic )
, mSourceCoordinateEpoch( other.mSourceCoordinateEpoch )
, mDestCoordinateEpoch( other.mDestCoordinateEpoch )
, mDefaultTime( other.mDefaultTime )
, mIsReversed( other.mIsReversed )
, mProjLock()
, mProjProjections()
@@ -150,6 +155,28 @@ bool QgsCoordinateTransformPrivate::initialize()
return true;
}

mSourceIsDynamic = mSourceCRS.isDynamic();
mSourceCoordinateEpoch = mSourceCRS.coordinateEpoch();
mDestIsDynamic = mDestCRS.isDynamic();
mDestCoordinateEpoch = mDestCRS.coordinateEpoch();

// Determine the default coordinate epoch.
// For time-dependent transformations, PROJ can currently only do
// staticCRS -> dynamicCRS or dynamicCRS -> staticCRS transformations, and
// in either case, the coordinate epoch of the dynamicCRS must be provided
// as the input time.
mDefaultTime = ( mSourceIsDynamic && !std::isnan( mSourceCoordinateEpoch ) && !mDestIsDynamic )
? mSourceCoordinateEpoch
: ( mDestIsDynamic && !std::isnan( mDestCoordinateEpoch ) && !mSourceIsDynamic )
? mDestCoordinateEpoch : std::numeric_limits< double >::quiet_NaN();

if ( mSourceIsDynamic && mDestIsDynamic
&& !std::isnan( mSourceCoordinateEpoch ) && mSourceCoordinateEpoch != mDestCoordinateEpoch )
{
// transforms from dynamic crs to dynamic crs with different coordinate epochs are not yet supported by PROJ

}

// init the projections (destination and source)
freeProj();

@@ -103,6 +103,12 @@ class QgsCoordinateTransformPrivate : public QSharedData
bool mShouldReverseCoordinateOperation = false;
bool mAllowFallbackTransforms = true;

bool mSourceIsDynamic = false;
bool mDestIsDynamic = false;
double mSourceCoordinateEpoch = std::numeric_limits< double >::quiet_NaN();
double mDestCoordinateEpoch = std::numeric_limits< double >::quiet_NaN();
double mDefaultTime = std::numeric_limits< double >::quiet_NaN();

//! True if the proj transform corresponds to the reverse direction, and must be flipped when transforming...
bool mIsReversed = false;

@@ -41,6 +41,10 @@ class TestQgsCoordinateTransform: public QObject
void scaleFactor_data();
void transform_data();
void transform();
#if PROJ_VERSION_MAJOR>7 || (PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 2)
void transformEpoch_data();
void transformEpoch();
#endif
void transformLKS();
void transformContextNormalize();
void transform2DPoint();
@@ -355,6 +359,102 @@ void TestQgsCoordinateTransform::transform()
QGSCOMPARENEAR( y, outY, precision );
}

#if PROJ_VERSION_MAJOR>7 || (PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 2)
void TestQgsCoordinateTransform::transformEpoch_data()
{
QTest::addColumn<QgsCoordinateReferenceSystem>( "sourceCrs" );
QTest::addColumn<double>( "sourceEpoch" );
QTest::addColumn<QgsCoordinateReferenceSystem>( "destCrs" );
QTest::addColumn<double>( "destEpoch" );
QTest::addColumn<double>( "x" );
QTest::addColumn<double>( "y" );
QTest::addColumn<int>( "direction" );
QTest::addColumn<double>( "outX" );
QTest::addColumn<double>( "outY" );
QTest::addColumn<double>( "precision" );

QTest::newRow( "GDA2020 to ITRF at central epoch -- no coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2020.0
<< 150.0 << -30.0 << static_cast< int >( QgsCoordinateTransform::ForwardTransform ) << 150.0 << -30.0 << 0.0000000001;

QTest::newRow( "GDA2020 to ITRF at central epoch (reverse) -- no coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2020.0
<< 150.0 << -30.0 << static_cast< int >( QgsCoordinateTransform::ReverseTransform ) << 150.0 << -30.0 << 0.0000000001;

QTest::newRow( "ITRF at central epoch to GDA2020 -- no coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2020.0
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< 150.0 << -30.0 << static_cast< int >( QgsCoordinateTransform::ForwardTransform ) << 150.0 << -30.0 << 0.0000000001;

QTest::newRow( "ITRF at central epoch to GDA2020 (reverse) -- no coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2020.0
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< 150.0 << -30.0 << static_cast< int >( QgsCoordinateTransform::ReverseTransform ) << 150.0 << -30.0 << 0.0000000001;

QTest::newRow( "GDA2020 to ITRF at 2030 -- coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2030.0
<< 150.0 << -30.0 << static_cast< int >( QgsCoordinateTransform::ForwardTransform ) << 150.0000022212 << -29.9999950478 << 0.0000001;

QTest::newRow( "GDA2020 to ITRF at 2030 (reverse) -- coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2030.0
<< 150.0000022212 << -29.9999950478 << static_cast< int >( QgsCoordinateTransform::ReverseTransform ) << 150.0 << -30.0 << 0.0000001;

QTest::newRow( "ITRF at 2030 to GDA2020-- coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2030.0
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< 150.0000022212 << -29.9999950478 << static_cast< int >( QgsCoordinateTransform::ForwardTransform ) << 150.0 << -30.0 << 0.0000001;

QTest::newRow( "ITRF at 2030 to GDA2020 (reverse) -- coordinate change expected" )
<< QgsCoordinateReferenceSystem::fromEpsgId( 9000 ) // ITRF2014
<< 2030.0
<< QgsCoordinateReferenceSystem::fromEpsgId( 7844 ) // GDA2020
<< std::numeric_limits< double >::quiet_NaN()
<< 150.0 << -30.0 << static_cast< int >( QgsCoordinateTransform::ReverseTransform ) << 150.0000022212 << -29.9999950478 << 0.0000001;
}

void TestQgsCoordinateTransform::transformEpoch()
{
QFETCH( QgsCoordinateReferenceSystem, sourceCrs );
QFETCH( double, sourceEpoch );
QFETCH( QgsCoordinateReferenceSystem, destCrs );
QFETCH( double, destEpoch );
QFETCH( double, x );
QFETCH( double, y );
QFETCH( int, direction );
QFETCH( double, outX );
QFETCH( double, outY );
QFETCH( double, precision );

sourceCrs.setCoordinateEpoch( sourceEpoch );
destCrs.setCoordinateEpoch( destEpoch );

double z = 0;
QgsCoordinateTransform ct( sourceCrs, destCrs, QgsProject::instance() );

ct.transformInPlace( x, y, z, static_cast< QgsCoordinateTransform::TransformDirection >( direction ) );
QGSCOMPARENEAR( x, outX, precision );
QGSCOMPARENEAR( y, outY, precision );
}
#endif

void TestQgsCoordinateTransform::transformBoundingBox()
{
//test transforming a bounding box which crosses the 180 degree longitude line

0 comments on commit bd26ae1

Please sign in to comment.