Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Tile Matrix in different CRS
Add a new static method to QgsTileMatrix to get tile matrix in different CRS with top left point and scale denominator for 0 zoom level.
  • Loading branch information
rldhont committed Sep 10, 2021
1 parent 9adc3aa commit 95328ee
Show file tree
Hide file tree
Showing 4 changed files with 231 additions and 11 deletions.
24 changes: 23 additions & 1 deletion python/core/auto_generated/qgstiles.sip.in
Expand Up @@ -108,9 +108,26 @@ Please note that we follow the XYZ convention of X/Y axes, i.e. top-left tile ha
%End
public:

static QgsTileMatrix fromWebMercator( int mZoomLevel );
static QgsTileMatrix fromWebMercator( int zoomLevel );
%Docstring
Returns a tile matrix for the usual web mercator
%End

static QgsTileMatrix fromCustomDef( int zoomLevel, const QgsCoordinateReferenceSystem &crs,
const QgsPointXY &z0TopLeftPoint, double z0Dimension,
int z0MatrixWidth = 1, int z0MatrixHeight = 1 );
%Docstring
Returns a tile matrix for a specific CRS, top left point, zoom level 0 dimension in CRS units
%End

static QgsTileMatrix fromTileMatrix( const int &zoomLevel, const QgsTileMatrix &tileMatrix );
%Docstring
Returns a tile matrix based on another one
%End

QgsCoordinateReferenceSystem crs() const;
%Docstring
Returns the authority identifier for the CRS of the tile matrix
%End

int zoomLevel() const;
Expand Down Expand Up @@ -156,6 +173,11 @@ Returns tile range that fully covers the given extent
QPointF mapToTileCoordinates( const QgsPointXY &mapPoint ) const;
%Docstring
Returns row/column coordinates (floating point number) from the given point in map coordinates
%End

bool isRootTileMatrix() const;
%Docstring
Returns the root status of the tile matrix (zoom level == 0)
%End

};
Expand Down
54 changes: 48 additions & 6 deletions src/core/qgstiles.cpp
Expand Up @@ -16,22 +16,64 @@
#include "qgstiles.h"

#include "qgslogger.h"
#include "qgscoordinatereferencesystem.h"
#include "qgssettings.h"

QgsTileMatrix QgsTileMatrix::fromWebMercator( int zoomLevel )
{
double z0xMin = -20037508.3427892, z0yMax = 20037508.3427892;

return fromCustomDef( zoomLevel, QgsCoordinateReferenceSystem( "EPSG:3857" ), QgsPointXY( z0xMin, z0yMax ), 2 * z0yMax );
}

QgsTileMatrix QgsTileMatrix::fromCustomDef( int zoomLevel, const QgsCoordinateReferenceSystem &crs,
const QgsPointXY &z0TopLeftPoint, double z0Dimension, int z0MatrixWidth, int z0MatrixHeight )
{
// Square extent calculation
double z0xMin = z0TopLeftPoint.x();
double z0yMax = z0TopLeftPoint.y();
double z0xMax = z0xMin + z0MatrixWidth * z0Dimension;
double z0yMin = z0yMax - z0MatrixHeight * z0Dimension;

// Constant for scale denominator calculation
const double tileSize = 256.0;
const double PIXELS_TO_M = 2.8 / 10000.0; // WMS/WMTS define "standardized rendering pixel size" as 0.28mm
const double UNIT_TO_M = QgsUnitTypes::fromUnitToUnitFactor( crs.mapUnits(), QgsUnitTypes::DistanceMeters );
// Scale denominator calculation
double scaleDenom0 = ( z0Dimension / tileSize ) * ( UNIT_TO_M / PIXELS_TO_M );

int numTiles = static_cast<int>( pow( 2, zoomLevel ) ); // assuming we won't ever go over 30 zoom levels
double z0xMin = -20037508.3427892, z0yMin = -20037508.3427892;
double z0xMax = 20037508.3427892, z0yMax = 20037508.3427892;
double s0 = 559082264.0287178; // scale denominator at zoom level 0 of GoogleCRS84Quad

QgsTileMatrix tm;
tm.mCrs = crs;
tm.mZoomLevel = zoomLevel;
tm.mMatrixWidth = numTiles;
tm.mMatrixHeight = numTiles;
tm.mMatrixWidth = z0MatrixWidth * numTiles;
tm.mMatrixHeight = z0MatrixHeight * numTiles;
tm.mTileXSpan = ( z0xMax - z0xMin ) / tm.mMatrixWidth;
tm.mTileYSpan = ( z0yMax - z0yMin ) / tm.mMatrixHeight;
tm.mExtent = QgsRectangle( z0xMin, z0yMin, z0xMax, z0yMax );
tm.mScaleDenom = s0 / pow( 2, zoomLevel );
tm.mScaleDenom = scaleDenom0 / pow( 2, zoomLevel );
return tm;
}


