Skip to content

Commit

Permalink
rasterize implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@vcloarec @PeterPetrik
vcloarec authored and PeterPetrik committed Nov 10, 2020
1 parent b2f49f0 commit da3c312
Show file tree
Hide file tree
Showing 5 changed files with 322 additions and 13 deletions.
214 changes: 203 additions & 11 deletions src/analysis/processing/qgsalgorithmexportmesh.cpp
View file
Expand Up @@ -21,7 +21,9 @@
#include "qgsmeshlayer.h"
#include "qgsmeshlayerutils.h"
#include "qgsmeshlayertemporalproperties.h"
#include "qgsmeshlayerinterpolator.h"
#include "qgspolygon.h"
#include "qgsrasterfilewriter.h"
#include "qgslinestring.h"

///@cond PRIVATE
Expand Down Expand Up @@ -390,27 +392,27 @@ QgsGeometry QgsExportMeshEdgesAlgorithm::meshElement( int index ) const
}


QString QgsExportMeshOnGrid::name() const {return QStringLiteral( "exportmeshongrid" );}
QString QgsExportMeshOnGridAlgorithm::name() const {return QStringLiteral( "exportmeshongrid" );}

QString QgsExportMeshOnGrid::displayName() const {return QStringLiteral( "Export mesh on grid" );}
QString QgsExportMeshOnGridAlgorithm::displayName() const {return QObject::tr( "Export mesh on grid" );}

QString QgsExportMeshOnGrid::group() const {return QStringLiteral( "Mesh" );}
QString QgsExportMeshOnGridAlgorithm::group() const {return QObject::tr( "Mesh" );}

QString QgsExportMeshOnGrid::groupId() const {return QStringLiteral( "mesh" );}
QString QgsExportMeshOnGridAlgorithm::groupId() const {return QStringLiteral( "mesh" );}

QString QgsExportMeshOnGrid::shortHelpString() const
QString QgsExportMeshOnGridAlgorithm::shortHelpString() const
{
return QObject::tr( "Exports mesh layer's dataset values to a gridded point vector layer, with the dataset values on this point as attribute values.\n"
"For data on volume (3D stacked dataset values), the exported dataset values are averaged on faces using the method defined in the mesh layer properties (default is Multi level averaging method).\n"
"1D meshes are not supported" );
}

QgsProcessingAlgorithm *QgsExportMeshOnGrid::createInstance() const
QgsProcessingAlgorithm *QgsExportMeshOnGridAlgorithm::createInstance() const
{
return new QgsExportMeshOnGrid();
return new QgsExportMeshOnGridAlgorithm();
}

void QgsExportMeshOnGrid::initAlgorithm( const QVariantMap &configuration )
void QgsExportMeshOnGridAlgorithm::initAlgorithm( const QVariantMap &configuration )
{
Q_UNUSED( configuration );

Expand Down Expand Up @@ -442,7 +444,7 @@ void QgsExportMeshOnGrid::initAlgorithm( const QVariantMap &configuration )
addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output vector layer" ), QgsProcessing::TypeVectorPoint ) );
}

bool QgsExportMeshOnGrid::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
bool QgsExportMeshOnGridAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, QStringLiteral( "INPUT" ), context );

Expand Down Expand Up @@ -521,7 +523,7 @@ bool QgsExportMeshOnGrid::prepareAlgorithm( const QVariantMap &parameters, QgsPr
return true;
}

