Skip to content

Commit

Permalink
Make root bounding volume accessible from tiled mesh data provider,
Browse files Browse the repository at this point in the history 
fix extent/crs handling for sphere/box bounding volume types

There's a tricky consideration here -- the default crs for sphere/
box bounding volumes is EPSG:4978. But we CANNOT transform from
EPSG:4978 without using valid z values (otherwise we get nonsense
results). This raises an issue with all the various QgsMapLayer
methods which use QgsMapLayer::crs to transform layer properties,
eg the layer's extent. If we advertise the layer (correctly) as
EPSG:4978, then the layer's extent will be incorrectly calculated
in 1000s of places in QGIS (and we are violating the api design
of these methods!)

So we handle this by always advertising these sources as EPSG:4979
(ie a geographic degrees based CRS of WGS84+height), and transform
the source's bounding volume to EPSG:4979 within the data provider,
considering correctly the source's z values, and use THIS as the
data provider's extent.

Then we add a new method to the QgsTiledMeshDataProvider interface
for meshCrs(), which returns the ACTUAL crs that the mesh
geometries are in (ie. EPSG:4978).

Care must be taken to use the correct choice of the advertised crs()
vs meshCrs() in methods which interact with tiled mesh data.
  • Loading branch information
@nyalldawson
nyalldawson committed Jul 17, 2023
1 parent 541e981 commit 7e35792
Show file tree
Hide file tree
Showing 7 changed files with 212 additions and 64 deletions.
25 changes: 25 additions & 0 deletions python/core/auto_generated/tiledmesh/qgstiledmeshdataprovider.sip.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@




class QgsTiledMeshDataProvider: QgsDataProvider
{
%Docstring(signature="appended")
Expand Down Expand Up @@ -54,6 +55,30 @@ Returns metadata in a format suitable for feeding directly
into a subset of the GUI properties "Metadata" tab.
%End

virtual const QgsCoordinateReferenceSystem meshCrs() const = 0;
%Docstring
Returns the original coordinate reference system for the tiled mesh data.

This may differ from the :py:func:`QgsDataProvider.crs()`, which is the best CRS representation
for the data provider for 2D use.

.. warning::

Care must be taken to ensure that :py:func:`~QgsTiledMeshDataProvider.meshCrs` is used instead of :py:func:`~QgsTiledMeshDataProvider.crs` whenever
transforming bounding volumes or geometries associated with the provider.
%End

virtual const QgsAbstractTiledMeshNodeBoundingVolume *boundingVolume() const = 0;
%Docstring
Returns the bounding volume for the data provider.

This corresponds to the root node bounding volume.

.. warning::

Coordinates in the returned volume are in the :py:func:`~QgsTiledMeshDataProvider.meshCrs` reference system, not the :py:func:`QgsDataProvider.crs()` system.
%End

};