QgsTileMatrix QgsTileMatrix::fromTileMatrix( const int &zoomLevel, const QgsTileMatrix &tileMatrix )
{
QgsTileMatrix tm;
int numTiles = static_cast<int>( pow( 2, zoomLevel ) ); // assuming we won't ever go over 30 zoom levels
int aZoomLevel = tileMatrix.zoomLevel();
int aNumTiles = static_cast<int>( pow( 2, aZoomLevel ) );
int aMatrixWidth = tileMatrix.matrixWidth();
int aMatrixHeight = tileMatrix.matrixHeight();
QgsRectangle aExtent = tileMatrix.extent();
tm.mCrs = tileMatrix.crs();
tm.mZoomLevel = zoomLevel;
tm.mMatrixWidth = aMatrixWidth * numTiles / aNumTiles;
tm.mMatrixHeight = aMatrixHeight * numTiles / aNumTiles;
tm.mTileXSpan = aExtent.width() / tm.mMatrixWidth;
tm.mTileYSpan = aExtent.height() / tm.mMatrixHeight;
tm.mExtent = aExtent;
tm.mScaleDenom = tileMatrix.scale() * pow( 2, aZoomLevel ) / pow( 2, zoomLevel );
return tm;
}

Expand Down
21 changes: 19 additions & 2 deletions src/core/qgstiles.h
Expand Up @@ -20,6 +20,7 @@
#include "qgis_sip.h"

#include "qgsrectangle.h"
#include "qgscoordinatereferencesystem.h"

