Skip to content
Permalink
Browse files

add option to export mesh dataset to raster block (for processing algs)

  • Loading branch information
PeterPetrik committed Jan 21, 2019
1 parent 8db14d8 commit 7423a74915effac8971948e4bf1950cafcbba44e
@@ -0,0 +1,56 @@
/************************************************************************
* This file has been generated automatically from *
* *
* src/core/mesh/qgsmeshlayerinterpolator.h *
* *
* Do not edit manually ! Edit header and run scripts/sipify.pl again *
************************************************************************/






%ModuleHeaderCode
#include "qgsmeshlayerinterpolator.h"
%End




namespace QgsMeshUtils
{

QgsRasterBlock *exportRasterBlock(
const QgsMeshLayer &layer,
const QgsMeshDatasetIndex &datasetIndex,
const QgsCoordinateReferenceSystem &destination,
const QgsCoordinateTransformContext &context,
double mapUnitsPerPixel,
const QgsRectangle &extent,
QgsRasterBlockFeedback *feedback = 0
) /Factory/;
%Docstring
Exports mesh layer's dataset values as raster block

:param layer: mesh layer
:param datasetIndex: index from layer defining group and dataset (time) to export
:param destination: destination/map CRS. Used to create triangular mesh from native mesh
:param context: Transform context to transform layer CRS to destination CRS
:param mapUnitsPerPixel: map units per pixel for block
:param extent: extent of block in destination CRS
:param feedback: optional raster feedback object for cancelation/preview

:return: raster block with Float.64 values. None on error

.. versionadded:: 3.6
%End
};

/************************************************************************
* This file has been generated automatically from *
* *
* src/core/mesh/qgsmeshlayerinterpolator.h *
* *
* Do not edit manually ! Edit header and run scripts/sipify.pl again *
************************************************************************/
@@ -236,6 +236,7 @@
%Include auto_generated/raster/qgssinglebandgrayrenderer.sip
%Include auto_generated/raster/qgssinglebandpseudocolorrenderer.sip
%Include auto_generated/raster/qgshillshaderenderer.sip
%Include auto_generated/mesh/qgsmeshlayerinterpolator.sip
%Include auto_generated/mesh/qgsmeshrenderersettings.sip
%Include auto_generated/scalebar/qgsdoubleboxscalebarrenderer.sip
%Include auto_generated/scalebar/qgsnumericscalebarrenderer.sip
@@ -18,15 +18,21 @@
///@cond PRIVATE

#include <memory>
#include <limits>

#include "qgsmeshlayerinterpolator.h"

#include "qgis.h"
#include "qgsrasterinterface.h"
#include "qgsmaptopixel.h"
#include "qgsvector.h"
#include "qgspoint.h"
#include "qgspointxy.h"
#include "qgsmeshlayerutils.h"
#include "qgsmeshlayer.h"
#include "qgscoordinatetransformcontext.h"
#include "qgscoordinatetransform.h"
#include "qgsmeshdataprovider.h"

QgsMeshLayerInterpolator::QgsMeshLayerInterpolator( const QgsTriangularMesh &m,
const QVector<double> &datasetValues, const QgsMeshDataBlock &activeFaceFlagValues,
@@ -63,6 +69,8 @@ int QgsMeshLayerInterpolator::bandCount() const
QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent, int width, int height, QgsRasterBlockFeedback *feedback )
{
std::unique_ptr<QgsRasterBlock> outputBlock( new QgsRasterBlock( Qgis::Float64, width, height ) );
const double noDataValue = std::numeric_limits<double>::min();
outputBlock->setNoDataValue( noDataValue );
outputBlock->setIsNoData(); // assume initially that all values are unset
double *data = reinterpret_cast<double *>( outputBlock->bits() );

@@ -142,3 +150,66 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
}

///@endcond