QVariantMap QgsExportMeshOnGrid::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
QVariantMap QgsExportMeshOnGridAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
if ( feedback )
{
Expand Down Expand Up @@ -659,7 +661,7 @@ QVariantMap QgsExportMeshOnGrid::processAlgorithm( const QVariantMap &parameters
return ret;
}

QSet<int> QgsExportMeshOnGrid::supportedDataType()
QSet<int> QgsExportMeshOnGridAlgorithm::supportedDataType()
{
return QSet<int>(
{
Expand All @@ -669,3 +671,193 @@ QSet<int> QgsExportMeshOnGrid::supportedDataType()
}

///@endcond PRIVATE

void QgsMeshRasterizeAlgorithm::initAlgorithm( const QVariantMap &configuration )
{
Q_UNUSED( configuration );

addParameter( new QgsProcessingParameterMeshLayer( QStringLiteral( "INPUT" ), QObject::tr( "Input Mesh Layer" ) ) );

addParameter( new QgsProcessingParameterMeshDatasetGroups(
QStringLiteral( "DATASET_GROUPS" ),
QObject::tr( "Dataset groups" ),
QStringLiteral( "INPUT" ),
supportedDataType() ) );

addParameter( new QgsProcessingParameterMeshDatasetTime(
QStringLiteral( "DATASET_TIME" ),
QObject::tr( "Dataset time" ),
QStringLiteral( "INPUT" ),
QStringLiteral( "DATASET_GROUPS" ) ) );

addParameter( new QgsProcessingParameterExtent( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ), QVariant(), true ) );

addParameter( new QgsProcessingParameterDistance( QStringLiteral( "PIXEL_SIZE" ), QObject::tr( "Pixel size" ), 1, QStringLiteral( "INPUT" ), false ) );

addParameter( new QgsProcessingParameterCrs( QStringLiteral( "CRS_OUTPUT" ), QObject::tr( "Output coordinate system" ), QVariant(), true ) );

addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Output raster layer" ) ) );
}

bool QgsMeshRasterizeAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, QStringLiteral( "INPUT" ), context );

if ( !meshLayer || !meshLayer->isValid() )
return false;

QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, QStringLiteral( "CRS_OUTPUT" ), context );
if ( !outputCrs.isValid() )
outputCrs = meshLayer->crs();
mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() );
if ( !meshLayer->nativeMesh() )
meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh

mTriangularMesh = *meshLayer->triangularMesh();
QgsMesh *nativeMesh = meshLayer->nativeMesh();

QList<int> datasetGroups =
QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( QStringLiteral( "DATASET_GROUPS" ) ) );

if ( feedback )
{
feedback->setProgressText( QObject::tr( "Preparing data" ) );
}

QDateTime layerReferenceTime = static_cast<QgsMeshLayerTemporalProperties *>( meshLayer->temporalProperties() )->referenceTime();

//! Extract the date time used to export dataset values under a relative time
QgsInterval relativeTime( 0 );
QVariant parameterTimeVariant = parameters.value( QStringLiteral( "DATASET_TIME" ) );

QString timeType = QgsProcessingParameterMeshDatasetTime::valueAsTimeType( parameterTimeVariant );

if ( timeType == QStringLiteral( "dataset-time-step" ) )
{
QgsMeshDatasetIndex datasetIndex = QgsProcessingParameterMeshDatasetTime::timeValueAsDatasetIndex( parameterTimeVariant );
relativeTime = meshLayer->datasetRelativeTime( datasetIndex );
}
else if ( timeType == QStringLiteral( "defined-date-time" ) )
{
QDateTime dateTime = QgsProcessingParameterMeshDatasetTime::timeValueAsDefinedDateTime( parameterTimeVariant );
if ( dateTime.isValid() )
relativeTime = QgsInterval( layerReferenceTime.secsTo( dateTime ) );
}
else if ( timeType == QStringLiteral( "current-context-time" ) )
{
QDateTime dateTime = context.currentTimeRange().begin();
if ( dateTime.isValid() )
relativeTime = QgsInterval( layerReferenceTime.secsTo( dateTime ) );
}

