QGIS API Documentation 4.1.0-Master (ca2ac17535b)
Loading...
Searching...
No Matches
qgsalgorithmexportmesh.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmexportmesh.cpp
3 ---------------------------
4 begin : October 2020
5 copyright : (C) 2020 by Vincent Cloarec
6 email : vcloarec at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19
20#include "qgslinestring.h"
21#include "qgsmeshcontours.h"
22#include "qgsmeshdataset.h"
23#include "qgsmeshlayer.h"
26#include "qgsmeshlayerutils.h"
27#include "qgspolygon.h"
29#include "qgsrasterfilewriter.h"
30
31#include <QString>
32#include <QTextStream>
33
34using namespace Qt::StringLiterals;
35
37
38
39static QgsFields createFields( const QList<QgsMeshDatasetGroupMetadata> &groupMetadataList, int vectorOption )
40{
41 QgsFields fields;
42 for ( const QgsMeshDatasetGroupMetadata &meta : groupMetadataList )
43 {
44 if ( meta.isVector() )
45 {
46 if ( vectorOption == 0 || vectorOption == 2 )
47 {
48 fields.append( QgsField( u"%1_x"_s.arg( meta.name() ), QMetaType::Type::Double ) );
49 fields.append( QgsField( u"%1_y"_s.arg( meta.name() ), QMetaType::Type::Double ) );
50 }
51
52 if ( vectorOption == 1 || vectorOption == 2 )
53 {
54 fields.append( QgsField( u"%1_mag"_s.arg( meta.name() ), QMetaType::Type::Double ) );
55 fields.append( QgsField( u"%1_dir"_s.arg( meta.name() ), QMetaType::Type::Double ) );
56 }
57 }
58 else
59 fields.append( QgsField( meta.name(), QMetaType::Type::Double ) );
60 }
61 return fields;
62}
63
64static QVector<double> vectorValue( const QgsMeshDatasetValue &value, int exportOption )
65{
66 QVector<double> ret( exportOption == 2 ? 4 : 2 );
67
68 if ( exportOption == 0 || exportOption == 2 )
69 {
70 ret[0] = value.x();
71 ret[1] = value.y();
72 }
73 if ( exportOption == 1 || exportOption == 2 )
74 {
75 double x = value.x();
76 double y = value.y();
77 double magnitude = sqrt( x * x + y * y );
78 double direction = ( asin( x / magnitude ) ) / M_PI * 180;
79 if ( y < 0 )
80 direction = 180 - direction;
81
82 if ( exportOption == 1 )
83 {
84 ret[0] = magnitude;
85 ret[1] = direction;
86 }
87 if ( exportOption == 2 )
88 {
89 ret[2] = magnitude;
90 ret[3] = direction;
91 }
92 }
93 return ret;
94}
95
96static void addAttributes( const QgsMeshDatasetValue &value, QgsAttributes &attributes, bool isVector, int vectorOption )
97{
98 if ( isVector )
99 {
100 QVector<double> vectorValues = vectorValue( value, vectorOption );
101 for ( double v : vectorValues )
102 {
103 if ( v == std::numeric_limits<double>::quiet_NaN() )
104 attributes.append( QVariant() );
105 else
106 attributes.append( v );
107 }
108 }
109 else
110 {
111 if ( value.scalar() == std::numeric_limits<double>::quiet_NaN() )
112 attributes.append( QVariant() );
113 else
114 attributes.append( value.scalar() );
115 }
116}
117
118static QgsMeshDatasetValue extractDatasetValue(
119 const QgsPointXY &point,
120 int nativeFaceIndex,
121 int triangularFaceIndex,
122 const QgsTriangularMesh &triangularMesh,
123 const QgsMeshDataBlock &activeFaces,
124 const QgsMeshDataBlock &datasetValues,
125 const QgsMeshDatasetGroupMetadata &metadata
126)
127{
128 bool faceActive = activeFaces.active( nativeFaceIndex );
130 if ( faceActive )
131 {
132 switch ( metadata.dataType() )
133 {
135 //not supported
136 break;
139 {
140 value = datasetValues.value( nativeFaceIndex );
141 }
142 break;
143
145 {
146 const QgsMeshFace &face = triangularMesh.triangles()[triangularFaceIndex];
147 const int v1 = face[0], v2 = face[1], v3 = face[2];
148 const QgsPoint p1 = triangularMesh.vertices()[v1], p2 = triangularMesh.vertices()[v2], p3 = triangularMesh.vertices()[v3];
149 const QgsMeshDatasetValue val1 = datasetValues.value( v1 );
150 const QgsMeshDatasetValue val2 = datasetValues.value( v2 );
151 const QgsMeshDatasetValue val3 = datasetValues.value( v3 );
152 const double x = QgsMeshLayerUtils::interpolateFromVerticesData( p1, p2, p3, val1.x(), val2.x(), val3.x(), point );
153 double y = std::numeric_limits<double>::quiet_NaN();
154 bool isVector = metadata.isVector();
155 if ( isVector )
156 y = QgsMeshLayerUtils::interpolateFromVerticesData( p1, p2, p3, val1.y(), val2.y(), val3.y(), point );
157
158 value = QgsMeshDatasetValue( x, y );
159 }
160 break;
161 }
162 }
163
164 return value;
165}
166
167QString QgsExportMeshOnElement::group() const
168{
169 return QObject::tr( "Mesh" );
170}
171
172QString QgsExportMeshOnElement::groupId() const
173{
174 return u"mesh"_s;
175}
176
177QString QgsExportMeshVerticesAlgorithm::shortHelpString() const
178{
179 return QObject::tr( "This algorithm exports a mesh layer's vertices to a point vector layer, with the dataset values on vertices as attribute values." );
180}
181
182QString QgsExportMeshVerticesAlgorithm::shortDescription() const
183{
184 return QObject::tr( "Exports mesh vertices to a point vector layer." );
185}
186
187QString QgsExportMeshVerticesAlgorithm::name() const
188{
189 return u"exportmeshvertices"_s;
190}
191
192QString QgsExportMeshVerticesAlgorithm::displayName() const
193{
194 return QObject::tr( "Export mesh vertices" );
195}
196
197QgsProcessingAlgorithm *QgsExportMeshVerticesAlgorithm::createInstance() const
198{
199 return new QgsExportMeshVerticesAlgorithm();
200}
201
202QgsGeometry QgsExportMeshVerticesAlgorithm::meshElement( int index ) const
203{
204 return QgsGeometry( new QgsPoint( mNativeMesh.vertex( index ) ) );
205}
206
207void QgsExportMeshOnElement::initAlgorithm( const QVariantMap &configuration )
208{
209 Q_UNUSED( configuration );
210
211 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
212
213
214 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType(), true ) );
215
216 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
217
218 addParameter( new QgsProcessingParameterCrs( u"CRS_OUTPUT"_s, QObject::tr( "Output coordinate system" ), QVariant(), true ) );
219
220 QStringList exportVectorOptions;
221 exportVectorOptions << QObject::tr( "Cartesian (x,y)" ) << QObject::tr( "Polar (magnitude,degree)" ) << QObject::tr( "Cartesian and Polar" );
222 addParameter( new QgsProcessingParameterEnum( u"VECTOR_OPTION"_s, QObject::tr( "Export vector option" ), exportVectorOptions, false, 0 ) );
223 addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT"_s, QObject::tr( "Output vector layer" ), sinkType() ) );
224}
225
226static QgsInterval datasetRelativetime( const QVariant parameterTimeVariant, QgsMeshLayer *meshLayer, const QgsProcessingContext &context )
227{
228 QgsInterval relativeTime( 0 );
229 QDateTime layerReferenceTime = static_cast<QgsMeshLayerTemporalProperties *>( meshLayer->temporalProperties() )->referenceTime();
230 QString timeType = QgsProcessingParameterMeshDatasetTime::valueAsTimeType( parameterTimeVariant );
231
232 if ( timeType == "dataset-time-step"_L1 )
233 {
235 relativeTime = meshLayer->datasetRelativeTime( datasetIndex );
236 }
237 else if ( timeType == "defined-date-time"_L1 )
238 {
239 QDateTime dateTime = QgsProcessingParameterMeshDatasetTime::timeValueAsDefinedDateTime( parameterTimeVariant );
240 if ( dateTime.isValid() )
241 relativeTime = QgsInterval( layerReferenceTime.secsTo( dateTime ) );
242 }
243 else if ( timeType == "current-context-time"_L1 )
244 {
245 QDateTime dateTime = context.currentTimeRange().begin();
246 if ( dateTime.isValid() )
247 relativeTime = QgsInterval( layerReferenceTime.secsTo( dateTime ) );
248 }
249
250 return relativeTime;
251}
252
253
254bool QgsExportMeshOnElement::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
255{
256 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
257
258 if ( !meshLayer || !meshLayer->isValid() )
259 return false;
260
261 if ( meshLayer->isEditable() )
262 throw QgsProcessingException( QObject::tr( "Input mesh layer in edit mode is not supported" ) );
263
264 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
265 if ( !outputCrs.isValid() )
266 outputCrs = meshLayer->crs();
267 mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() );
268 if ( !meshLayer->nativeMesh() )
269 meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh
270
271 mNativeMesh = *meshLayer->nativeMesh();
272
273 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
274
275 if ( feedback )
276 {
277 feedback->setProgressText( QObject::tr( "Preparing data" ) );
278 }
279
280 // Extract the date time used to export dataset values under a relative time
281 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
282 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
283
284 switch ( meshElementType() )
285 {
286 case QgsMesh::Face:
287 mElementCount = mNativeMesh.faceCount();
288 break;
289 case QgsMesh::Vertex:
290 mElementCount = mNativeMesh.vertexCount();
291 break;
292 case QgsMesh::Edge:
293 mElementCount = mNativeMesh.edgeCount();
294 break;
295 }
296
297 for ( int i = 0; i < datasetGroups.count(); ++i )
298 {
299 int groupIndex = datasetGroups.at( i );
300 QgsMeshDatasetIndex datasetIndex = meshLayer->datasetIndexAtRelativeTime( relativeTime, groupIndex );
301
302 DataGroup dataGroup;
303 dataGroup.metadata = meshLayer->datasetGroupMetadata( datasetIndex );
304 if ( supportedDataType().contains( dataGroup.metadata.dataType() ) )
305 {
306 dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, mElementCount );
307 mDataPerGroup.append( dataGroup );
308 }
309 if ( feedback )
310 feedback->setProgress( 100 * i / datasetGroups.count() );
311 }
312
313 mExportVectorOption = parameterAsInt( parameters, u"VECTOR_OPTION"_s, context );
314
315 return true;
316}
317
318QVariantMap QgsExportMeshOnElement::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
319{
320 if ( feedback )
321 {
322 if ( feedback->isCanceled() )
323 return QVariantMap();
324 feedback->setProgress( 0 );
325 feedback->setProgressText( QObject::tr( "Creating output vector layer" ) );
326 }
327
328 QList<QgsMeshDatasetGroupMetadata> metaList;
329 metaList.reserve( mDataPerGroup.size() );
330 for ( const DataGroup &dataGroup : std::as_const( mDataPerGroup ) )
331 metaList.append( dataGroup.metadata );
332 QgsFields fields = createFields( metaList, mExportVectorOption );
333
334 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
335 QString identifier;
336 std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, u"OUTPUT"_s, context, identifier, fields, sinkGeometryType(), outputCrs ) );
337 if ( !sink )
338 return QVariantMap();
339
340 if ( feedback )
341 {
342 if ( feedback->isCanceled() )
343 return QVariantMap();
344 feedback->setProgress( 0 );
345 feedback->setProgressText( QObject::tr( "Creating points for each vertices" ) );
346 }
347
348 for ( int i = 0; i < mElementCount; ++i )
349 {
350 QgsAttributes attributes;
351 for ( const DataGroup &dataGroup : std::as_const( mDataPerGroup ) )
352 {
353 const QgsMeshDatasetValue &value = dataGroup.datasetValues.value( i );
354 addAttributes( value, attributes, dataGroup.metadata.isVector(), mExportVectorOption );
355 }
356
357 QgsFeature feat;
358 QgsGeometry geom = meshElement( i );
359 try
360 {
361 geom.transform( mTransform );
362 }
363 catch ( QgsCsException & )
364 {
365 geom = meshElement( i );
366 if ( feedback )
367 feedback->reportError( QObject::tr( "Could not transform point to destination CRS" ) );
368 }
369 feat.setGeometry( geom );
370 feat.setAttributes( attributes );
371
372 if ( !sink->addFeature( feat, QgsFeatureSink::FastInsert ) )
373 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, u"OUTPUT"_s ) );
374 else
375 feedback->featureAddedToSink( u"OUTPUT"_s );
376
377 if ( feedback )
378 {
379 if ( feedback->isCanceled() )
380 return QVariantMap();
381 feedback->setProgress( 100 * i / mElementCount );
382 }
383 }
384
385 sink->finalize();
386 feedback->featureSinkFinalized( u"OUTPUT"_s );
387
388 QVariantMap ret;
389 ret[u"OUTPUT"_s] = identifier;
390
391 return ret;
392}
393
394QString QgsExportMeshFacesAlgorithm::shortHelpString() const
395{
396 return QObject::tr( "This algorithm exports a mesh layer's faces to a polygon vector layer, with the dataset values on faces as attribute values." );
397}
398
399QString QgsExportMeshFacesAlgorithm::shortDescription() const
400{
401 return QObject::tr( "Exports mesh faces to a polygon vector layer." );
402}
403
404QString QgsExportMeshFacesAlgorithm::name() const
405{
406 return u"exportmeshfaces"_s;
407}
408
409QString QgsExportMeshFacesAlgorithm::displayName() const
410{
411 return QObject::tr( "Export mesh faces" );
412}
413
414QgsProcessingAlgorithm *QgsExportMeshFacesAlgorithm::createInstance() const
415{
416 return new QgsExportMeshFacesAlgorithm();
417}
418
419QgsGeometry QgsExportMeshFacesAlgorithm::meshElement( int index ) const
420{
421 const QgsMeshFace &face = mNativeMesh.face( index );
422 QVector<QgsPoint> vertices( face.size() );
423 for ( int i = 0; i < face.size(); ++i )
424 vertices[i] = mNativeMesh.vertex( face.at( i ) );
425 auto polygon = std::make_unique<QgsPolygon>();
426 polygon->setExteriorRing( new QgsLineString( vertices ) );
427 return QgsGeometry( polygon.release() );
428}
429
430QString QgsExportMeshEdgesAlgorithm::shortHelpString() const
431{
432 return QObject::tr( "This algorithm exports a mesh layer's edges to a line vector layer, with the dataset values on edges as attribute values." );
433}
434
435QString QgsExportMeshEdgesAlgorithm::shortDescription() const
436{
437 return QObject::tr( "Exports mesh edges to a line vector layer." );
438}
439
440QString QgsExportMeshEdgesAlgorithm::name() const
441{
442 return u"exportmeshedges"_s;
443}
444
445QString QgsExportMeshEdgesAlgorithm::displayName() const
446{
447 return QObject::tr( "Export mesh edges" );
448}
449
450QgsProcessingAlgorithm *QgsExportMeshEdgesAlgorithm::createInstance() const
451{
452 return new QgsExportMeshEdgesAlgorithm();
453}
454
455QgsGeometry QgsExportMeshEdgesAlgorithm::meshElement( int index ) const
456{
457 const QgsMeshEdge &edge = mNativeMesh.edge( index );
458 QVector<QgsPoint> vertices( 2 );
459 vertices[0] = mNativeMesh.vertex( edge.first );
460 vertices[1] = mNativeMesh.vertex( edge.second );
461 return QgsGeometry( new QgsLineString( vertices ) );
462}
463
464
465QString QgsExportMeshOnGridAlgorithm::name() const
466{
467 return u"exportmeshongrid"_s;
468}
469
470QString QgsExportMeshOnGridAlgorithm::displayName() const
471{
472 return QObject::tr( "Export mesh on grid" );
473}
474
475QString QgsExportMeshOnGridAlgorithm::group() const
476{
477 return QObject::tr( "Mesh" );
478}
479
480QString QgsExportMeshOnGridAlgorithm::groupId() const
481{
482 return u"mesh"_s;
483}
484
485QString QgsExportMeshOnGridAlgorithm::shortHelpString() const
486{
487 return QObject::tr(
488 "This algorithm exports a mesh layer's dataset values to a gridded point vector layer, with the dataset values on each point as attribute values.\n"
489 "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 "
490 "method).\n"
491 "1D meshes are not supported."
492 );
493}
494
495QString QgsExportMeshOnGridAlgorithm::shortDescription() const
496{
497 return QObject::tr( "Exports mesh dataset values to a gridded point vector layer." );
498}
499
500QgsProcessingAlgorithm *QgsExportMeshOnGridAlgorithm::createInstance() const
501{
502 return new QgsExportMeshOnGridAlgorithm();
503}
504
505void QgsExportMeshOnGridAlgorithm::initAlgorithm( const QVariantMap &configuration )
506{
507 Q_UNUSED( configuration );
508
509 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
510
511 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
512
513 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
514
515 addParameter( new QgsProcessingParameterExtent( u"EXTENT"_s, QObject::tr( "Extent" ), QVariant(), true ) );
516
517 addParameter( new QgsProcessingParameterDistance( u"GRID_SPACING"_s, QObject::tr( "Grid spacing" ), 10, u"INPUT"_s, false ) );
518
519 addParameter( new QgsProcessingParameterCrs( u"CRS_OUTPUT"_s, QObject::tr( "Output coordinate system" ), QVariant(), true ) );
520
521 QStringList exportVectorOptions;
522 exportVectorOptions << QObject::tr( "Cartesian (x,y)" ) << QObject::tr( "Polar (magnitude,degree)" ) << QObject::tr( "Cartesian and Polar" );
523 addParameter( new QgsProcessingParameterEnum( u"VECTOR_OPTION"_s, QObject::tr( "Export vector option" ), exportVectorOptions, false, 0 ) );
524 addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT"_s, QObject::tr( "Output vector layer" ), Qgis::ProcessingSourceType::VectorPoint ) );
525}
526
527static void extractDatasetValues(
528 const QList<int> &datasetGroups,
529 QgsMeshLayer *meshLayer,
530 const QgsMesh &nativeMesh,
531 const QgsInterval &relativeTime,
532 const QSet<int> supportedDataType,
533 QList<DataGroup> &datasetPerGroup,
534 QgsProcessingFeedback *feedback
535)
536{
537 for ( int i = 0; i < datasetGroups.count(); ++i )
538 {
539 int groupIndex = datasetGroups.at( i );
540 QgsMeshDatasetIndex datasetIndex = meshLayer->datasetIndexAtRelativeTime( relativeTime, groupIndex );
541
542 DataGroup dataGroup;
543 dataGroup.metadata = meshLayer->datasetGroupMetadata( datasetIndex );
544 if ( supportedDataType.contains( dataGroup.metadata.dataType() ) )
545 {
546 int valueCount = dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices ? nativeMesh.vertices.count() : nativeMesh.faceCount();
547 dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, valueCount );
548 dataGroup.activeFaces = meshLayer->areFacesActive( datasetIndex, 0, nativeMesh.faceCount() );
549 if ( dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVolumes )
550 {
551 dataGroup.dataset3dStakedValue = meshLayer->dataset3dValues( datasetIndex, 0, valueCount );
552 }
553 datasetPerGroup.append( dataGroup );
554 }
555 if ( feedback )
556 feedback->setProgress( 100 * i / datasetGroups.count() );
557 }
558}
559
560bool QgsExportMeshOnGridAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
561{
562 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
563
564 if ( !meshLayer || !meshLayer->isValid() )
565 return false;
566
567 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
568 if ( !outputCrs.isValid() )
569 outputCrs = meshLayer->crs();
570 mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() );
571 if ( !meshLayer->nativeMesh() )
572 meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh
573
574 const QgsMesh &nativeMesh = *meshLayer->nativeMesh();
575
576 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
577
578 if ( feedback )
579 {
580 feedback->setProgressText( QObject::tr( "Preparing data" ) );
581 }
582
583 // Extract the date time used to export dataset values under a relative time
584 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
585 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
586
587 extractDatasetValues( datasetGroups, meshLayer, nativeMesh, relativeTime, supportedDataType(), mDataPerGroup, feedback );
588 mTriangularMesh.update( meshLayer->nativeMesh(), mTransform );
589
590 mExportVectorOption = parameterAsInt( parameters, u"VECTOR_OPTION"_s, context );
591
592 return true;
593}
594
595QVariantMap QgsExportMeshOnGridAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
596{
597 if ( feedback )
598 {
599 if ( feedback->isCanceled() )
600 return QVariantMap();
601 feedback->setProgress( 0 );
602 feedback->setProgressText( QObject::tr( "Creating output vector layer" ) );
603 }
604
605 //First, if present, average 3D staked dataset value to 2D face value
606 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
607 for ( DataGroup &dataGroup : mDataPerGroup )
608 {
609 if ( dataGroup.dataset3dStakedValue.isValid() )
610 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
611 }
612
613 QList<QgsMeshDatasetGroupMetadata> metaList;
614 metaList.reserve( mDataPerGroup.size() );
615 for ( const DataGroup &dataGroup : std::as_const( mDataPerGroup ) )
616 metaList.append( dataGroup.metadata );
617 QgsFields fields = createFields( metaList, mExportVectorOption );
618
619 //create sink
620 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
621 QString identifier;
622 std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, u"OUTPUT"_s, context, identifier, fields, Qgis::WkbType::Point, outputCrs ) );
623 if ( !sink )
624 return QVariantMap();
625
626 if ( feedback )
627 {
628 if ( feedback->isCanceled() )
629 return QVariantMap();
630 feedback->setProgress( 0 );
631 feedback->setProgressText( QObject::tr( "Creating gridded points" ) );
632 }
633
634 // grid definition
635 const double gridSpacing = parameterAsDouble( parameters, u"GRID_SPACING"_s, context );
636 if ( qgsDoubleNear( gridSpacing, 0 ) )
637 {
638 throw QgsProcessingException( QObject::tr( "Grid spacing cannot be 0" ) );
639 }
640
641 QgsRectangle extent = parameterAsExtent( parameters, u"EXTENT"_s, context );
642 if ( extent.isEmpty() )
643 extent = mTriangularMesh.extent();
644 int pointXCount = int( extent.width() / gridSpacing ) + 1;
645 int pointYCount = int( extent.height() / gridSpacing ) + 1;
646
647 for ( int ix = 0; ix < pointXCount; ++ix )
648 {
649 for ( int iy = 0; iy < pointYCount; ++iy )
650 {
651 QgsPoint point( extent.xMinimum() + ix * gridSpacing, extent.yMinimum() + iy * gridSpacing );
652 int triangularFaceIndex = mTriangularMesh.faceIndexForPoint_v2( point );
653 if ( triangularFaceIndex >= 0 )
654 {
655 //extract dataset values for the point
656 QgsAttributes attributes;
657 int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces().at( triangularFaceIndex );
658 for ( int i = 0; i < mDataPerGroup.count(); ++i )
659 {
660 const DataGroup &dataGroup = mDataPerGroup.at( i );
661 bool faceActive = dataGroup.activeFaces.active( nativeFaceIndex );
662 if ( !faceActive )
663 continue;
664 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
665
666 if ( dataGroup.metadata.isVector() )
667 {
668 QVector<double> vector = vectorValue( dataGroup.datasetValues.value( i ), mExportVectorOption );
669 for ( double v : vector )
670 {
671 attributes.append( v );
672 }
673 }
674 else
675 attributes.append( value.scalar() );
676 }
677 QgsFeature feat;
678 QgsGeometry geom( point.clone() );
679 try
680 {
681 geom.transform( mTransform );
682 }
683 catch ( QgsCsException & )
684 {
685 geom = QgsGeometry( point.clone() );
686 feedback->reportError( QObject::tr( "Could not transform point to destination CRS" ) );
687 }
688 feat.setGeometry( geom );
689 feat.setAttributes( attributes );
690
691 if ( !sink->addFeature( feat, QgsFeatureSink::FastInsert ) )
692 {
693 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QString() ) );
694 }
695 else
696 {
697 feedback->featureAddedToSink( u"OUTPUT"_s );
698 }
699 }
700 }
701 }
702
703 sink->finalize();
704 feedback->featureSinkFinalized( u"OUTPUT"_s );
705
706 QVariantMap ret;
707 ret[u"OUTPUT"_s] = identifier;
708
709 return ret;
710}
711
712QSet<int> QgsExportMeshOnGridAlgorithm::supportedDataType()
713{
715}
716
717QString QgsMeshRasterizeAlgorithm::name() const
718{
719 return u"meshrasterize"_s;
720}
721
722QString QgsMeshRasterizeAlgorithm::displayName() const
723{
724 return QObject::tr( "Rasterize mesh dataset" );
725}
726
727QString QgsMeshRasterizeAlgorithm::group() const
728{
729 return QObject::tr( "Mesh" );
730}
731
732QString QgsMeshRasterizeAlgorithm::groupId() const
733{
734 return u"mesh"_s;
735}
736
737QString QgsMeshRasterizeAlgorithm::shortHelpString() const
738{
739 return QObject::tr(
740 "This algorithm creates a raster layer from a mesh dataset.\n"
741 "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 "
742 "method).\n"
743 "1D meshes are not supported."
744 );
745}
746
747QString QgsMeshRasterizeAlgorithm::shortDescription() const
748{
749 return QObject::tr( "Creates a raster layer from a mesh dataset." );
750}
751
752QgsProcessingAlgorithm *QgsMeshRasterizeAlgorithm::createInstance() const
753{
754 return new QgsMeshRasterizeAlgorithm();
755}
756
757void QgsMeshRasterizeAlgorithm::initAlgorithm( const QVariantMap &configuration )
758{
759 Q_UNUSED( configuration );
760
761 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
762
763 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType(), true ) );
764
765 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
766
767 addParameter( new QgsProcessingParameterExtent( u"EXTENT"_s, QObject::tr( "Extent" ), QVariant(), true ) );
768 addParameter( new QgsProcessingParameterDistance( u"PIXEL_SIZE"_s, QObject::tr( "Pixel size" ), 1, u"INPUT"_s, false ) );
769 addParameter( new QgsProcessingParameterCrs( u"CRS_OUTPUT"_s, QObject::tr( "Output coordinate system" ), QVariant(), true ) );
770
771 // backwards compatibility parameter
772 // TODO QGIS 5: remove parameter and related logic
773 auto createOptsParam = std::make_unique<QgsProcessingParameterString>( u"CREATE_OPTIONS"_s, QObject::tr( "Creation options" ), QVariant(), false, true );
774 createOptsParam->setMetadata( QVariantMap( { { u"widget_wrapper"_s, QVariantMap( { { u"widget_type"_s, u"rasteroptions"_s } } ) } } ) );
775 createOptsParam->setFlags( createOptsParam->flags() | Qgis::ProcessingParameterFlag::Hidden );
776 addParameter( createOptsParam.release() );
777
778 auto creationOptsParam = std::make_unique<QgsProcessingParameterString>( u"CREATION_OPTIONS"_s, QObject::tr( "Creation options" ), QVariant(), false, true );
779 creationOptsParam->setMetadata( QVariantMap( { { u"widget_wrapper"_s, QVariantMap( { { u"widget_type"_s, u"rasteroptions"_s } } ) } } ) );
780 creationOptsParam->setFlags( creationOptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
781 addParameter( creationOptsParam.release() );
782
783 addParameter( new QgsProcessingParameterRasterDestination( u"OUTPUT"_s, QObject::tr( "Output raster layer" ) ) );
784}
785
786bool QgsMeshRasterizeAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
787{
788 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
789
790 if ( !meshLayer || !meshLayer->isValid() )
791 return false;
792
793 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
794 if ( !outputCrs.isValid() )
795 outputCrs = meshLayer->crs();
796 mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() );
797 if ( !meshLayer->nativeMesh() )
798 meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh
799
800 mTriangularMesh.update( meshLayer->nativeMesh(), mTransform );
801
802 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
803
804 if ( feedback )
805 {
806 feedback->setProgressText( QObject::tr( "Preparing data" ) );
807 }
808
809 // Extract the date time used to export dataset values under a relative time
810 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
811 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
812
813 extractDatasetValues( datasetGroups, meshLayer, *meshLayer->nativeMesh(), relativeTime, supportedDataType(), mDataPerGroup, feedback );
814
815 mLayerRendererSettings = meshLayer->rendererSettings();
816
817 return true;
818}
819
820QVariantMap QgsMeshRasterizeAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
821{
822 if ( feedback )
823 {
824 if ( feedback->isCanceled() )
825 return QVariantMap();
826 feedback->setProgress( 0 );
827 feedback->setProgressText( QObject::tr( "Creating raster layer" ) );
828 }
829
830 //First, if present, average 3D staked dataset value to 2D face value
831 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
832 for ( DataGroup &dataGroup : mDataPerGroup )
833 {
834 if ( dataGroup.dataset3dStakedValue.isValid() )
835 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
836 }
837
838 // create raster
839 const double pixelSize = parameterAsDouble( parameters, u"PIXEL_SIZE"_s, context );
840 if ( qgsDoubleNear( pixelSize, 0 ) )
841 {
842 throw QgsProcessingException( QObject::tr( "Pixel size cannot be 0" ) );
843 }
844
845 QgsRectangle extent = parameterAsExtent( parameters, u"EXTENT"_s, context );
846 if ( extent.isEmpty() )
847 extent = mTriangularMesh.extent();
848
849 int width = extent.width() / pixelSize;
850 int height = extent.height() / pixelSize;
851
852 QString creationOptions = parameterAsString( parameters, u"CREATION_OPTIONS"_s, context ).trimmed();
853 // handle backwards compatibility parameter CREATE_OPTIONS
854 const QString optionsString = parameterAsString( parameters, u"CREATE_OPTIONS"_s, context );
855 if ( !optionsString.isEmpty() )
856 creationOptions = optionsString;
857
858 const QString fileName = parameterAsOutputLayer( parameters, u"OUTPUT"_s, context );
859 const QString outputFormat = parameterAsOutputRasterFormat( parameters, u"OUTPUT"_s, context );
860 QgsRasterFileWriter rasterFileWriter( fileName );
861 rasterFileWriter.setOutputProviderKey( u"gdal"_s );
862 if ( !creationOptions.isEmpty() )
863 {
864 rasterFileWriter.setCreationOptions( creationOptions.split( '|' ) );
865 }
866 rasterFileWriter.setOutputFormat( outputFormat );
867
868 std::unique_ptr<QgsRasterDataProvider> rasterDataProvider( rasterFileWriter.createMultiBandRaster( Qgis::DataType::Float64, width, height, extent, mTransform.destinationCrs(), mDataPerGroup.count() ) );
869 rasterDataProvider->setEditable( true );
870
871 const bool hasReportsDuringClose = rasterDataProvider->hasReportsDuringClose();
872 const double maxProgressDuringBlockWriting = hasReportsDuringClose ? 50.0 : 100.0;
873
874 for ( int i = 0; i < mDataPerGroup.count(); ++i )
875 {
876 const DataGroup &dataGroup = mDataPerGroup.at( i );
877 QgsRasterBlockFeedback rasterBlockFeedBack;
878 if ( feedback )
879 QObject::connect( &rasterBlockFeedBack, &QgsFeedback::canceled, feedback, &QgsFeedback::cancel );
880
881 if ( dataGroup.datasetValues.isValid() )
882 {
883 std::unique_ptr<QgsRasterBlock> block(
884 QgsMeshUtils::exportRasterBlock( mTriangularMesh, dataGroup.datasetValues, dataGroup.activeFaces, dataGroup.metadata.dataType(), mTransform, pixelSize, extent, &rasterBlockFeedBack )
885 );
886
887 if ( !rasterDataProvider->writeBlock( block.get(), i + 1 ) )
888 {
889 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( rasterDataProvider->error().summary() ) );
890 }
891 rasterDataProvider->setNoDataValue( i + 1, block->noDataValue() );
892 }
893 else
894 rasterDataProvider->setNoDataValue( i + 1, std::numeric_limits<double>::quiet_NaN() );
895
896 if ( feedback )
897 {
898 if ( feedback->isCanceled() )
899 return QVariantMap();
900 feedback->setProgress( maxProgressDuringBlockWriting * i / mDataPerGroup.count() );
901 }
902 }
903
904 rasterDataProvider->setEditable( false );
905
906 if ( feedback )
907 feedback->setProgress( maxProgressDuringBlockWriting );
908
909 if ( feedback && hasReportsDuringClose )
910 {
911 std::unique_ptr<QgsFeedback> scaledFeedback( QgsFeedback::createScaledFeedback( feedback, maxProgressDuringBlockWriting, 100.0 ) );
912 if ( !rasterDataProvider->closeWithProgress( scaledFeedback.get() ) )
913 {
914 if ( feedback->isCanceled() )
915 return {};
916 throw QgsProcessingException( QObject::tr( "Could not write raster dataset" ) );
917 }
918 }
919
920 QVariantMap ret;
921 ret[u"OUTPUT"_s] = fileName;
922
923 return ret;
924}
925
926QSet<int> QgsMeshRasterizeAlgorithm::supportedDataType()
927{
929}
930
931QString QgsMeshContoursAlgorithm::name() const
932{
933 return u"meshcontours"_s;
934}
935
936QString QgsMeshContoursAlgorithm::displayName() const
937{
938 return QObject::tr( "Export contours" );
939}
940
941QString QgsMeshContoursAlgorithm::group() const
942{
943 return QObject::tr( "Mesh" );
944}
945
946QString QgsMeshContoursAlgorithm::groupId() const
947{
948 return u"mesh"_s;
949}
950
951QString QgsMeshContoursAlgorithm::shortHelpString() const
952{
953 return QObject::tr( "This algorithm creates contours as a vector layer from a mesh scalar dataset." );
954}
955
956QString QgsMeshContoursAlgorithm::shortDescription() const
957{
958 return QObject::tr( "Creates contours as vector layer from mesh scalar dataset." );
959}
960
961QgsProcessingAlgorithm *QgsMeshContoursAlgorithm::createInstance() const
962{
963 return new QgsMeshContoursAlgorithm();
964}
965
966void QgsMeshContoursAlgorithm::initAlgorithm( const QVariantMap &configuration )
967{
968 Q_UNUSED( configuration );
969
970 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
971
972 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
973
974 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
975
976 addParameter( new QgsProcessingParameterNumber( u"INCREMENT"_s, QObject::tr( "Increment between contour levels" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
977
978 addParameter( new QgsProcessingParameterNumber( u"MINIMUM"_s, QObject::tr( "Minimum contour level" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
979 addParameter( new QgsProcessingParameterNumber( u"MAXIMUM"_s, QObject::tr( "Maximum contour level" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
980
981 auto contourLevelList = std::make_unique<QgsProcessingParameterString>( u"CONTOUR_LEVEL_LIST"_s, QObject::tr( "List of contours level" ), QVariant(), false, true );
982 contourLevelList->setHelp( QObject::tr( "Comma separated list of values to export. If filled, the increment, minimum and maximum settings are ignored." ) );
983 addParameter( contourLevelList.release() );
984
985 addParameter( new QgsProcessingParameterCrs( u"CRS_OUTPUT"_s, QObject::tr( "Output coordinate system" ), QVariant(), true ) );
986
987
988 addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT_LINES"_s, QObject::tr( "Exported contour lines" ), Qgis::ProcessingSourceType::VectorLine ) );
989 addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT_POLYGONS"_s, QObject::tr( "Exported contour polygons" ), Qgis::ProcessingSourceType::VectorPolygon ) );
990}
991
992bool QgsMeshContoursAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
993{
994 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
995
996 if ( !meshLayer || !meshLayer->isValid() )
997 return false;
998
999 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
1000 if ( !outputCrs.isValid() )
1001 outputCrs = meshLayer->crs();
1002 mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() );
1003 if ( !meshLayer->nativeMesh() )
1004 meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh
1005
1006 mTriangularMesh.update( meshLayer->nativeMesh(), mTransform );
1007 mNativeMesh = *meshLayer->nativeMesh();
1008
1009 // Prepare levels
1010 mLevels.clear();
1011 // First, try with the levels list
1012 QString levelsString = parameterAsString( parameters, u"CONTOUR_LEVEL_LIST"_s, context );
1013 if ( !levelsString.isEmpty() )
1014 {
1015 QStringList levelStringList = levelsString.split( ',' );
1016 if ( !levelStringList.isEmpty() )
1017 {
1018 for ( const QString &stringVal : levelStringList )
1019 {
1020 bool ok;
1021 double val = stringVal.toDouble( &ok );
1022 if ( ok )
1023 mLevels.append( val );
1024 else
1025 throw QgsProcessingException( QObject::tr( "Invalid format for level values, must be numbers separated with comma" ) );
1026
1027 if ( mLevels.count() >= 2 )
1028 if ( mLevels.last() <= mLevels.at( mLevels.count() - 2 ) )
1029 throw QgsProcessingException( QObject::tr( "Invalid format for level values, must be different numbers and in increasing order" ) );
1030 }
1031 }
1032 }
1033
1034 if ( mLevels.isEmpty() )
1035 {
1036 double minimum = parameterAsDouble( parameters, u"MINIMUM"_s, context );
1037 double maximum = parameterAsDouble( parameters, u"MAXIMUM"_s, context );
1038 double interval = parameterAsDouble( parameters, u"INCREMENT"_s, context );
1039
1040 if ( interval <= 0 )
1041 throw QgsProcessingException( QObject::tr( "Invalid interval value, must be greater than zero" ) );
1042
1043 if ( minimum >= maximum )
1044 throw QgsProcessingException( QObject::tr( "Invalid minimum and maximum values, minimum must be lesser than maximum" ) );
1045
1046 if ( interval > ( maximum - minimum ) )
1047 throw QgsProcessingException( QObject::tr( "Invalid minimum, maximum and interval values, difference between minimum and maximum must be greater or equal than interval" ) );
1048
1049 int intervalCount = ( maximum - minimum ) / interval;
1050
1051 mLevels.reserve( intervalCount );
1052 for ( int i = 0; i < intervalCount; ++i )
1053 {
1054 mLevels.append( minimum + i * interval );
1055 }
1056 }
1057
1058 // Prepare data
1059 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
1060
1061 if ( feedback )
1062 {
1063 feedback->setProgressText( QObject::tr( "Preparing data" ) );
1064 }
1065
1066 // Extract the date time used to export dataset values under a relative time
1067 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
1068 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
1069
1070 mDateTimeString = meshLayer->formatTime( relativeTime.hours() );
1071
1072 extractDatasetValues( datasetGroups, meshLayer, mNativeMesh, relativeTime, supportedDataType(), mDataPerGroup, feedback );
1073
1074 mLayerRendererSettings = meshLayer->rendererSettings();
1075
1076 return true;
1077}
1078
1079QVariantMap QgsMeshContoursAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1080{
1081 //First, if present, average 3D staked dataset value to 2D face value
1082 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
1083 for ( DataGroup &dataGroup : mDataPerGroup )
1084 {
1085 if ( dataGroup.dataset3dStakedValue.isValid() )
1086 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
1087 }
1088
1089 // Create vector layers
1090 QgsFields polygonFields;
1091 QgsFields lineFields;
1092 polygonFields.append( QgsField( QObject::tr( "group" ), QMetaType::Type::QString ) );
1093 polygonFields.append( QgsField( QObject::tr( "time" ), QMetaType::Type::QString ) );
1094 polygonFields.append( QgsField( QObject::tr( "min_value" ), QMetaType::Type::Double ) );
1095 polygonFields.append( QgsField( QObject::tr( "max_value" ), QMetaType::Type::Double ) );
1096 lineFields.append( QgsField( QObject::tr( "group" ), QMetaType::Type::QString ) );
1097 lineFields.append( QgsField( QObject::tr( "time" ), QMetaType::Type::QString ) );
1098 lineFields.append( QgsField( QObject::tr( "value" ), QMetaType::Type::Double ) );
1099
1100 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
1101
1102 QString lineIdentifier;
1103 QString polygonIdentifier;
1104 std::unique_ptr<QgsFeatureSink> sinkPolygons( parameterAsSink( parameters, u"OUTPUT_POLYGONS"_s, context, polygonIdentifier, polygonFields, Qgis::WkbType::PolygonZ, outputCrs ) );
1105 std::unique_ptr<QgsFeatureSink> sinkLines( parameterAsSink( parameters, u"OUTPUT_LINES"_s, context, lineIdentifier, lineFields, Qgis::WkbType::LineStringZ, outputCrs ) );
1106
1107 if ( !sinkLines || !sinkPolygons )
1108 return QVariantMap();
1109
1110
1111 for ( int i = 0; i < mDataPerGroup.count(); ++i )
1112 {
1113 DataGroup dataGroup = mDataPerGroup.at( i );
1114 bool scalarDataOnVertices = dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
1115 int count = scalarDataOnVertices ? mNativeMesh.vertices.count() : mNativeMesh.faces.count();
1116
1117 QVector<double> values;
1118 if ( dataGroup.datasetValues.isValid() )
1119 {
1120 // vals could be scalar or vectors, for contour rendering we want always magnitude
1121 values = QgsMeshLayerUtils::calculateMagnitudes( dataGroup.datasetValues );
1122 }
1123 else
1124 {
1125 values = QVector<double>( count, std::numeric_limits<double>::quiet_NaN() );
1126 }
1127
1128 if ( ( !scalarDataOnVertices ) )
1129 {
1130 values = QgsMeshLayerUtils::interpolateFromFacesData( values, mNativeMesh, &dataGroup.activeFaces, QgsMeshRendererScalarSettings::NeighbourAverage );
1131 }
1132
1133 QgsMeshContours contoursExported( mTriangularMesh, mNativeMesh, values, dataGroup.activeFaces );
1134
1135 QgsAttributes firstAttributes;
1136 firstAttributes.append( dataGroup.metadata.name() );
1137 firstAttributes.append( mDateTimeString );
1138
1139 for ( double level : std::as_const( mLevels ) )
1140 {
1141 QgsGeometry line = contoursExported.exportLines( level, feedback );
1142 if ( feedback->isCanceled() )
1143 return QVariantMap();
1144 if ( line.isEmpty() )
1145 continue;
1146 QgsAttributes lineAttributes = firstAttributes;
1147 lineAttributes.append( level );
1148
1149 QgsFeature lineFeat;
1150 lineFeat.setGeometry( line );
1151 lineFeat.setAttributes( lineAttributes );
1152
1153 if ( !sinkLines->addFeature( lineFeat, QgsFeatureSink::FastInsert ) )
1154 throw QgsProcessingException( writeFeatureError( sinkLines.get(), parameters, u"OUTPUT_LINES"_s ) );
1155 else
1156 feedback->featureAddedToSink( u"OUTPUT_LINES"_s );
1157 }
1158
1159 for ( int l = 0; l < mLevels.count() - 1; ++l )
1160 {
1161 QgsGeometry polygon = contoursExported.exportPolygons( mLevels.at( l ), mLevels.at( l + 1 ), feedback );
1162 if ( feedback->isCanceled() )
1163 return QVariantMap();
1164
1165 if ( polygon.isEmpty() )
1166 continue;
1167 QgsAttributes polygonAttributes = firstAttributes;
1168 polygonAttributes.append( mLevels.at( l ) );
1169 polygonAttributes.append( mLevels.at( l + 1 ) );
1170
1171 QgsFeature polygonFeature;
1172 polygonFeature.setGeometry( polygon );
1173 polygonFeature.setAttributes( polygonAttributes );
1174 if ( !sinkPolygons->addFeature( polygonFeature ) )
1175 {
1176 throw QgsProcessingException( writeFeatureError( sinkPolygons.get(), parameters, QString() ) );
1177 }
1178 else
1179 {
1180 feedback->featureAddedToSink( u"OUTPUT_POLYGONS"_s );
1181 }
1182 }
1183
1184 if ( feedback )
1185 {
1186 feedback->setProgress( 100 * i / mDataPerGroup.count() );
1187 }
1188 }
1189
1190 if ( sinkPolygons )
1191 {
1192 sinkPolygons->finalize();
1193 feedback->featureSinkFinalized( u"OUTPUT_POLYGONS"_s );
1194 }
1195 if ( sinkLines )
1196 {
1197 sinkLines->finalize();
1198 feedback->featureSinkFinalized( u"OUTPUT_LINES"_s );
1199 }
1200
1201 QVariantMap ret;
1202 ret[u"OUTPUT_LINES"_s] = lineIdentifier;
1203 ret[u"OUTPUT_POLYGONS"_s] = polygonIdentifier;
1204
1205 return ret;
1206}
1207
1208QString QgsMeshExportCrossSection::name() const
1209{
1210 return u"meshexportcrosssection"_s;
1211}
1212
1213QString QgsMeshExportCrossSection::displayName() const
1214{
1215 return QObject::tr( "Export cross section dataset values on lines from mesh" );
1216}
1217
1218QString QgsMeshExportCrossSection::group() const
1219{
1220 return QObject::tr( "Mesh" );
1221}
1222
1223QString QgsMeshExportCrossSection::groupId() const
1224{
1225 return u"mesh"_s;
1226}
1227
1228QString QgsMeshExportCrossSection::shortHelpString() const
1229{
1230 return QObject::tr(
1231 "This algorithm extracts mesh's dataset values from line contained in a vector layer.\n"
1232 "Each line is discretized with a resolution distance parameter for extraction of values on its vertices."
1233 );
1234}
1235
1236QString QgsMeshExportCrossSection::shortDescription() const
1237{
1238 return QObject::tr( "Extracts a mesh dataset's values from lines contained in a vector layer." );
1239}
1240
1241QgsProcessingAlgorithm *QgsMeshExportCrossSection::createInstance() const
1242{
1243 return new QgsMeshExportCrossSection();
1244}
1245
1246void QgsMeshExportCrossSection::initAlgorithm( const QVariantMap &configuration )
1247{
1248 Q_UNUSED( configuration );
1249
1250 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
1251
1252 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
1253
1254 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
1255
1256 QList<int> datatype;
1257 datatype << static_cast<int>( Qgis::ProcessingSourceType::VectorLine );
1258 addParameter( new QgsProcessingParameterFeatureSource( u"INPUT_LINES"_s, QObject::tr( "Lines for data export" ), datatype, QVariant(), false ) );
1259
1260 addParameter( new QgsProcessingParameterDistance( u"RESOLUTION"_s, QObject::tr( "Line segmentation resolution" ), 10.0, u"INPUT_LINES"_s, false, 0 ) );
1261
1262 addParameter( new QgsProcessingParameterNumber( u"COORDINATES_DIGITS"_s, QObject::tr( "Digits count for coordinates" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1263
1264 addParameter( new QgsProcessingParameterNumber( u"DATASET_DIGITS"_s, QObject::tr( "Digits count for dataset value" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1265
1266 addParameter( new QgsProcessingParameterFileDestination( u"OUTPUT"_s, QObject::tr( "Exported data CSV file" ), QObject::tr( "CSV file (*.csv)" ) ) );
1267}
1268
1269bool QgsMeshExportCrossSection::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1270{
1271 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
1272
1273 if ( !meshLayer || !meshLayer->isValid() )
1274 return false;
1275
1276 mMeshLayerCrs = meshLayer->crs();
1277 mTriangularMesh.update( meshLayer->nativeMesh() );
1278 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
1279
1280 if ( feedback )
1281 {
1282 feedback->setProgressText( QObject::tr( "Preparing data" ) );
1283 }
1284
1285 // Extract the date time used to export dataset values under a relative time
1286 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
1287 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
1288
1289 extractDatasetValues( datasetGroups, meshLayer, *meshLayer->nativeMesh(), relativeTime, supportedDataType(), mDataPerGroup, feedback );
1290
1291 mLayerRendererSettings = meshLayer->rendererSettings();
1292
1293 return true;
1294}
1295
1296QVariantMap QgsMeshExportCrossSection::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1297{
1298 if ( feedback )
1299 feedback->setProgress( 0 );
1300 //First, if present, average 3D staked dataset value to 2D face value
1301 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
1302 for ( DataGroup &dataGroup : mDataPerGroup )
1303 {
1304 if ( dataGroup.dataset3dStakedValue.isValid() )
1305 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
1306 }
1307 double resolution = parameterAsDouble( parameters, u"RESOLUTION"_s, context );
1308 int datasetDigits = parameterAsInt( parameters, u"DATASET_DIGITS"_s, context );
1309 int coordDigits = parameterAsInt( parameters, u"COORDINATES_DIGITS"_s, context );
1310
1311 std::unique_ptr<QgsProcessingFeatureSource> featureSource( parameterAsSource( parameters, u"INPUT_LINES"_s, context ) );
1312 if ( !featureSource )
1313 throw QgsProcessingException( QObject::tr( "Input lines vector layer required" ) );
1314
1315 QgsCoordinateTransform transform( featureSource->sourceCrs(), mMeshLayerCrs, context.transformContext() );
1316
1317 QString outputFileName = parameterAsFileOutput( parameters, u"OUTPUT"_s, context );
1318 QFile file( outputFileName );
1319 if ( !file.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
1320 throw QgsProcessingException( QObject::tr( "Unable to create the output file" ) );
1321
1322 QTextStream textStream( &file );
1323 QStringList header;
1324 header << u"fid"_s << u"x"_s << u"y"_s << QObject::tr( "offset" );
1325 for ( const DataGroup &datagroup : std::as_const( mDataPerGroup ) )
1326 header << datagroup.metadata.name();
1327 textStream << header.join( ',' ) << u"\n"_s;
1328
1329 long long featCount = featureSource->featureCount();
1330 long long featCounter = 0;
1331 QgsFeatureIterator featIt = featureSource->getFeatures();
1332 QgsFeature feat;
1333 while ( featIt.nextFeature( feat ) )
1334 {
1335 QgsFeatureId fid = feat.id();
1336 QgsGeometry line = feat.geometry();
1337 try
1338 {
1339 line.transform( transform );
1340 }
1341 catch ( QgsCsException & )
1342 {
1343 line = feat.geometry();
1344 feedback->reportError( QObject::tr( "Could not transform line to mesh CRS" ) );
1345 }
1346
1347 if ( line.isEmpty() )
1348 continue;
1349 double offset = 0;
1350 while ( offset <= line.length() )
1351 {
1352 if ( feedback->isCanceled() )
1353 return QVariantMap();
1354
1355 QStringList textLine;
1356 QgsPointXY point = line.interpolate( offset ).asPoint();
1357 int triangularFaceIndex = mTriangularMesh.faceIndexForPoint_v2( point );
1358 textLine << QString::number( fid ) << QString::number( point.x(), 'f', coordDigits ) << QString::number( point.y(), 'f', coordDigits ) << QString::number( offset, 'f', coordDigits );
1359 if ( triangularFaceIndex >= 0 )
1360 {
1361 //extract dataset values for the point
1362 QgsAttributes attributes;
1363 int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces().at( triangularFaceIndex );
1364 for ( int i = 0; i < mDataPerGroup.count(); ++i )
1365 {
1366 const DataGroup &dataGroup = mDataPerGroup.at( i );
1367 bool faceActive = dataGroup.activeFaces.active( nativeFaceIndex );
1368 if ( !faceActive )
1369 continue;
1370 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
1371
1372 if ( abs( value.x() ) == std::numeric_limits<double>::quiet_NaN() )
1373 textLine << QString( ' ' );
1374 else
1375 textLine << QString::number( value.scalar(), 'f', datasetDigits );
1376 }
1377 }
1378 else
1379 for ( int i = 0; i < mDataPerGroup.count(); ++i )
1380 textLine << QString( ' ' );
1381
1382 textStream << textLine.join( ',' ) << u"\n"_s;
1383
1384 offset += resolution;
1385 }
1386
1387 if ( feedback )
1388 {
1389 feedback->setProgress( 100.0 * featCounter / featCount );
1390 if ( feedback->isCanceled() )
1391 return QVariantMap();
1392 }
1393 }
1394
1395 file.close();
1396
1397 QVariantMap ret;
1398 ret[u"OUTPUT"_s] = outputFileName;
1399 return ret;
1400}
1401
1402QString QgsMeshExportTimeSeries::name() const
1403{
1404 return u"meshexporttimeseries"_s;
1405}
1406
1407QString QgsMeshExportTimeSeries::displayName() const
1408{
1409 return QObject::tr( "Export time series values from points of a mesh dataset" );
1410}
1411
1412QString QgsMeshExportTimeSeries::group() const
1413{
1414 return QObject::tr( "Mesh" );
1415}
1416
1417QString QgsMeshExportTimeSeries::groupId() const
1418{
1419 return u"mesh"_s;
1420}
1421
1422QString QgsMeshExportTimeSeries::shortHelpString() const
1423{
1424 return QObject::tr(
1425 "This algorithm extracts mesh's dataset time series values from points contained in a vector layer.\n"
1426 "If the time step is kept to its default value (0 hours), the time step used is the one of the two first datasets of the first selected dataset group."
1427 );
1428}
1429
1430QString QgsMeshExportTimeSeries::shortDescription() const
1431{
1432 return QObject::tr( "Extracts a mesh dataset's time series values from points contained in a vector layer." );
1433}
1434
1435QgsProcessingAlgorithm *QgsMeshExportTimeSeries::createInstance() const
1436{
1437 return new QgsMeshExportTimeSeries();
1438}
1439
1440void QgsMeshExportTimeSeries::initAlgorithm( const QVariantMap &configuration )
1441{
1442 Q_UNUSED( configuration );
1443
1444 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
1445
1446 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
1447
1448 addParameter( new QgsProcessingParameterMeshDatasetTime( u"STARTING_TIME"_s, QObject::tr( "Starting time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
1449
1450 addParameter( new QgsProcessingParameterMeshDatasetTime( u"FINISHING_TIME"_s, QObject::tr( "Finishing time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
1451
1452 addParameter( new QgsProcessingParameterNumber( u"TIME_STEP"_s, QObject::tr( "Time step (hours)" ), Qgis::ProcessingNumberParameterType::Double, 0, true, 0 ) );
1453
1454 QList<int> datatype;
1455 datatype << static_cast<int>( Qgis::ProcessingSourceType::VectorPoint );
1456 addParameter( new QgsProcessingParameterFeatureSource( u"INPUT_POINTS"_s, QObject::tr( "Points for data export" ), datatype, QVariant(), false ) );
1457
1458 addParameter( new QgsProcessingParameterNumber( u"COORDINATES_DIGITS"_s, QObject::tr( "Digits count for coordinates" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1459
1460 addParameter( new QgsProcessingParameterNumber( u"DATASET_DIGITS"_s, QObject::tr( "Digits count for dataset value" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1461
1462 addParameter( new QgsProcessingParameterFileDestination( u"OUTPUT"_s, QObject::tr( "Exported data CSV file" ), QObject::tr( "CSV file (*.csv)" ) ) );
1463}
1464
1465bool QgsMeshExportTimeSeries::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1466{
1467 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
1468
1469 if ( !meshLayer || !meshLayer->isValid() )
1470 return false;
1471
1472 mMeshLayerCrs = meshLayer->crs();
1473 mTriangularMesh.update( meshLayer->nativeMesh() );
1474
1475 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
1476
1477 if ( feedback )
1478 {
1479 feedback->setProgressText( QObject::tr( "Preparing data" ) );
1480 }
1481
1482 // Extract the date times used to export dataset values
1483 QVariant parameterStartTimeVariant = parameters.value( u"STARTING_TIME"_s );
1484 QgsInterval relativeStartTime = datasetRelativetime( parameterStartTimeVariant, meshLayer, context );
1485
1486 QVariant parameterEndTimeVariant = parameters.value( u"FINISHING_TIME"_s );
1487 QgsInterval relativeEndTime = datasetRelativetime( parameterEndTimeVariant, meshLayer, context );
1488
1489 // calculate time steps
1490 qint64 timeStepInterval = parameterAsDouble( parameters, u"TIME_STEP"_s, context ) * 1000 * 3600;
1491 if ( timeStepInterval == 0 )
1492 {
1493 //take the first time step of the first temporal dataset group
1494 for ( int groupIndex : datasetGroups )
1495 {
1496 QgsMeshDatasetGroupMetadata meta = meshLayer->datasetGroupMetadata( QgsMeshDatasetIndex( groupIndex, 0 ) );
1497 if ( !meta.isTemporal() && meshLayer->datasetCount( QgsMeshDatasetIndex( groupIndex, 0 ) ) < 2 )
1498 continue;
1499 else
1500 {
1501 timeStepInterval = meshLayer->datasetRelativeTimeInMilliseconds( QgsMeshDatasetIndex( groupIndex, 1 ) ) - meshLayer->datasetRelativeTimeInMilliseconds( QgsMeshDatasetIndex( groupIndex, 0 ) );
1502 break;
1503 }
1504 }
1505 }
1506
1507 mRelativeTimeSteps.clear();
1508 mTimeStepString.clear();
1509 if ( timeStepInterval != 0 )
1510 {
1511 mRelativeTimeSteps.append( relativeStartTime.seconds() * 1000 );
1512 while ( mRelativeTimeSteps.last() < relativeEndTime.seconds() * 1000 )
1513 mRelativeTimeSteps.append( mRelativeTimeSteps.last() + timeStepInterval );
1514
1515 for ( qint64 relativeTimeStep : std::as_const( mRelativeTimeSteps ) )
1516 {
1517 mTimeStepString.append( meshLayer->formatTime( relativeTimeStep / 3600.0 / 1000.0 ) );
1518 }
1519 }
1520
1521 //Extract needed dataset values
1522 for ( int i = 0; i < datasetGroups.count(); ++i )
1523 {
1524 int groupIndex = datasetGroups.at( i );
1525 QgsMeshDatasetGroupMetadata meta = meshLayer->datasetGroupMetadata( QgsMeshDatasetIndex( groupIndex, 0 ) );
1526 if ( supportedDataType().contains( meta.dataType() ) )
1527 {
1528 mGroupIndexes.append( groupIndex );
1529 mGroupsMetadata[groupIndex] = meta;
1530 int valueCount = meta.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices ? mTriangularMesh.vertices().count() : meshLayer->nativeMesh()->faceCount();
1531
1532 if ( !mRelativeTimeSteps.isEmpty() )
1533 {
1534 //QMap<qint64, DataGroup> temporalGroup;
1535 QgsMeshDatasetIndex lastDatasetIndex;
1536 for ( qint64 relativeTimeStep : std::as_const( mRelativeTimeSteps ) )
1537 {
1538 QMap<int, int> &groupIndexToData = mRelativeTimeToData[relativeTimeStep];
1539 QgsInterval timeStepInterval( relativeTimeStep / 1000.0 );
1540 QgsMeshDatasetIndex datasetIndex = meshLayer->datasetIndexAtRelativeTime( timeStepInterval, groupIndex );
1541 if ( !datasetIndex.isValid() )
1542 continue;
1543 if ( datasetIndex != lastDatasetIndex )
1544 {
1545 DataGroup dataGroup;
1546 dataGroup.metadata = meta;
1547 dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, valueCount );
1548 dataGroup.activeFaces = meshLayer->areFacesActive( datasetIndex, 0, meshLayer->nativeMesh()->faceCount() );
1549 if ( dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVolumes )
1550 {
1551 dataGroup.dataset3dStakedValue = meshLayer->dataset3dValues( datasetIndex, 0, valueCount );
1552 }
1553 mDatasets.append( dataGroup );
1554 lastDatasetIndex = datasetIndex;
1555 }
1556 groupIndexToData[groupIndex] = mDatasets.count() - 1;
1557 }
1558 }
1559 else
1560 {
1561 // we have only static dataset group
1562 QMap<int, int> &groupIndexToData = mRelativeTimeToData[0];
1563 QgsMeshDatasetIndex datasetIndex( groupIndex, 0 );
1564 DataGroup dataGroup;
1565 dataGroup.metadata = meta;
1566 dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, valueCount );
1567 dataGroup.activeFaces = meshLayer->areFacesActive( datasetIndex, 0, meshLayer->nativeMesh()->faceCount() );
1568 if ( dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVolumes )
1569 {
1570 dataGroup.dataset3dStakedValue = meshLayer->dataset3dValues( datasetIndex, 0, valueCount );
1571 }
1572 mDatasets.append( dataGroup );
1573 groupIndexToData[groupIndex] = mDatasets.count() - 1;
1574 }
1575 }
1576
1577 if ( feedback )
1578 feedback->setProgress( 100 * i / datasetGroups.count() );
1579 }
1580
1581 mLayerRendererSettings = meshLayer->rendererSettings();
1582
1583 return true;
1584}
1585
1586
1587QVariantMap QgsMeshExportTimeSeries::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1588{
1589 if ( feedback )
1590 feedback->setProgress( 0 );
1591 //First, if present, average 3D staked dataset value to 2D face value
1592 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
1593
1594 for ( DataGroup &dataGroup : mDatasets )
1595 {
1596 if ( dataGroup.dataset3dStakedValue.isValid() )
1597 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
1598 }
1599
1600 int datasetDigits = parameterAsInt( parameters, u"DATASET_DIGITS"_s, context );
1601 int coordDigits = parameterAsInt( parameters, u"COORDINATES_DIGITS"_s, context );
1602
1603 std::unique_ptr<QgsProcessingFeatureSource> featureSource( parameterAsSource( parameters, u"INPUT_POINTS"_s, context ) );
1604 if ( !featureSource )
1605 throw QgsProcessingException( QObject::tr( "Input points vector layer required" ) );
1606
1607 QgsCoordinateTransform transform( featureSource->sourceCrs(), mMeshLayerCrs, context.transformContext() );
1608
1609 QString outputFileName = parameterAsFileOutput( parameters, u"OUTPUT"_s, context );
1610 QFile file( outputFileName );
1611 if ( !file.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
1612 throw QgsProcessingException( QObject::tr( "Unable to create the output file" ) );
1613
1614 QTextStream textStream( &file );
1615 QStringList header;
1616 header << u"fid"_s << u"x"_s << u"y"_s << QObject::tr( "time" );
1617
1618 for ( int gi : std::as_const( mGroupIndexes ) )
1619 header << mGroupsMetadata.value( gi ).name();
1620
1621 textStream << header.join( ',' ) << u"\n"_s;
1622
1623 long long featCount = featureSource->featureCount();
1624 long long featCounter = 0;
1625 QgsFeatureIterator featIt = featureSource->getFeatures();
1626 QgsFeature feat;
1627 while ( featIt.nextFeature( feat ) )
1628 {
1629 QgsFeatureId fid = feat.id();
1630 QgsGeometry geom = feat.geometry();
1631 try
1632 {
1633 geom.transform( transform );
1634 }
1635 catch ( QgsCsException & )
1636 {
1637 geom = feat.geometry();
1638 feedback->reportError( QObject::tr( "Could not transform line to mesh CRS" ) );
1639 }
1640
1641 if ( geom.isEmpty() )
1642 continue;
1643
1644 QgsPointXY point = geom.asPoint();
1645 int triangularFaceIndex = mTriangularMesh.faceIndexForPoint_v2( point );
1646
1647 if ( triangularFaceIndex >= 0 )
1648 {
1649 int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces().at( triangularFaceIndex );
1650 if ( !mRelativeTimeSteps.isEmpty() )
1651 {
1652 for ( int timeIndex = 0; timeIndex < mRelativeTimeSteps.count(); ++timeIndex )
1653 {
1654 qint64 timeStep = mRelativeTimeSteps.at( timeIndex );
1655 QStringList textLine;
1656 textLine << QString::number( fid ) << QString::number( point.x(), 'f', coordDigits ) << QString::number( point.y(), 'f', coordDigits ) << mTimeStepString.at( timeIndex );
1657
1658 if ( mRelativeTimeToData.contains( timeStep ) )
1659 {
1660 const QMap<int, int> &groupToData = mRelativeTimeToData.value( timeStep );
1661 for ( int groupIndex : std::as_const( mGroupIndexes ) )
1662 {
1663 if ( !groupToData.contains( groupIndex ) )
1664 continue;
1665 int dataIndex = groupToData.value( groupIndex );
1666 if ( dataIndex < 0 || dataIndex > mDatasets.count() - 1 )
1667 continue;
1668
1669 const DataGroup &dataGroup = mDatasets.at( dataIndex );
1670 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
1671 if ( abs( value.x() ) == std::numeric_limits<double>::quiet_NaN() )
1672 textLine << QString( ' ' );
1673 else
1674 textLine << QString::number( value.scalar(), 'f', datasetDigits );
1675 }
1676 }
1677 textStream << textLine.join( ',' ) << u"\n"_s;
1678 }
1679 }
1680 else
1681 {
1682 QStringList textLine;
1683 textLine << QString::number( fid ) << QString::number( point.x(), 'f', coordDigits ) << QString::number( point.y(), 'f', coordDigits ) << QObject::tr( "static dataset" );
1684 const QMap<int, int> &groupToData = mRelativeTimeToData.value( 0 );
1685 for ( int groupIndex : std::as_const( mGroupIndexes ) )
1686 {
1687 if ( !groupToData.contains( groupIndex ) )
1688 continue;
1689 int dataIndex = groupToData.value( groupIndex );
1690 if ( dataIndex < 0 || dataIndex > mDatasets.count() - 1 )
1691 continue;
1692 const DataGroup &dataGroup = mDatasets.at( dataIndex );
1693 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
1694 if ( abs( value.x() ) == std::numeric_limits<double>::quiet_NaN() )
1695 textLine << QString( ' ' );
1696 else
1697 textLine << QString::number( value.scalar(), 'f', datasetDigits );
1698 }
1699 textStream << textLine.join( ',' ) << u"\n"_s;
1700 }
1701 }
1702 featCounter++;
1703 if ( feedback )
1704 {
1705 feedback->setProgress( 100.0 * featCounter / featCount );
1706 if ( feedback->isCanceled() )
1707 return QVariantMap();
1708 }
1709 }
1710
1711 file.close();
1712
1713 QVariantMap ret;
1714 ret[u"OUTPUT"_s] = outputFileName;
1715 return ret;
1716}
1717
@ VectorPoint
Vector point layers.
Definition qgis.h:3715
@ VectorPolygon
Vector polygon layers.
Definition qgis.h:3717
@ VectorLine
Vector line layers.
Definition qgis.h:3716
@ Float64
Sixty four bit floating point (double).
Definition qgis.h:402
@ Point
Point.
Definition qgis.h:296
@ LineStringZ
LineStringZ.
Definition qgis.h:314
@ PolygonZ
PolygonZ.
Definition qgis.h:315
@ Hidden
Parameter is hidden and should not be shown to users.
Definition qgis.h:3948
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
Definition qgis.h:3947
@ Double
Double/float values.
Definition qgis.h:3988
A vector of attributes.
Represents a coordinate reference system (CRS).
bool isValid() const
Returns whether this CRS is correctly initialized and usable.
Handles coordinate transforms between two coordinate systems.
Custom exception class for Coordinate Reference System related exceptions.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:60
QgsFeatureId id
Definition qgsfeature.h:63
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
QgsGeometry geometry
Definition qgsfeature.h:66
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:56
void canceled()
Internal routines can connect to this signal if they use event loop.
void cancel()
Tells the internal routines that the current operation should be canceled. This should be run by the ...
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:65
static std::unique_ptr< QgsFeedback > createScaledFeedback(QgsFeedback *parentFeedback, double startPercentage, double endPercentage)
Returns a feedback object whose [0, 100] progression range will be mapped to parentFeedback [startPer...
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:56
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:75
A geometry is the spatial representation of a feature.
double length() const
Returns the planar, 2-dimensional length of geometry.
Qgis::GeometryOperationResult transform(const QgsCoordinateTransform &ct, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward, bool transformZ=false)
Transforms this geometry as described by the coordinate transform ct.
QgsGeometry interpolate(double distance) const
Returns an interpolated point on the geometry at the specified distance.
QgsPointXY asPoint() const
Returns the contents of the geometry as a 2-dimensional point.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
A representation of the interval between two datetime values.
Definition qgsinterval.h:52
double seconds() const
Returns the interval duration in seconds.
double hours() const
Returns the interval duration in hours.
Line string geometry type, with support for z-dimension and m-values.
QgsCoordinateReferenceSystem crs
Definition qgsmaplayer.h:90
Abstract class for interpolating 3d stacked mesh data to 2d data.
QgsMeshDataBlock calculate(const QgsMesh3DDataBlock &block3d, QgsFeedback *feedback=nullptr) const
Calculated 2d block values from 3d stacked mesh values.
Exporter of contours lines or polygons from a mesh layer.
A block of integers/doubles from a mesh dataset.
QgsMeshDatasetValue value(int index) const
Returns a value represented by the index For active flag the behavior is undefined.
bool active(int index) const
Returns a value for active flag by the index For scalar and vector 2d the behavior is undefined.
int count() const
Number of items stored in the block.
A collection of dataset group metadata such as whether the data is vector or scalar,...
bool isTemporal() const
Returns whether the dataset group is temporal (contains time-related dataset).
bool isVector() const
Returns whether dataset group has vector data.
DataType dataType() const
Returns whether dataset group data is defined on vertices or faces or volumes.
@ DataOnEdges
Data is defined on edges.
@ DataOnFaces
Data is defined on faces.
@ DataOnVertices
Data is defined on vertices.
@ DataOnVolumes
Data is defined on volumes.
An index that identifies the dataset group (e.g.
bool isValid() const
Returns whether index is valid, ie at least groups is set.
Represents a single mesh dataset value.
double y() const
Returns y value.
double scalar() const
Returns magnitude of vector for vector data or scalar value for scalar data.
double x() const
Returns x value.
Implementation of map layer temporal properties for mesh layers.
QDateTime referenceTime() const
Returns the reference time.
Represents a mesh layer supporting display of data on structured or unstructured meshes.
int datasetCount(const QgsMeshDatasetIndex &index) const
Returns the dataset count in the dataset groups.
QgsMeshRendererSettings rendererSettings() const
Returns renderer settings.
void updateTriangularMesh(const QgsCoordinateTransform &transform=QgsCoordinateTransform())
Gets native mesh and updates (creates if it doesn't exist) the base triangular mesh.
QgsMesh * nativeMesh()
Returns native mesh (nullptr before rendering or calling to updateMesh).
QgsMeshDatasetIndex datasetIndexAtRelativeTime(const QgsInterval &relativeTime, int datasetGroupIndex) const
Returns dataset index from datasets group depending on the relative time from the layer reference tim...
QgsMeshDataBlock datasetValues(const QgsMeshDatasetIndex &index, int valueIndex, int count) const
Returns N vector/scalar values from the index from the dataset.
bool isEditable() const override
Returns true if the layer can be edited.
QgsMeshDataBlock areFacesActive(const QgsMeshDatasetIndex &index, int faceIndex, int count) const
Returns whether the faces are active for particular dataset.
QgsInterval datasetRelativeTime(const QgsMeshDatasetIndex &index)
Returns the relative time of the dataset from the reference time of its group.
QgsMapLayerTemporalProperties * temporalProperties() override
Returns the layer's temporal properties.
qint64 datasetRelativeTimeInMilliseconds(const QgsMeshDatasetIndex &index)
Returns the relative time (in milliseconds) of the dataset from the reference time of its group.
QgsMesh3DDataBlock dataset3dValues(const QgsMeshDatasetIndex &index, int faceIndex, int count) const
Returns N vector/scalar values from the face index from the dataset for 3d stacked meshes.
QString formatTime(double hours)
Returns (date) time in hours formatted to human readable form.
QgsMeshDatasetGroupMetadata datasetGroupMetadata(const QgsMeshDatasetIndex &index) const
Returns the dataset groups metadata.
@ NeighbourAverage
Does a simple average of values defined for all surrounding faces/vertices.
static QgsRasterBlock * exportRasterBlock(const QgsMeshLayer &layer, const QgsMeshDatasetIndex &datasetIndex, const QgsCoordinateReferenceSystem &destinationCrs, const QgsCoordinateTransformContext &transformContext, double mapUnitsPerPixel, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Exports mesh layer's dataset values as raster block.
Represents a 2D point.
Definition qgspointxy.h:62
double y
Definition qgspointxy.h:66
double x
Definition qgspointxy.h:65
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:53
Abstract base class for processing algorithms.
Contains information about the context in which a processing algorithm is executed.
QgsDateTimeRange currentTimeRange() const
Returns the current time range to use for temporal operations.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
void featureAddedToSink(const QString &output)
Reports that a feature was added to the the sink associated with the specified algorithm output.
void featureSinkFinalized(const QString &output)
Reports that a feature sink has been finalized.
virtual void reportError(const QString &error, bool fatalError=false)
Reports that the algorithm encountered an error while executing.
virtual void setProgressText(const QString &text)
Sets a progress report text string.
A coordinate reference system parameter for processing algorithms.
A double numeric parameter for distance values.
An enum based parameter for processing algorithms, allowing for selection from predefined values.
A rectangular map extent parameter for processing algorithms.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A generic file based destination parameter, for specifying the destination path for a file (non-map l...
A parameter for processing algorithms that need a list of mesh dataset groups.
static QList< int > valueAsDatasetGroup(const QVariant &value)
Returns the value as a list if dataset group indexes.
A parameter for processing algorithms that need a list of mesh dataset index from time parameter.
static QString valueAsTimeType(const QVariant &value)
Returns the dataset value time type as a string : current-context-time : the time is store in the pro...
static QgsMeshDatasetIndex timeValueAsDatasetIndex(const QVariant &value)
Returns the value as a QgsMeshDatasetIndex if the value has "dataset-time-step" type.
static QDateTime timeValueAsDefinedDateTime(const QVariant &value)
Returns the value as a QDateTime if the value has "defined-date-time" type.
A mesh layer parameter for processing algorithms.
A numeric parameter for processing algorithms.
A raster layer destination parameter, for specifying the destination path for a raster layer created ...
Feedback object tailored for raster block reading.
The raster file writer which allows you to save a raster to a new file.
A rectangle specified with double values.
double xMinimum
double yMinimum
T begin() const
Returns the beginning of the range.
Definition qgsrange.h:408
A triangular/derived mesh with vertices in map coordinates.
const QVector< QgsMeshFace > & triangles() const
Returns triangles.
const QVector< QgsMeshVertex > & vertices() const
Returns vertices in map coordinate system.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference).
Definition qgis.h:7236
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
QVector< int > QgsMeshFace
List of vertex indexes.
QPair< int, int > QgsMeshEdge
Edge is a straight line seqment between 2 points.
Mesh - vertices, edges and faces.
QVector< QgsMeshVertex > vertices
void clear()
Remove all vertices, edges and faces.
int faceCount() const
Returns number of faces.