QgsRasterBlock *QgsMeshUtils::exportRasterBlock(
const QgsMeshLayer &layer,
const QgsMeshDatasetIndex &datasetIndex,
const QgsCoordinateReferenceSystem &destination,
const QgsCoordinateTransformContext &context,
double mapUnitsPerPixel,
const QgsRectangle &extent,
QgsRasterBlockFeedback *feedback )
{
if ( !layer.dataProvider() )
return nullptr;

if ( !datasetIndex.isValid() )
return nullptr;

int widthPixel = static_cast<int>( extent.width() / mapUnitsPerPixel );
int heightPixel = static_cast<int>( extent.height() / mapUnitsPerPixel );

const QgsPointXY center = extent.center();
QgsMapToPixel mapToPixel( mapUnitsPerPixel,
center.x(),
center.y(),
widthPixel,
heightPixel,
0 );
QgsCoordinateTransform transform( layer.crs(), destination, context );

QgsRenderContext renderContext;
renderContext.setCoordinateTransform( transform );
renderContext.setMapToPixel( mapToPixel );
renderContext.setExtent( extent );

std::unique_ptr<QgsMesh> nativeMesh = qgis::make_unique<QgsMesh>();
layer.dataProvider()->populateMesh( nativeMesh.get() );
std::unique_ptr<QgsTriangularMesh> triangularMesh = qgis::make_unique<QgsTriangularMesh>();
triangularMesh->update( nativeMesh.get(), &renderContext );

const QgsMeshDatasetGroupMetadata metadata = layer.dataProvider()->datasetGroupMetadata( datasetIndex );
bool scalarDataOnVertices = metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;

QgsMeshDataBlock vals = layer.dataProvider()->datasetValues(
datasetIndex,
0,
scalarDataOnVertices ? nativeMesh->vertices.count() : nativeMesh->faces.count() );

QVector<double> datasetValues = QgsMeshLayerUtils::calculateMagnitudes( vals );
QgsMeshDataBlock activeFaceFlagValues = layer.dataProvider()->areFacesActive(
datasetIndex,
0,
nativeMesh->faces.count() );

QgsMeshLayerInterpolator iterpolator(
*( triangularMesh.get() ),
datasetValues,
activeFaceFlagValues,
scalarDataOnVertices,
renderContext,
QSize( widthPixel, heightPixel )
);

return iterpolator.block( 0, extent, widthPixel, heightPixel, feedback );
}
@@ -20,20 +20,27 @@

class QgsMeshLayer;
class QgsSymbol;

#define SIP_NO_FILE

#include <QSize>
class QgsCoordinateReferenceSystem;
class QgsCoordinateTransformContext;
class QgsMeshDatasetIndex;

#include "qgis.h"
#include "qgis_sip.h"

#include <QSize>
#include "qgsmaplayerrenderer.h"
#include "qgsrendercontext.h"
#include "qgstriangularmesh.h"
#include "qgsrasterinterface.h"
#include "qgssinglebandpseudocolorrenderer.h"
#include "qgsrastershader.h"

#ifdef SIP_RUN
% ModuleHeaderCode
#include "qgsmeshlayerinterpolator.h"
% End
#endif

///@cond PRIVATE