for ( int i = 0; i < datasetGroups.count(); ++i )
{
int groupIndex = datasetGroups.at( i );
QgsMeshDatasetIndex datasetIndex = meshLayer->datasetIndexAtRelativeTime( relativeTime, groupIndex );

DataGroup dataGroup;
dataGroup.metadata = meshLayer->datasetGroupMetadata( datasetIndex );
if ( supportedDataType().contains( dataGroup.metadata.dataType() ) )
{
int valueCount = dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices ?
mTriangularMesh.vertices().count() : nativeMesh->faceCount();
dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, valueCount );
dataGroup.activeFaces = meshLayer->areFacesActive( datasetIndex, 0, nativeMesh->faceCount() );
if ( dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVolumes )
{
dataGroup.dataset3dStakedValue = meshLayer->dataset3dValues( datasetIndex, 0, valueCount );
}
mDataPerGroup.append( dataGroup );
}
if ( feedback )
feedback->setProgress( 100 * i / datasetGroups.count() );
}

return true;
}

QVariantMap QgsMeshRasterizeAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
if ( feedback )
{
if ( feedback->isCanceled() )
return QVariantMap();
feedback->setProgress( 0 );
feedback->setProgressText( QObject::tr( "Creating raster layer" ) );
}

//First, if present, average 3D staked dataset value to 2D face value
const QgsMesh3dAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
for ( DataGroup &dataGroup : mDataPerGroup )
{
if ( dataGroup.dataset3dStakedValue.isValid() )
dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
}

// create raster
double pixelSize = parameterAsDouble( parameters, QStringLiteral( "PIXEL_SIZE" ), context );
QgsRectangle extent = parameterAsExtent( parameters, QStringLiteral( "EXTENT" ), context );

int width = extent.width() / pixelSize;
int height = extent.height() / pixelSize;

QString fileName = parameterAsFileOutput( parameters, QStringLiteral( "OUTPUT" ), context );
QFileInfo fileInfo( fileName );
QString outputFormat = QgsRasterFileWriter::driverForExtension( fileInfo.suffix() );
QgsRasterFileWriter rasterFileWriter( fileName );
rasterFileWriter.setOutputProviderKey( QStringLiteral( "gdal" ) );
rasterFileWriter.setOutputFormat( outputFormat );

std::unique_ptr<QgsRasterDataProvider> rasterDataProvider(
rasterFileWriter.createMultiBandRaster( Qgis::Float64, width, height, extent, mTransform.destinationCrs(), mDataPerGroup.count() ) );
rasterDataProvider->setEditable( true );

for ( int i = 0; i < mDataPerGroup.count(); ++i )
{
const DataGroup &dataGroup = mDataPerGroup.at( i );
QgsRasterBlockFeedback rasterBlockFeedBack;
if ( feedback )
QObject::connect( &rasterBlockFeedBack, &QgsFeedback::canceled, feedback, &QgsFeedback::cancel );

QgsRasterBlock *block = QgsMeshUtils::exportRasterBlock(
mTriangularMesh,
dataGroup.datasetValues,
dataGroup.activeFaces,
dataGroup.metadata.dataType(),
mTransform,
pixelSize,
extent,
&rasterBlockFeedBack );



rasterDataProvider->writeBlock( block, i + 1 );
rasterDataProvider->setNoDataValue( i + 1, block->noDataValue() );

if ( feedback )
{
if ( feedback->isCanceled() )
return QVariantMap();
feedback->setProgress( 100 * i / mDataPerGroup.count() );
}
}

rasterDataProvider->setEditable( false );

if ( feedback )
feedback->setProgress( 100 );

QVariantMap ret;
ret[QStringLiteral( "OUTPUT" )] = fileName;

return ret;
}

QSet<int> QgsMeshRasterizeAlgorithm::supportedDataType()
{
return QSet<int>(
{
QgsMeshDatasetGroupMetadata::DataOnVertices,
QgsMeshDatasetGroupMetadata::DataOnFaces,
QgsMeshDatasetGroupMetadata::DataOnVolumes} );
}
54 changes: 53 additions & 1 deletion src/analysis/processing/qgsalgorithmexportmesh.h
View file
Expand Up @@ -125,7 +125,7 @@ class QgsExportMeshEdgesAlgorithm : public QgsExportMeshOnElement
};