/**
* \ingroup core
Expand Down Expand Up @@ -104,7 +105,18 @@ class CORE_EXPORT QgsTileMatrix
public:

//! Returns a tile matrix for the usual web mercator
static QgsTileMatrix fromWebMercator( int mZoomLevel );
static QgsTileMatrix fromWebMercator( int zoomLevel );

//! Returns a tile matrix for a specific CRS, top left point, zoom level 0 dimension in CRS units
static QgsTileMatrix fromCustomDef( int zoomLevel, const QgsCoordinateReferenceSystem &crs,
const QgsPointXY &z0TopLeftPoint, double z0Dimension,
int z0MatrixWidth = 1, int z0MatrixHeight = 1 );

//! Returns a tile matrix based on another one
static QgsTileMatrix fromTileMatrix( const int &zoomLevel, const QgsTileMatrix &tileMatrix );

//! Returns the authority identifier for the CRS of the tile matrix
QgsCoordinateReferenceSystem crs() const { return mCrs; }

//! Returns zoom level of the tile matrix
int zoomLevel() const { return mZoomLevel; }
Expand Down Expand Up @@ -133,9 +145,14 @@ class CORE_EXPORT QgsTileMatrix
//! Returns row/column coordinates (floating point number) from the given point in map coordinates
QPointF mapToTileCoordinates( const QgsPointXY &mapPoint ) const;

//! Returns the root status of the tile matrix (zoom level == 0)
bool isRootTileMatrix() const { return mZoomLevel == 0; }

private:
//! Crs associated with the tile matrix
QgsCoordinateReferenceSystem mCrs;
//! Zoom level index associated with the tile matrix
int mZoomLevel;
int mZoomLevel = -1;
//! Number of columns of the tile matrix
int mMatrixWidth;
//! Number of rows of the tile matrix
Expand Down
143 changes: 141 additions & 2 deletions tests/src/core/testqgstiles.cpp
Expand Up @@ -44,6 +44,8 @@ class TestQgsTiles : public QObject
void cleanup() {} // will be called after every testfunction.

void test_matrixFromWebMercator();
void test_matrixFromCustomDef();
void test_matrixFromTileMatrix();
};


Expand All @@ -65,24 +67,26 @@ void TestQgsTiles::cleanupTestCase()
void TestQgsTiles::test_matrixFromWebMercator()
{
QgsTileMatrix tm = QgsTileMatrix::fromWebMercator( 0 );
QCOMPARE( tm.crs().authid(), "EPSG:3857" );
QCOMPARE( tm.zoomLevel(), 0 );
QCOMPARE( tm.matrixWidth(), 1 );
QCOMPARE( tm.matrixHeight(), 1 );
QVERIFY( qgsDoubleNear( tm.extent().xMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().xMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 559082264.0287178, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 559082264.0287178, 0.00001 ) );

tm = QgsTileMatrix::fromWebMercator( 3 );
QCOMPARE( tm.crs().authid(), "EPSG:3857" );
QCOMPARE( tm.zoomLevel(), 3 );
QCOMPARE( tm.matrixWidth(), 8 );
QCOMPARE( tm.matrixHeight(), 8 );
QVERIFY( qgsDoubleNear( tm.extent().xMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().xMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 559082264.0287178 / 8.0, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 559082264.0287178 / 8.0, 0.00001 ) );

QgsRectangle te = tm.tileExtent( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( te.xMinimum(), -5009377.0856973, 0.0000001 ) );
Expand All @@ -108,6 +112,141 @@ void TestQgsTiles::test_matrixFromWebMercator()
QCOMPARE( tr.endRow(), 3 );
}

void TestQgsTiles::test_matrixFromCustomDef()
{
// Try to build the Web Mercator tile matrix
QgsTileMatrix tm = QgsTileMatrix::fromCustomDef( 0, QgsCoordinateReferenceSystem( "EPSG:3857" ), QgsPointXY( -20037508.3427892, 20037508.3427892 ), 40075016.6855784 );
QCOMPARE( tm.crs().authid(), "EPSG:3857" );
QCOMPARE( tm.zoomLevel(), 0 );
QCOMPARE( tm.matrixWidth(), 1 );
QCOMPARE( tm.matrixHeight(), 1 );
QVERIFY( qgsDoubleNear( tm.extent().xMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().xMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 559082264.0287178, 0.00001 ) );

tm = QgsTileMatrix::fromCustomDef( 3, QgsCoordinateReferenceSystem( "EPSG:3857" ), QgsPointXY( -20037508.3427892, 20037508.3427892 ), 40075016.6855784 );
QCOMPARE( tm.crs().authid(), "EPSG:3857" );
QCOMPARE( tm.zoomLevel(), 3 );
QCOMPARE( tm.matrixWidth(), 8 );
QCOMPARE( tm.matrixHeight(), 8 );
QVERIFY( qgsDoubleNear( tm.extent().xMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().xMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 559082264.0287178 / 8.0, 0.00001 ) );

QgsRectangle te = tm.tileExtent( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( te.xMinimum(), -5009377.0856973, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.yMinimum(), 0, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.xMaximum(), 0, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.yMaximum(), 5009377.0856973, 0.0000001 ) );

QgsPointXY tc = tm.tileCenter( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( tc.x(), te.xMinimum() + te.width() / 2, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.x(), -2504688.54284865, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.y(), te.yMinimum() + te.height() / 2, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.y(), 2504688.54284865, 0.0000001 ) );

QgsVectorLayer *vl = new QgsVectorLayer( mDataDir + "/points.shp", "points", "ogr" );
QgsCoordinateReferenceSystem destCrs( "EPSG:3857" );
QgsCoordinateTransform ct( vl->crs(), destCrs, QgsCoordinateTransformContext() );
QgsRectangle r = ct.transformBoundingBox( vl->extent() );
QgsTileRange tr = tm.tileRangeFromExtent( r );
QVERIFY( tr.isValid() );
QCOMPARE( tr.startColumn(), 1 );
QCOMPARE( tr.endColumn(), 2 );
QCOMPARE( tr.startRow(), 2 );
QCOMPARE( tr.endRow(), 3 );

// Try to build the World Global tile matrix
tm = QgsTileMatrix::fromCustomDef( 0, QgsCoordinateReferenceSystem( "EPSG:4326" ), QgsPointXY( -180.0, 90.0 ), 180.0, 2, 1 );
QCOMPARE( tm.crs().authid(), "EPSG:4326" );
QCOMPARE( tm.zoomLevel(), 0 );
QCOMPARE( tm.matrixWidth(), 2 );
QCOMPARE( tm.matrixHeight(), 1 );
QVERIFY( qgsDoubleNear( tm.extent().xMinimum(), -180.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMinimum(), -90.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.extent().xMaximum(), 180.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMaximum(), 90.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 279541132.014, 0.001 ) );

tm = QgsTileMatrix::fromCustomDef( 3, QgsCoordinateReferenceSystem( "EPSG:4326" ), QgsPointXY( -180.0, 90.0 ), 180.0, 2, 1 );
QCOMPARE( tm.crs().authid(), "EPSG:4326" );
QCOMPARE( tm.zoomLevel(), 3 );
QCOMPARE( tm.matrixWidth(), 16 );
QCOMPARE( tm.matrixHeight(), 8 );
QVERIFY( qgsDoubleNear( tm.extent().xMinimum(), -180.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMinimum(), -90.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.extent().xMaximum(), 180.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.extent().yMaximum(), 90.0, 0.1 ) );
QVERIFY( qgsDoubleNear( tm.scale(), 279541132.014 / 8.0, 0.001 ) );

te = tm.tileExtent( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( te.xMinimum(), -112.5, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.yMinimum(), 0, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.xMaximum(), -90, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.yMaximum(), 22.5, 0.0000001 ) );

tc = tm.tileCenter( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( tc.x(), te.xMinimum() + te.width() / 2, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.x(), -101.25, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.y(), te.yMinimum() + te.height() / 2, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.y(), 11.25, 0.0000001 ) );

destCrs = QgsCoordinateReferenceSystem( "EPSG:4326" );
ct = QgsCoordinateTransform( vl->crs(), destCrs, QgsCoordinateTransformContext() );
r = ct.transformBoundingBox( vl->extent() );
tr = tm.tileRangeFromExtent( r );
QVERIFY( tr.isValid() );
QCOMPARE( tr.startColumn(), 2 );
QCOMPARE( tr.endColumn(), 4 );
QCOMPARE( tr.startRow(), 1 );
QCOMPARE( tr.endRow(), 2 );
}

void TestQgsTiles::test_matrixFromTileMatrix()
{
QgsTileMatrix tm = QgsTileMatrix::fromWebMercator( 0 );

// Build level 3 from 0
QgsTileMatrix tm3 = QgsTileMatrix::fromTileMatrix( 3, tm );
QCOMPARE( tm3.crs().authid(), "EPSG:3857" );
QCOMPARE( tm3.zoomLevel(), 3 );
QCOMPARE( tm3.matrixWidth(), 8 );
QCOMPARE( tm3.matrixHeight(), 8 );
QVERIFY( qgsDoubleNear( tm3.extent().xMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm3.extent().yMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm3.extent().xMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm3.extent().yMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm3.scale(), 559082264.0287178 / 8.0, 0.00001 ) );

QgsRectangle te = tm3.tileExtent( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( te.xMinimum(), -5009377.0856973, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.yMinimum(), 0, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.xMaximum(), 0, 0.0000001 ) );
QVERIFY( qgsDoubleNear( te.yMaximum(), 5009377.0856973, 0.0000001 ) );

QgsPointXY tc = tm3.tileCenter( QgsTileXYZ( 3, 3, 3 ) );
QVERIFY( qgsDoubleNear( tc.x(), te.xMinimum() + te.width() / 2, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.x(), -2504688.54284865, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.y(), te.yMinimum() + te.height() / 2, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tc.y(), 2504688.54284865, 0.0000001 ) );

// Build level 0 from 3
QgsTileMatrix tm0 = QgsTileMatrix::fromTileMatrix( 0, tm3 );
QCOMPARE( tm0.crs().authid(), "EPSG:3857" );
QCOMPARE( tm0.zoomLevel(), 0 );
QCOMPARE( tm0.matrixWidth(), 1 );
QCOMPARE( tm0.matrixHeight(), 1 );
QVERIFY( qgsDoubleNear( tm0.extent().xMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm0.extent().yMinimum(), -20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm0.extent().xMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm0.extent().yMaximum(), 20037508.3427892, 0.0000001 ) );
QVERIFY( qgsDoubleNear( tm0.scale(), 559082264.0287178, 0.00001 ) );
}


QGSTEST_MAIN( TestQgsTiles )
#include "testqgstiles.moc"

0 comments on commit 95328ee

Please sign in to comment.