/************************************************************************
Expand Down
100 changes: 72 additions & 28 deletions src/core/tiledmesh/qgscesiumtilesdataprovider.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@
#include "qgssphere.h"
#include "qgslogger.h"
#include "qgsorientedbox3d.h"
#include "qgstiledmeshboundingvolume.h"
#include "qgscoordinatetransform.h"

#include <QUrl>
#include <QIcon>
Expand All @@ -45,58 +47,88 @@

QgsCesiumTilesDataProviderSharedData::QgsCesiumTilesDataProviderSharedData() = default;

void QgsCesiumTilesDataProviderSharedData::setTilesetContent( const QString &tileset )
void QgsCesiumTilesDataProviderSharedData::setTilesetContent( const QString &tileset, const QgsCoordinateTransformContext &transformContext )
{
mTileset = json::parse( tileset.toStdString() );

// parse root
{
const auto &root = mTileset[ "root" ];
// parse root bounding volume

// TODO -- read crs from metadata tags. Need to find real world examples of this. And can metadata crs override
// the EPSG:4979 requirement from a region bounding volume??

{
const auto &rootBoundingVolume = root[ "boundingVolume" ];
if ( rootBoundingVolume.contains( "region" ) )
{
const QgsBox3d rootRegion = QgsCesiumUtils::parseRegion( rootBoundingVolume[ "region" ] );
if ( !rootRegion.isNull() )
{
mBoundingVolume = std::make_unique< QgsTiledMeshNodeBoundingVolumeRegion >( rootRegion );
mZRange = QgsDoubleRange( rootRegion.zMinimum(), rootRegion.zMaximum() );
// The latitude and longitude values are given in radians!
// TODO -- is this ALWAYS the case? What if there's a region root bounding volume, but a transform object present? What if there's crs metadata specifying a different crs?
mCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4979" ) );
const double xMin = rootRegion.xMinimum() * 180 / M_PI;
const double xMax = rootRegion.xMaximum() * 180 / M_PI;
const double yMin = rootRegion.yMinimum() * 180 / M_PI;
const double yMax = rootRegion.yMaximum() * 180 / M_PI;
mExtent = QgsRectangle( xMin, yMin, xMax, yMax );
mLayerCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4979" ) );
mMeshCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4979" ) );
mExtent = rootRegion.toRectangle();
}
}
else if ( rootBoundingVolume.contains( "box" ) )
{
const QgsOrientedBox3D bbox = QgsCesiumUtils::parseBox( rootBoundingVolume["box"] );
if ( !bbox.isNull() )
{
const QgsBox3d rootRegion = bbox.extent();
mZRange = QgsDoubleRange( rootRegion.zMinimum(), rootRegion.zMaximum() );
const double xMin = rootRegion.xMinimum();
const double xMax = rootRegion.xMaximum();
const double yMin = rootRegion.yMinimum();
const double yMax = rootRegion.yMaximum();
mExtent = QgsRectangle( xMin, yMin, xMax, yMax );
// layer must advertise as EPSG:4979, as the various QgsMapLayer
// methods which utilise QgsMapLayer::crs() (such as layer extent transformation)
// are all purely 2D and can't handle the cesium data source z value
// range in EPSG:4978 !
mLayerCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4979" ) );
mMeshCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4978" ) );

const QgsCoordinateTransform transform( mMeshCrs, mLayerCrs, transformContext );

mBoundingVolume = std::make_unique< QgsTiledMeshNodeBoundingVolumeBox >( bbox );
try
{
const QgsBox3d rootRegion = mBoundingVolume->bounds( transform );
mZRange = QgsDoubleRange( rootRegion.zMinimum(), rootRegion.zMaximum() );

std::unique_ptr< QgsAbstractGeometry > extent2D( mBoundingVolume->as2DGeometry( transform ) );
mExtent = extent2D->boundingBox();
}
catch ( QgsCsException & )
{
QgsDebugError( QStringLiteral( "Caught transform exception when transforming boundingVolume" ) );
}
}
}
else if ( rootBoundingVolume.contains( "sphere" ) )
{
const QgsSphere sphere = QgsCesiumUtils::parseSphere( rootBoundingVolume["sphere"] );
if ( !sphere.isNull() )
{
const QgsBox3d rootRegion = sphere.boundingBox();
mZRange = QgsDoubleRange( rootRegion.zMinimum(), rootRegion.zMaximum() );
const double xMin = rootRegion.xMinimum();
const double xMax = rootRegion.xMaximum();
const double yMin = rootRegion.yMinimum();
const double yMax = rootRegion.yMaximum();
mExtent = QgsRectangle( xMin, yMin, xMax, yMax );
// layer must advertise as EPSG:4979, as the various QgsMapLayer
// methods which utilise QgsMapLayer::crs() (such as layer extent transformation)
// are all purely 2D and can't handle the cesium data source z value
// range in EPSG:4978 !
mLayerCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4979" ) );
mMeshCrs = QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4978" ) );

const QgsCoordinateTransform transform( mMeshCrs, mLayerCrs, transformContext );

mBoundingVolume = std::make_unique< QgsTiledMeshNodeBoundingVolumeSphere >( sphere );
try
{
const QgsBox3d rootRegion = mBoundingVolume->bounds( transform );
mZRange = QgsDoubleRange( rootRegion.zMinimum(), rootRegion.zMaximum() );

std::unique_ptr< QgsAbstractGeometry > extent2D( mBoundingVolume->as2DGeometry( transform ) );
mExtent = extent2D->boundingBox();
}
catch ( QgsCsException & )
{
QgsDebugError( QStringLiteral( "Caught transform exception when transforming boundingVolume" ) );
}
}
}
else
Expand All @@ -105,8 +137,6 @@ void QgsCesiumTilesDataProviderSharedData::setTilesetContent( const QString &til
}
}
}
// TODO -- read crs from metadata tags. Need to find real world examples of this. And can metadata crs override
// the EPSG:4979 requirement from a region bounding volume??
}


Expand Down Expand Up @@ -175,7 +205,7 @@ bool QgsCesiumTilesDataProvider::init()
}

const QgsNetworkReplyContent content = networkRequest.reply();
mShared->setTilesetContent( content.content() );
mShared->setTilesetContent( content.content(), transformContext() );
}
else
{
Expand All @@ -186,7 +216,7 @@ bool QgsCesiumTilesDataProvider::init()
if ( file.open( QIODevice::ReadOnly | QIODevice::Text ) )
{
const QByteArray raw = file.readAll();
mShared->setTilesetContent( raw );
mShared->setTilesetContent( raw, transformContext() );
}
else
{
Expand All @@ -206,7 +236,7 @@ QgsCoordinateReferenceSystem QgsCesiumTilesDataProvider::crs() const
{
QGIS_PROTECT_QOBJECT_THREAD_ACCESS

return mShared->mCrs;
return mShared->mLayerCrs;
}

QgsRectangle QgsCesiumTilesDataProvider::extent() const
Expand Down Expand Up @@ -297,6 +327,20 @@ QString QgsCesiumTilesDataProvider::htmlMetadata() const
return metadata;
}

const QgsCoordinateReferenceSystem QgsCesiumTilesDataProvider::meshCrs() const
{
QGIS_PROTECT_QOBJECT_THREAD_ACCESS

return mShared ? mShared->mMeshCrs : QgsCoordinateReferenceSystem();
}

const QgsAbstractTiledMeshNodeBoundingVolume *QgsCesiumTilesDataProvider::boundingVolume() const
{
QGIS_PROTECT_QOBJECT_THREAD_ACCESS

return mShared ? mShared->mBoundingVolume.get() : nullptr;
}


//
// QgsCesiumTilesProviderMetadata
Expand Down
13 changes: 11 additions & 2 deletions src/core/tiledmesh/qgscesiumtilesdataprovider.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,15 +28,21 @@

#define SIP_NO_FILE

class QgsAbstractTiledMeshNodeBoundingVolume;
class QgsCoordinateTransformContext;

///@cond PRIVATE

class QgsCesiumTilesDataProviderSharedData
{
public:
QgsCesiumTilesDataProviderSharedData();
void setTilesetContent( const QString &tileset );
void setTilesetContent( const QString &tileset, const QgsCoordinateTransformContext &transformContext );

QgsCoordinateReferenceSystem mLayerCrs;
QgsCoordinateReferenceSystem mMeshCrs;
std::unique_ptr< QgsAbstractTiledMeshNodeBoundingVolume > mBoundingVolume;

QgsCoordinateReferenceSystem mCrs;
QgsRectangle mExtent;
nlohmann::json mTileset;
QgsDoubleRange mZRange;
Expand Down Expand Up @@ -68,6 +74,9 @@ class CORE_EXPORT QgsCesiumTilesDataProvider final: public QgsTiledMeshDataProvi
QString name() const final;
QString description() const final;
QString htmlMetadata() const final;
const QgsCoordinateReferenceSystem meshCrs() const final;
const QgsAbstractTiledMeshNodeBoundingVolume *boundingVolume() const override;

private:

bool init();
Expand Down
12 changes: 8 additions & 4 deletions src/core/tiledmesh/qgscesiumutils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,16 @@ QgsBox3d QgsCesiumUtils::parseRegion( const json &region )
{
try
{
const double west = region[0].get<double>();
const double south = region[1].get<double>();
const double east = region[2].get<double>();
const double north = region[3].get<double>();
// The latitude and longitude values are given in radians!
// TODO -- is this ALWAYS the case? What if there's a region root bounding volume, but a transform object present? What if there's crs metadata specifying a different crs?

const double west = region[0].get<double>() * 180 / M_PI;
const double south = region[1].get<double>() * 180 / M_PI;
const double east = region[2].get<double>() * 180 / M_PI;
const double north = region[3].get<double>() * 180 / M_PI;
double minHeight = region[4].get<double>();
double maxHeight = region[5].get<double>();

return QgsBox3d( west, south, minHeight, east, north, maxHeight );
}
catch ( nlohmann::json::exception & )
Expand Down
22 changes: 22 additions & 0 deletions src/core/tiledmesh/qgstiledmeshdataprovider.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@
#include "qgsdataprovider.h"
#include "qgis.h"

class QgsAbstractTiledMeshNodeBoundingVolume;

/**
* \ingroup core
* \brief Base class for data providers for QgsTiledMeshLayer
Expand Down Expand Up @@ -68,6 +70,26 @@ class CORE_EXPORT QgsTiledMeshDataProvider: public QgsDataProvider
*/
virtual QString htmlMetadata() const;

/**
* Returns the original coordinate reference system for the tiled mesh data.
*
* This may differ from the QgsDataProvider::crs(), which is the best CRS representation
* for the data provider for 2D use.
*
* \warning Care must be taken to ensure that meshCrs() is used instead of crs() whenever
* transforming bounding volumes or geometries associated with the provider.
*/
virtual const QgsCoordinateReferenceSystem meshCrs() const = 0;

/**
* Returns the bounding volume for the data provider.
*
* This corresponds to the root node bounding volume.
*
* \warning Coordinates in the returned volume are in the meshCrs() reference system, not the QgsDataProvider::crs() system.
*/
virtual const QgsAbstractTiledMeshNodeBoundingVolume *boundingVolume() const = 0;

};

#endif // QGSTILEDMESHDATAPROVIDER_H
Loading

0 comments on commit 7e35792

Please sign in to comment.