class QgsExportMeshOnGrid : public QgsProcessingAlgorithm
class QgsExportMeshOnGridAlgorithm : public QgsProcessingAlgorithm
{

public:
Expand Down Expand Up @@ -160,6 +160,58 @@ class QgsExportMeshOnGrid : public QgsProcessingAlgorithm
QgsMeshRendererSettings mLayerRendererSettings;
};

class QgsMeshRasterizeAlgorithm : public QgsProcessingAlgorithm
{

public:
QString name() const override
{
return QStringLiteral( "meshrasterize" );
}
QString displayName() const override
{
return QObject::tr( "Rasterize mesh dataset" );
}
QString group() const override
{
return QObject::tr( "Mesh" );
}
QString groupId() const override
{
return QStringLiteral( "mesh" );
}
QString shortHelpString() const override
{
return QObject::tr( "Create a raster layer from a mesh dataset" );
}

protected:
QgsProcessingAlgorithm *createInstance() const override
{
return new QgsMeshRasterizeAlgorithm();
}
void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
QVariantMap processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;

private:

QSet<int> supportedDataType();

QgsTriangularMesh mTriangularMesh;
struct DataGroup
{
QgsMeshDatasetGroupMetadata metadata;
QgsMeshDataBlock datasetValues;
QgsMeshDataBlock activeFaces;
QgsMesh3dDataBlock dataset3dStakedValue; //will be filled only if data are 3d stacked
};

QList<DataGroup> mDataPerGroup;
QgsCoordinateTransform mTransform;
QgsMeshRendererSettings mLayerRendererSettings;
};

///@endcond PRIVATE

#endif // QGSALGORITHMEXPORTMESHVERTICES_H
3 changes: 2 additions & 1 deletion src/analysis/processing/qgsnativealgorithms.cpp
View file
Expand Up @@ -287,7 +287,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsExportMeshVerticesAlgorithm );
addAlgorithm( new QgsExportMeshFacesAlgorithm );
addAlgorithm( new QgsExportMeshEdgesAlgorithm );
addAlgorithm( new QgsExportMeshOnGrid );
addAlgorithm( new QgsExportMeshOnGridAlgorithm );
addAlgorithm( new QgsExtendLinesAlgorithm() );
addAlgorithm( new QgsExtentFromLayerAlgorithm() );
addAlgorithm( new QgsExtentToLayerAlgorithm() );
Expand Down Expand Up @@ -341,6 +341,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
addAlgorithm( new QgsMeanCoordinatesAlgorithm() );
addAlgorithm( new QgsMergeLinesAlgorithm() );
addAlgorithm( new QgsMergeVectorAlgorithm() );
addAlgorithm( new QgsMeshRasterizeAlgorithm );
addAlgorithm( new QgsMinimumEnclosingCircleAlgorithm() );
addAlgorithm( new QgsMultipartToSinglepartAlgorithm() );
addAlgorithm( new QgsMultiRingConstantBufferAlgorithm() );
Expand Down
41 changes: 41 additions & 0 deletions src/core/mesh/qgsmeshlayerinterpolator.cpp
View file
Expand Up @@ -241,3 +241,44 @@ QgsRasterBlock *QgsMeshUtils::exportRasterBlock(

return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
}

QgsRasterBlock *QgsMeshUtils::exportRasterBlock(
const QgsTriangularMesh &triangularMesh,
const QgsMeshDataBlock &datasetValues,
const QgsMeshDataBlock &activeFlags,
const QgsMeshDatasetGroupMetadata::DataType dataType,
const QgsCoordinateTransform &transform,
double mapUnitsPerPixel,
const QgsRectangle &extent,
QgsRasterBlockFeedback *feedback )
{

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 );

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

QVector<double> magnitudes = QgsMeshLayerUtils::calculateMagnitudes( datasetValues );

QgsMeshLayerInterpolator interpolator(
triangularMesh,
magnitudes,
activeFlags,
dataType,
renderContext,
QSize( widthPixel, heightPixel )
);

return interpolator.block( 0, extent, widthPixel, heightPixel, feedback );
}

0 comments on commit da3c312

Please sign in to comment.