/**
@@ -43,7 +50,7 @@ class QgsSymbol;
* \note not available in Python bindings
* \since QGIS 3.2
*/
class QgsMeshLayerInterpolator : public QgsRasterInterface
class QgsMeshLayerInterpolator : public QgsRasterInterface SIP_SKIP
{
public:
//! Ctor
@@ -71,4 +78,32 @@ class QgsMeshLayerInterpolator : public QgsRasterInterface

///@endcond

namespace QgsMeshUtils
{

/**
* Exports mesh layer's dataset values as raster block
*
* \param layer mesh layer
* \param datasetIndex index from layer defining group and dataset (time) to export
* \param destination destination/map CRS. Used to create triangular mesh from native mesh
* \param context Transform context to transform layer CRS to destination CRS
* \param mapUnitsPerPixel map units per pixel for block
* \param extent extent of block in destination CRS
* \param feedback optional raster feedback object for cancelation/preview
* \returns raster block with Float::64 values. nullptr on error
*
* \since QGIS 3.6
*/
CORE_EXPORT QgsRasterBlock *exportRasterBlock(
const QgsMeshLayer &layer,
const QgsMeshDatasetIndex &datasetIndex,
const QgsCoordinateReferenceSystem &destination,
const QgsCoordinateTransformContext &context,
double mapUnitsPerPixel,
const QgsRectangle &extent,
QgsRasterBlockFeedback *feedback = nullptr
) SIP_FACTORY;
};

#endif // QGSMESHLAYERINTERPOLATOR_H
@@ -159,6 +159,9 @@ SET(TESTS
testqgsmaptopixelgeometrysimplifier.cpp
testqgsmaptopixel.cpp
testqgsmarkerlinesymbol.cpp
testqgsmeshlayer.cpp
testqgsmeshlayerinterpolator.cpp
testqgsmeshlayerrenderer.cpp
testqgsnetworkcontentfetcher.cpp
testqgsogcutils.cpp
testqgsogrutils.cpp
@@ -206,8 +209,6 @@ SET(TESTS
testqgsvectorlayerutils.cpp
testqgsziputils.cpp
testziplayer.cpp
testqgsmeshlayer.cpp
testqgsmeshlayerrenderer.cpp
testqgslayerdefinition.cpp
testqgssqliteutils.cpp
testqgsmimedatautils.cpp
@@ -0,0 +1,93 @@
/***************************************************************************
testqgsmeshlayerinterpolator.cpp
--------------------------------
Date : January 2019
Copyright : (C) 2019 by Peter Petrik
Email : zilolv at gmail dot com
***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#include "qgstest.h"
#include <QObject>
#include <QString>
#include <QApplication>

//qgis includes...
#include <qgsrasterblock.h>
#include <qgsmeshlayer.h>
#include <qgsmeshlayerinterpolator.h>
#include <qgsapplication.h>
#include <qgscoordinatereferencesystem.h>
#include <qgsproject.h>

/**
* \ingroup UnitTests
* This is a unit test for the TestQgsMeshLayerInterpolator class
*/
class TestQgsMeshLayerInterpolator: public QObject
{
Q_OBJECT
private slots:
void initTestCase();// will be called before the first testfunction is executed.
void cleanupTestCase();// will be called after the last testfunction was executed.
void init() {} // will be called before each testfunction is executed.
void cleanup() {} // will be called after every testfunction.

void testExportRasterBand();
private:
QString mTestDataDir;
};

//runs before all tests
void TestQgsMeshLayerInterpolator::initTestCase()
{
// init QGIS's paths - true means that all path will be inited from prefix
QgsApplication::init();
QgsApplication::initQgis();
QgsApplication::showSettings();
mTestDataDir = QStringLiteral( TEST_DATA_DIR ); //defined in CmakeLists.txt
}
//runs after all tests
void TestQgsMeshLayerInterpolator::cleanupTestCase()
{
QgsApplication::exitQgis();
}

void TestQgsMeshLayerInterpolator::testExportRasterBand()
{
QgsMeshLayer memoryLayer( mTestDataDir + "/mesh/quad_and_triangle.2dm",
"Triangle and Quad Mdal",
"mdal" );
QVERIFY( memoryLayer.isValid() );
QgsMeshDatasetIndex index( 0, 0 ); // bed elevation
memoryLayer.setCrs( QgsCoordinateReferenceSystem::fromEpsgId( 27700 ) );

QgsRasterBlock *block = QgsMeshUtils::exportRasterBlock(
memoryLayer,
index,
memoryLayer.crs(),
QgsProject::instance()->transformContext(),
100,
memoryLayer.extent()
);

QCOMPARE( block->width(), 20 );
QCOMPARE( block->height(), 10 );
QVERIFY( block->isValid() );
QVERIFY( !block->isEmpty() );
QCOMPARE( block->dataType(), Qgis::Float64 );
QVERIFY( block->hasNoDataValue() );
QVERIFY( block->hasNoData() );

QCOMPARE( block->value( 0, 0 ), 10.0 );
QCOMPARE( block->value( 5, 5 ), 35.0 );
QVERIFY( block->isNoData( 10, 10 ) );
}

QGSTEST_MAIN( TestQgsMeshLayerInterpolator )
#include "testqgsmeshlayerinterpolator.moc"

0 comments on commit 7423a74

Please sign in to comment.
You can’t perform that action at this time.