QGIS API Documentation 4.3.0-Master (18b5e825726)
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, feedback );
836 }
837
838 if ( feedback && feedback->isCanceled() )
839 return {};
840
841 // create raster
842 const double pixelSize = parameterAsDouble( parameters, u"PIXEL_SIZE"_s, context );
843 if ( qgsDoubleNear( pixelSize, 0 ) )
844 {
845 throw QgsProcessingException( QObject::tr( "Pixel size cannot be 0" ) );
846 }
847
848 QgsRectangle extent = parameterAsExtent( parameters, u"EXTENT"_s, context );
849 if ( extent.isEmpty() )
850 extent = mTriangularMesh.extent();
851
852 int width = extent.width() / pixelSize;
853 int height = extent.height() / pixelSize;
854
855 QString creationOptions = parameterAsString( parameters, u"CREATION_OPTIONS"_s, context ).trimmed();
856 // handle backwards compatibility parameter CREATE_OPTIONS
857 const QString optionsString = parameterAsString( parameters, u"CREATE_OPTIONS"_s, context );
858 if ( !optionsString.isEmpty() )
859 creationOptions = optionsString;
860
861 const QString fileName = parameterAsOutputLayer( parameters, u"OUTPUT"_s, context );
862 const QString outputFormat = parameterAsOutputRasterFormat( parameters, u"OUTPUT"_s, context );
863 QgsRasterFileWriter rasterFileWriter( fileName );
864 rasterFileWriter.setOutputProviderKey( u"gdal"_s );
865 if ( !creationOptions.isEmpty() )
866 {
867 rasterFileWriter.setCreationOptions( creationOptions.split( '|' ) );
868 }
869 rasterFileWriter.setOutputFormat( outputFormat );
870
871 std::unique_ptr<QgsRasterDataProvider> rasterDataProvider( rasterFileWriter.createMultiBandRaster( Qgis::DataType::Float64, width, height, extent, mTransform.destinationCrs(), mDataPerGroup.count() ) );
872 if ( !rasterDataProvider )
873 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( fileName ) );
874 if ( !rasterDataProvider->isEditable() && !rasterDataProvider->setEditable( true ) )
875 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( rasterDataProvider->error().summary() ) );
876
877 const bool hasReportsDuringClose = rasterDataProvider->hasReportsDuringClose();
878 const double maxProgressDuringBlockWriting = hasReportsDuringClose ? 50.0 : 100.0;
879
880 for ( int i = 0; i < mDataPerGroup.count(); ++i )
881 {
882 const DataGroup &dataGroup = mDataPerGroup.at( i );
883 QgsRasterBlockFeedback rasterBlockFeedBack;
884 if ( feedback )
885 QObject::connect( feedback, &QgsFeedback::canceled, &rasterBlockFeedBack, &QgsRasterBlockFeedback::cancel, Qt::DirectConnection );
886
887 if ( dataGroup.datasetValues.isValid() )
888 {
889 std::unique_ptr<QgsRasterBlock> block(
890 QgsMeshUtils::exportRasterBlock( mTriangularMesh, dataGroup.datasetValues, dataGroup.activeFaces, dataGroup.metadata.dataType(), mTransform, pixelSize, extent, &rasterBlockFeedBack )
891 );
892
893 if ( feedback && feedback->isCanceled() )
894 return {};
895
896 if ( !rasterDataProvider->writeBlock( block.get(), i + 1 ) )
897 {
898 throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( rasterDataProvider->error().summary() ) );
899 }
900 rasterDataProvider->setNoDataValue( i + 1, block->noDataValue() );
901 }
902 else
903 rasterDataProvider->setNoDataValue( i + 1, std::numeric_limits<double>::quiet_NaN() );
904
905 if ( feedback )
906 {
907 if ( feedback->isCanceled() )
908 return QVariantMap();
909 feedback->setProgress( maxProgressDuringBlockWriting * i / mDataPerGroup.count() );
910 }
911 }
912
913 rasterDataProvider->setEditable( false );
914
915 if ( feedback )
916 feedback->setProgress( maxProgressDuringBlockWriting );
917
918 if ( feedback && hasReportsDuringClose )
919 {
920 std::unique_ptr<QgsFeedback> scaledFeedback( QgsFeedback::createScaledFeedback( feedback, maxProgressDuringBlockWriting, 100.0 ) );
921 if ( !rasterDataProvider->closeWithProgress( scaledFeedback.get() ) )
922 {
923 if ( feedback->isCanceled() )
924 return {};
925 throw QgsProcessingException( QObject::tr( "Could not write raster dataset" ) );
926 }
927 }
928
929 QVariantMap ret;
930 ret[u"OUTPUT"_s] = fileName;
931
932 return ret;
933}
934
935QSet<int> QgsMeshRasterizeAlgorithm::supportedDataType()
936{
938}
939
940QString QgsMeshContoursAlgorithm::name() const
941{
942 return u"meshcontours"_s;
943}
944
945QString QgsMeshContoursAlgorithm::displayName() const
946{
947 return QObject::tr( "Export contours" );
948}
949
950QString QgsMeshContoursAlgorithm::group() const
951{
952 return QObject::tr( "Mesh" );
953}
954
955QString QgsMeshContoursAlgorithm::groupId() const
956{
957 return u"mesh"_s;
958}
959
960QString QgsMeshContoursAlgorithm::shortHelpString() const
961{
962 return QObject::tr( "This algorithm creates contours as a vector layer from a mesh scalar dataset." );
963}
964
965QString QgsMeshContoursAlgorithm::shortDescription() const
966{
967 return QObject::tr( "Creates contours as vector layer from mesh scalar dataset." );
968}
969
970QgsProcessingAlgorithm *QgsMeshContoursAlgorithm::createInstance() const
971{
972 return new QgsMeshContoursAlgorithm();
973}
974
975void QgsMeshContoursAlgorithm::initAlgorithm( const QVariantMap &configuration )
976{
977 Q_UNUSED( configuration );
978
979 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
980
981 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
982
983 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
984
985 addParameter( new QgsProcessingParameterNumber( u"INCREMENT"_s, QObject::tr( "Increment between contour levels" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
986
987 addParameter( new QgsProcessingParameterNumber( u"MINIMUM"_s, QObject::tr( "Minimum contour level" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
988 addParameter( new QgsProcessingParameterNumber( u"MAXIMUM"_s, QObject::tr( "Maximum contour level" ), Qgis::ProcessingNumberParameterType::Double, QVariant(), true ) );
989
990 auto contourLevelList = std::make_unique<QgsProcessingParameterString>( u"CONTOUR_LEVEL_LIST"_s, QObject::tr( "List of contours level" ), QVariant(), false, true );
991 contourLevelList->setHelp( QObject::tr( "Comma separated list of values to export. If filled, the increment, minimum and maximum settings are ignored." ) );
992 addParameter( contourLevelList.release() );
993
994 addParameter( new QgsProcessingParameterCrs( u"CRS_OUTPUT"_s, QObject::tr( "Output coordinate system" ), QVariant(), true ) );
995
996
997 addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT_LINES"_s, QObject::tr( "Exported contour lines" ), Qgis::ProcessingSourceType::VectorLine ) );
998 addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT_POLYGONS"_s, QObject::tr( "Exported contour polygons" ), Qgis::ProcessingSourceType::VectorPolygon ) );
999}
1000
1001bool QgsMeshContoursAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1002{
1003 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
1004
1005 if ( !meshLayer || !meshLayer->isValid() )
1006 return false;
1007
1008 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
1009 if ( !outputCrs.isValid() )
1010 outputCrs = meshLayer->crs();
1011 mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() );
1012 if ( !meshLayer->nativeMesh() )
1013 meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh
1014
1015 mTriangularMesh.update( meshLayer->nativeMesh(), mTransform );
1016 mNativeMesh = *meshLayer->nativeMesh();
1017
1018 // Prepare levels
1019 mLevels.clear();
1020 // First, try with the levels list
1021 QString levelsString = parameterAsString( parameters, u"CONTOUR_LEVEL_LIST"_s, context );
1022 if ( !levelsString.isEmpty() )
1023 {
1024 QStringList levelStringList = levelsString.split( ',' );
1025 if ( !levelStringList.isEmpty() )
1026 {
1027 for ( const QString &stringVal : levelStringList )
1028 {
1029 bool ok;
1030 double val = stringVal.toDouble( &ok );
1031 if ( ok )
1032 mLevels.append( val );
1033 else
1034 throw QgsProcessingException( QObject::tr( "Invalid format for level values, must be numbers separated with comma" ) );
1035
1036 if ( mLevels.count() >= 2 )
1037 if ( mLevels.last() <= mLevels.at( mLevels.count() - 2 ) )
1038 throw QgsProcessingException( QObject::tr( "Invalid format for level values, must be different numbers and in increasing order" ) );
1039 }
1040 }
1041 }
1042
1043 if ( mLevels.isEmpty() )
1044 {
1045 double minimum = parameterAsDouble( parameters, u"MINIMUM"_s, context );
1046 double maximum = parameterAsDouble( parameters, u"MAXIMUM"_s, context );
1047 double interval = parameterAsDouble( parameters, u"INCREMENT"_s, context );
1048
1049 if ( interval <= 0 )
1050 throw QgsProcessingException( QObject::tr( "Invalid interval value, must be greater than zero" ) );
1051
1052 if ( minimum >= maximum )
1053 throw QgsProcessingException( QObject::tr( "Invalid minimum and maximum values, minimum must be lesser than maximum" ) );
1054
1055 if ( interval > ( maximum - minimum ) )
1056 throw QgsProcessingException( QObject::tr( "Invalid minimum, maximum and interval values, difference between minimum and maximum must be greater or equal than interval" ) );
1057
1058 int intervalCount = ( maximum - minimum ) / interval;
1059
1060 mLevels.reserve( intervalCount );
1061 for ( int i = 0; i < intervalCount; ++i )
1062 {
1063 mLevels.append( minimum + i * interval );
1064 }
1065 }
1066
1067 // Prepare data
1068 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
1069
1070 if ( feedback )
1071 {
1072 feedback->setProgressText( QObject::tr( "Preparing data" ) );
1073 }
1074
1075 // Extract the date time used to export dataset values under a relative time
1076 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
1077 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
1078
1079 mDateTimeString = meshLayer->formatTime( relativeTime.hours() );
1080
1081 extractDatasetValues( datasetGroups, meshLayer, mNativeMesh, relativeTime, supportedDataType(), mDataPerGroup, feedback );
1082
1083 mLayerRendererSettings = meshLayer->rendererSettings();
1084
1085 return true;
1086}
1087
1088QVariantMap QgsMeshContoursAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1089{
1090 //First, if present, average 3D staked dataset value to 2D face value
1091 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
1092 for ( DataGroup &dataGroup : mDataPerGroup )
1093 {
1094 if ( dataGroup.dataset3dStakedValue.isValid() )
1095 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
1096 }
1097
1098 // Create vector layers
1099 QgsFields polygonFields;
1100 QgsFields lineFields;
1101 polygonFields.append( QgsField( QObject::tr( "group" ), QMetaType::Type::QString ) );
1102 polygonFields.append( QgsField( QObject::tr( "time" ), QMetaType::Type::QString ) );
1103 polygonFields.append( QgsField( QObject::tr( "min_value" ), QMetaType::Type::Double ) );
1104 polygonFields.append( QgsField( QObject::tr( "max_value" ), QMetaType::Type::Double ) );
1105 lineFields.append( QgsField( QObject::tr( "group" ), QMetaType::Type::QString ) );
1106 lineFields.append( QgsField( QObject::tr( "time" ), QMetaType::Type::QString ) );
1107 lineFields.append( QgsField( QObject::tr( "value" ), QMetaType::Type::Double ) );
1108
1109 QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, u"CRS_OUTPUT"_s, context );
1110
1111 QString lineIdentifier;
1112 QString polygonIdentifier;
1113 std::unique_ptr<QgsFeatureSink> sinkPolygons( parameterAsSink( parameters, u"OUTPUT_POLYGONS"_s, context, polygonIdentifier, polygonFields, Qgis::WkbType::PolygonZ, outputCrs ) );
1114 std::unique_ptr<QgsFeatureSink> sinkLines( parameterAsSink( parameters, u"OUTPUT_LINES"_s, context, lineIdentifier, lineFields, Qgis::WkbType::LineStringZ, outputCrs ) );
1115
1116 if ( !sinkLines || !sinkPolygons )
1117 return QVariantMap();
1118
1119
1120 for ( int i = 0; i < mDataPerGroup.count(); ++i )
1121 {
1122 DataGroup dataGroup = mDataPerGroup.at( i );
1123 bool scalarDataOnVertices = dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
1124 int count = scalarDataOnVertices ? mNativeMesh.vertices.count() : mNativeMesh.faces.count();
1125
1126 QVector<double> values;
1127 if ( dataGroup.datasetValues.isValid() )
1128 {
1129 // vals could be scalar or vectors, for contour rendering we want always magnitude
1130 values = QgsMeshLayerUtils::calculateMagnitudes( dataGroup.datasetValues );
1131 }
1132 else
1133 {
1134 values = QVector<double>( count, std::numeric_limits<double>::quiet_NaN() );
1135 }
1136
1137 if ( ( !scalarDataOnVertices ) )
1138 {
1139 values = QgsMeshLayerUtils::interpolateFromFacesData( values, mNativeMesh, &dataGroup.activeFaces, QgsMeshRendererScalarSettings::NeighbourAverage );
1140 }
1141
1142 QgsMeshContours contoursExported( mTriangularMesh, mNativeMesh, values, dataGroup.activeFaces );
1143
1144 QgsAttributes firstAttributes;
1145 firstAttributes.append( dataGroup.metadata.name() );
1146 firstAttributes.append( mDateTimeString );
1147
1148 for ( double level : std::as_const( mLevels ) )
1149 {
1150 QgsGeometry line = contoursExported.exportLines( level, feedback );
1151 if ( feedback->isCanceled() )
1152 return QVariantMap();
1153 if ( line.isEmpty() )
1154 continue;
1155 QgsAttributes lineAttributes = firstAttributes;
1156 lineAttributes.append( level );
1157
1158 QgsFeature lineFeat;
1159 lineFeat.setGeometry( line );
1160 lineFeat.setAttributes( lineAttributes );
1161
1162 if ( !sinkLines->addFeature( lineFeat, QgsFeatureSink::FastInsert ) )
1163 throw QgsProcessingException( writeFeatureError( sinkLines.get(), parameters, u"OUTPUT_LINES"_s ) );
1164 else
1165 feedback->featureAddedToSink( u"OUTPUT_LINES"_s );
1166 }
1167
1168 for ( int l = 0; l < mLevels.count() - 1; ++l )
1169 {
1170 QgsGeometry polygon = contoursExported.exportPolygons( mLevels.at( l ), mLevels.at( l + 1 ), feedback );
1171 if ( feedback->isCanceled() )
1172 return QVariantMap();
1173
1174 if ( polygon.isEmpty() )
1175 continue;
1176 QgsAttributes polygonAttributes = firstAttributes;
1177 polygonAttributes.append( mLevels.at( l ) );
1178 polygonAttributes.append( mLevels.at( l + 1 ) );
1179
1180 QgsFeature polygonFeature;
1181 polygonFeature.setGeometry( polygon );
1182 polygonFeature.setAttributes( polygonAttributes );
1183 if ( !sinkPolygons->addFeature( polygonFeature ) )
1184 {
1185 throw QgsProcessingException( writeFeatureError( sinkPolygons.get(), parameters, QString() ) );
1186 }
1187 else
1188 {
1189 feedback->featureAddedToSink( u"OUTPUT_POLYGONS"_s );
1190 }
1191 }
1192
1193 if ( feedback )
1194 {
1195 feedback->setProgress( 100 * i / mDataPerGroup.count() );
1196 }
1197 }
1198
1199 if ( sinkPolygons )
1200 {
1201 sinkPolygons->finalize();
1202 feedback->featureSinkFinalized( u"OUTPUT_POLYGONS"_s );
1203 }
1204 if ( sinkLines )
1205 {
1206 sinkLines->finalize();
1207 feedback->featureSinkFinalized( u"OUTPUT_LINES"_s );
1208 }
1209
1210 QVariantMap ret;
1211 ret[u"OUTPUT_LINES"_s] = lineIdentifier;
1212 ret[u"OUTPUT_POLYGONS"_s] = polygonIdentifier;
1213
1214 return ret;
1215}
1216
1217QString QgsMeshExportCrossSection::name() const
1218{
1219 return u"meshexportcrosssection"_s;
1220}
1221
1222QString QgsMeshExportCrossSection::displayName() const
1223{
1224 return QObject::tr( "Export cross section dataset values on lines from mesh" );
1225}
1226
1227QString QgsMeshExportCrossSection::group() const
1228{
1229 return QObject::tr( "Mesh" );
1230}
1231
1232QString QgsMeshExportCrossSection::groupId() const
1233{
1234 return u"mesh"_s;
1235}
1236
1237QString QgsMeshExportCrossSection::shortHelpString() const
1238{
1239 return QObject::tr(
1240 "This algorithm extracts mesh's dataset values from line contained in a vector layer.\n"
1241 "Each line is discretized with a resolution distance parameter for extraction of values on its vertices."
1242 );
1243}
1244
1245QString QgsMeshExportCrossSection::shortDescription() const
1246{
1247 return QObject::tr( "Extracts a mesh dataset's values from lines contained in a vector layer." );
1248}
1249
1250QgsProcessingAlgorithm *QgsMeshExportCrossSection::createInstance() const
1251{
1252 return new QgsMeshExportCrossSection();
1253}
1254
1255void QgsMeshExportCrossSection::initAlgorithm( const QVariantMap &configuration )
1256{
1257 Q_UNUSED( configuration );
1258
1259 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
1260
1261 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
1262
1263 addParameter( new QgsProcessingParameterMeshDatasetTime( u"DATASET_TIME"_s, QObject::tr( "Dataset time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
1264
1265 QList<int> datatype;
1266 datatype << static_cast<int>( Qgis::ProcessingSourceType::VectorLine );
1267 addParameter( new QgsProcessingParameterFeatureSource( u"INPUT_LINES"_s, QObject::tr( "Lines for data export" ), datatype, QVariant(), false ) );
1268
1269 addParameter( new QgsProcessingParameterDistance( u"RESOLUTION"_s, QObject::tr( "Line segmentation resolution" ), 10.0, u"INPUT_LINES"_s, false, 0 ) );
1270
1271 addParameter( new QgsProcessingParameterNumber( u"COORDINATES_DIGITS"_s, QObject::tr( "Digits count for coordinates" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1272
1273 addParameter( new QgsProcessingParameterNumber( u"DATASET_DIGITS"_s, QObject::tr( "Digits count for dataset value" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1274
1275 addParameter( new QgsProcessingParameterFileDestination( u"OUTPUT"_s, QObject::tr( "Exported data CSV file" ), QObject::tr( "CSV file (*.csv)" ) ) );
1276}
1277
1278bool QgsMeshExportCrossSection::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1279{
1280 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
1281
1282 if ( !meshLayer || !meshLayer->isValid() )
1283 return false;
1284
1285 mMeshLayerCrs = meshLayer->crs();
1286 mTriangularMesh.update( meshLayer->nativeMesh() );
1287 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
1288
1289 if ( feedback )
1290 {
1291 feedback->setProgressText( QObject::tr( "Preparing data" ) );
1292 }
1293
1294 // Extract the date time used to export dataset values under a relative time
1295 QVariant parameterTimeVariant = parameters.value( u"DATASET_TIME"_s );
1296 QgsInterval relativeTime = datasetRelativetime( parameterTimeVariant, meshLayer, context );
1297
1298 extractDatasetValues( datasetGroups, meshLayer, *meshLayer->nativeMesh(), relativeTime, supportedDataType(), mDataPerGroup, feedback );
1299
1300 mLayerRendererSettings = meshLayer->rendererSettings();
1301
1302 return true;
1303}
1304
1305QVariantMap QgsMeshExportCrossSection::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1306{
1307 if ( feedback )
1308 feedback->setProgress( 0 );
1309 //First, if present, average 3D staked dataset value to 2D face value
1310 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
1311 for ( DataGroup &dataGroup : mDataPerGroup )
1312 {
1313 if ( dataGroup.dataset3dStakedValue.isValid() )
1314 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
1315 }
1316 double resolution = parameterAsDouble( parameters, u"RESOLUTION"_s, context );
1317 int datasetDigits = parameterAsInt( parameters, u"DATASET_DIGITS"_s, context );
1318 int coordDigits = parameterAsInt( parameters, u"COORDINATES_DIGITS"_s, context );
1319
1320 std::unique_ptr<QgsProcessingFeatureSource> featureSource( parameterAsSource( parameters, u"INPUT_LINES"_s, context ) );
1321 if ( !featureSource )
1322 throw QgsProcessingException( QObject::tr( "Input lines vector layer required" ) );
1323
1324 QgsCoordinateTransform transform( featureSource->sourceCrs(), mMeshLayerCrs, context.transformContext() );
1325
1326 QString outputFileName = parameterAsFileOutput( parameters, u"OUTPUT"_s, context );
1327 QFile file( outputFileName );
1328 if ( !file.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
1329 throw QgsProcessingException( QObject::tr( "Unable to create the output file" ) );
1330
1331 QTextStream textStream( &file );
1332 QStringList header;
1333 header << u"fid"_s << u"x"_s << u"y"_s << QObject::tr( "offset" );
1334 for ( const DataGroup &datagroup : std::as_const( mDataPerGroup ) )
1335 header << datagroup.metadata.name();
1336 textStream << header.join( ',' ) << u"\n"_s;
1337
1338 long long featCount = featureSource->featureCount();
1339 long long featCounter = 0;
1340 QgsFeatureIterator featIt = featureSource->getFeatures();
1341 QgsFeature feat;
1342 while ( featIt.nextFeature( feat ) )
1343 {
1344 QgsFeatureId fid = feat.id();
1345 QgsGeometry line = feat.geometry();
1346 try
1347 {
1348 line.transform( transform );
1349 }
1350 catch ( QgsCsException & )
1351 {
1352 line = feat.geometry();
1353 feedback->reportError( QObject::tr( "Could not transform line to mesh CRS" ) );
1354 }
1355
1356 if ( line.isEmpty() )
1357 continue;
1358 double offset = 0;
1359 while ( offset <= line.length() )
1360 {
1361 if ( feedback->isCanceled() )
1362 return QVariantMap();
1363
1364 QStringList textLine;
1365 QgsPointXY point = line.interpolate( offset ).asPoint();
1366 int triangularFaceIndex = mTriangularMesh.faceIndexForPoint_v2( point );
1367 textLine << QString::number( fid ) << QString::number( point.x(), 'f', coordDigits ) << QString::number( point.y(), 'f', coordDigits ) << QString::number( offset, 'f', coordDigits );
1368 if ( triangularFaceIndex >= 0 )
1369 {
1370 //extract dataset values for the point
1371 QgsAttributes attributes;
1372 int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces().at( triangularFaceIndex );
1373 for ( int i = 0; i < mDataPerGroup.count(); ++i )
1374 {
1375 const DataGroup &dataGroup = mDataPerGroup.at( i );
1376 bool faceActive = dataGroup.activeFaces.active( nativeFaceIndex );
1377 if ( !faceActive )
1378 continue;
1379 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
1380
1381 if ( abs( value.x() ) == std::numeric_limits<double>::quiet_NaN() )
1382 textLine << QString( ' ' );
1383 else
1384 textLine << QString::number( value.scalar(), 'f', datasetDigits );
1385 }
1386 }
1387 else
1388 for ( int i = 0; i < mDataPerGroup.count(); ++i )
1389 textLine << QString( ' ' );
1390
1391 textStream << textLine.join( ',' ) << u"\n"_s;
1392
1393 offset += resolution;
1394 }
1395
1396 if ( feedback )
1397 {
1398 feedback->setProgress( 100.0 * featCounter / featCount );
1399 if ( feedback->isCanceled() )
1400 return QVariantMap();
1401 }
1402 }
1403
1404 file.close();
1405
1406 QVariantMap ret;
1407 ret[u"OUTPUT"_s] = outputFileName;
1408 return ret;
1409}
1410
1411QString QgsMeshExportTimeSeries::name() const
1412{
1413 return u"meshexporttimeseries"_s;
1414}
1415
1416QString QgsMeshExportTimeSeries::displayName() const
1417{
1418 return QObject::tr( "Export time series values from points of a mesh dataset" );
1419}
1420
1421QString QgsMeshExportTimeSeries::group() const
1422{
1423 return QObject::tr( "Mesh" );
1424}
1425
1426QString QgsMeshExportTimeSeries::groupId() const
1427{
1428 return u"mesh"_s;
1429}
1430
1431QString QgsMeshExportTimeSeries::shortHelpString() const
1432{
1433 return QObject::tr(
1434 "This algorithm extracts mesh's dataset time series values from points contained in a vector layer.\n"
1435 "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."
1436 );
1437}
1438
1439QString QgsMeshExportTimeSeries::shortDescription() const
1440{
1441 return QObject::tr( "Extracts a mesh dataset's time series values from points contained in a vector layer." );
1442}
1443
1444QgsProcessingAlgorithm *QgsMeshExportTimeSeries::createInstance() const
1445{
1446 return new QgsMeshExportTimeSeries();
1447}
1448
1449void QgsMeshExportTimeSeries::initAlgorithm( const QVariantMap &configuration )
1450{
1451 Q_UNUSED( configuration );
1452
1453 addParameter( new QgsProcessingParameterMeshLayer( u"INPUT"_s, QObject::tr( "Input mesh layer" ) ) );
1454
1455 addParameter( new QgsProcessingParameterMeshDatasetGroups( u"DATASET_GROUPS"_s, QObject::tr( "Dataset groups" ), u"INPUT"_s, supportedDataType() ) );
1456
1457 addParameter( new QgsProcessingParameterMeshDatasetTime( u"STARTING_TIME"_s, QObject::tr( "Starting time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
1458
1459 addParameter( new QgsProcessingParameterMeshDatasetTime( u"FINISHING_TIME"_s, QObject::tr( "Finishing time" ), u"INPUT"_s, u"DATASET_GROUPS"_s ) );
1460
1461 addParameter( new QgsProcessingParameterNumber( u"TIME_STEP"_s, QObject::tr( "Time step (hours)" ), Qgis::ProcessingNumberParameterType::Double, 0, true, 0 ) );
1462
1463 QList<int> datatype;
1464 datatype << static_cast<int>( Qgis::ProcessingSourceType::VectorPoint );
1465 addParameter( new QgsProcessingParameterFeatureSource( u"INPUT_POINTS"_s, QObject::tr( "Points for data export" ), datatype, QVariant(), false ) );
1466
1467 addParameter( new QgsProcessingParameterNumber( u"COORDINATES_DIGITS"_s, QObject::tr( "Digits count for coordinates" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1468
1469 addParameter( new QgsProcessingParameterNumber( u"DATASET_DIGITS"_s, QObject::tr( "Digits count for dataset value" ), Qgis::ProcessingNumberParameterType::Integer, 2 ) );
1470
1471 addParameter( new QgsProcessingParameterFileDestination( u"OUTPUT"_s, QObject::tr( "Exported data CSV file" ), QObject::tr( "CSV file (*.csv)" ) ) );
1472}
1473
1474bool QgsMeshExportTimeSeries::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1475{
1476 QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, u"INPUT"_s, context );
1477
1478 if ( !meshLayer || !meshLayer->isValid() )
1479 return false;
1480
1481 mMeshLayerCrs = meshLayer->crs();
1482 mTriangularMesh.update( meshLayer->nativeMesh() );
1483
1484 QList<int> datasetGroups = QgsProcessingParameterMeshDatasetGroups::valueAsDatasetGroup( parameters.value( u"DATASET_GROUPS"_s ) );
1485
1486 if ( feedback )
1487 {
1488 feedback->setProgressText( QObject::tr( "Preparing data" ) );
1489 }
1490
1491 // Extract the date times used to export dataset values
1492 QVariant parameterStartTimeVariant = parameters.value( u"STARTING_TIME"_s );
1493 QgsInterval relativeStartTime = datasetRelativetime( parameterStartTimeVariant, meshLayer, context );
1494
1495 QVariant parameterEndTimeVariant = parameters.value( u"FINISHING_TIME"_s );
1496 QgsInterval relativeEndTime = datasetRelativetime( parameterEndTimeVariant, meshLayer, context );
1497
1498 // calculate time steps
1499 qint64 timeStepInterval = parameterAsDouble( parameters, u"TIME_STEP"_s, context ) * 1000 * 3600;
1500 if ( timeStepInterval == 0 )
1501 {
1502 //take the first time step of the first temporal dataset group
1503 for ( int groupIndex : datasetGroups )
1504 {
1505 QgsMeshDatasetGroupMetadata meta = meshLayer->datasetGroupMetadata( QgsMeshDatasetIndex( groupIndex, 0 ) );
1506 if ( !meta.isTemporal() && meshLayer->datasetCount( QgsMeshDatasetIndex( groupIndex, 0 ) ) < 2 )
1507 continue;
1508 else
1509 {
1510 timeStepInterval = meshLayer->datasetRelativeTimeInMilliseconds( QgsMeshDatasetIndex( groupIndex, 1 ) ) - meshLayer->datasetRelativeTimeInMilliseconds( QgsMeshDatasetIndex( groupIndex, 0 ) );
1511 break;
1512 }
1513 }
1514 }
1515
1516 mRelativeTimeSteps.clear();
1517 mTimeStepString.clear();
1518 if ( timeStepInterval != 0 )
1519 {
1520 mRelativeTimeSteps.append( relativeStartTime.seconds() * 1000 );
1521 while ( mRelativeTimeSteps.last() < relativeEndTime.seconds() * 1000 )
1522 mRelativeTimeSteps.append( mRelativeTimeSteps.last() + timeStepInterval );
1523
1524 for ( qint64 relativeTimeStep : std::as_const( mRelativeTimeSteps ) )
1525 {
1526 mTimeStepString.append( meshLayer->formatTime( relativeTimeStep / 3600.0 / 1000.0 ) );
1527 }
1528 }
1529
1530 //Extract needed dataset values
1531 for ( int i = 0; i < datasetGroups.count(); ++i )
1532 {
1533 int groupIndex = datasetGroups.at( i );
1534 QgsMeshDatasetGroupMetadata meta = meshLayer->datasetGroupMetadata( QgsMeshDatasetIndex( groupIndex, 0 ) );
1535 if ( supportedDataType().contains( meta.dataType() ) )
1536 {
1537 mGroupIndexes.append( groupIndex );
1538 mGroupsMetadata[groupIndex] = meta;
1539 int valueCount = meta.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices ? mTriangularMesh.vertices().count() : meshLayer->nativeMesh()->faceCount();
1540
1541 if ( !mRelativeTimeSteps.isEmpty() )
1542 {
1543 //QMap<qint64, DataGroup> temporalGroup;
1544 QgsMeshDatasetIndex lastDatasetIndex;
1545 for ( qint64 relativeTimeStep : std::as_const( mRelativeTimeSteps ) )
1546 {
1547 QMap<int, int> &groupIndexToData = mRelativeTimeToData[relativeTimeStep];
1548 QgsInterval timeStepInterval( relativeTimeStep / 1000.0 );
1549 QgsMeshDatasetIndex datasetIndex = meshLayer->datasetIndexAtRelativeTime( timeStepInterval, groupIndex );
1550 if ( !datasetIndex.isValid() )
1551 continue;
1552 if ( datasetIndex != lastDatasetIndex )
1553 {
1554 DataGroup dataGroup;
1555 dataGroup.metadata = meta;
1556 dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, valueCount );
1557 dataGroup.activeFaces = meshLayer->areFacesActive( datasetIndex, 0, meshLayer->nativeMesh()->faceCount() );
1558 if ( dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVolumes )
1559 {
1560 dataGroup.dataset3dStakedValue = meshLayer->dataset3dValues( datasetIndex, 0, valueCount );
1561 }
1562 mDatasets.append( dataGroup );
1563 lastDatasetIndex = datasetIndex;
1564 }
1565 groupIndexToData[groupIndex] = mDatasets.count() - 1;
1566 }
1567 }
1568 else
1569 {
1570 // we have only static dataset group
1571 QMap<int, int> &groupIndexToData = mRelativeTimeToData[0];
1572 QgsMeshDatasetIndex datasetIndex( groupIndex, 0 );
1573 DataGroup dataGroup;
1574 dataGroup.metadata = meta;
1575 dataGroup.datasetValues = meshLayer->datasetValues( datasetIndex, 0, valueCount );
1576 dataGroup.activeFaces = meshLayer->areFacesActive( datasetIndex, 0, meshLayer->nativeMesh()->faceCount() );
1577 if ( dataGroup.metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVolumes )
1578 {
1579 dataGroup.dataset3dStakedValue = meshLayer->dataset3dValues( datasetIndex, 0, valueCount );
1580 }
1581 mDatasets.append( dataGroup );
1582 groupIndexToData[groupIndex] = mDatasets.count() - 1;
1583 }
1584 }
1585
1586 if ( feedback )
1587 feedback->setProgress( 100 * i / datasetGroups.count() );
1588 }
1589
1590 mLayerRendererSettings = meshLayer->rendererSettings();
1591
1592 return true;
1593}
1594
1595
1596QVariantMap QgsMeshExportTimeSeries::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
1597{
1598 if ( feedback )
1599 feedback->setProgress( 0 );
1600 //First, if present, average 3D staked dataset value to 2D face value
1601 const QgsMesh3DAveragingMethod *avgMethod = mLayerRendererSettings.averagingMethod();
1602
1603 for ( DataGroup &dataGroup : mDatasets )
1604 {
1605 if ( dataGroup.dataset3dStakedValue.isValid() )
1606 dataGroup.datasetValues = avgMethod->calculate( dataGroup.dataset3dStakedValue );
1607 }
1608
1609 int datasetDigits = parameterAsInt( parameters, u"DATASET_DIGITS"_s, context );
1610 int coordDigits = parameterAsInt( parameters, u"COORDINATES_DIGITS"_s, context );
1611
1612 std::unique_ptr<QgsProcessingFeatureSource> featureSource( parameterAsSource( parameters, u"INPUT_POINTS"_s, context ) );
1613 if ( !featureSource )
1614 throw QgsProcessingException( QObject::tr( "Input points vector layer required" ) );
1615
1616 QgsCoordinateTransform transform( featureSource->sourceCrs(), mMeshLayerCrs, context.transformContext() );
1617
1618 QString outputFileName = parameterAsFileOutput( parameters, u"OUTPUT"_s, context );
1619 QFile file( outputFileName );
1620 if ( !file.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
1621 throw QgsProcessingException( QObject::tr( "Unable to create the output file" ) );
1622
1623 QTextStream textStream( &file );
1624 QStringList header;
1625 header << u"fid"_s << u"x"_s << u"y"_s << QObject::tr( "time" );
1626
1627 for ( int gi : std::as_const( mGroupIndexes ) )
1628 header << mGroupsMetadata.value( gi ).name();
1629
1630 textStream << header.join( ',' ) << u"\n"_s;
1631
1632 long long featCount = featureSource->featureCount();
1633 long long featCounter = 0;
1634 QgsFeatureIterator featIt = featureSource->getFeatures();
1635 QgsFeature feat;
1636 while ( featIt.nextFeature( feat ) )
1637 {
1638 QgsFeatureId fid = feat.id();
1639 QgsGeometry geom = feat.geometry();
1640 try
1641 {
1642 geom.transform( transform );
1643 }
1644 catch ( QgsCsException & )
1645 {
1646 geom = feat.geometry();
1647 feedback->reportError( QObject::tr( "Could not transform line to mesh CRS" ) );
1648 }
1649
1650 if ( geom.isEmpty() )
1651 continue;
1652
1653 QgsPointXY point = geom.asPoint();
1654 int triangularFaceIndex = mTriangularMesh.faceIndexForPoint_v2( point );
1655
1656 if ( triangularFaceIndex >= 0 )
1657 {
1658 int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces().at( triangularFaceIndex );
1659 if ( !mRelativeTimeSteps.isEmpty() )
1660 {
1661 for ( int timeIndex = 0; timeIndex < mRelativeTimeSteps.count(); ++timeIndex )
1662 {
1663 qint64 timeStep = mRelativeTimeSteps.at( timeIndex );
1664 QStringList textLine;
1665 textLine << QString::number( fid ) << QString::number( point.x(), 'f', coordDigits ) << QString::number( point.y(), 'f', coordDigits ) << mTimeStepString.at( timeIndex );
1666
1667 if ( mRelativeTimeToData.contains( timeStep ) )
1668 {
1669 const QMap<int, int> &groupToData = mRelativeTimeToData.value( timeStep );
1670 for ( int groupIndex : std::as_const( mGroupIndexes ) )
1671 {
1672 if ( !groupToData.contains( groupIndex ) )
1673 continue;
1674 int dataIndex = groupToData.value( groupIndex );
1675 if ( dataIndex < 0 || dataIndex > mDatasets.count() - 1 )
1676 continue;
1677
1678 const DataGroup &dataGroup = mDatasets.at( dataIndex );
1679 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
1680 if ( abs( value.x() ) == std::numeric_limits<double>::quiet_NaN() )
1681 textLine << QString( ' ' );
1682 else
1683 textLine << QString::number( value.scalar(), 'f', datasetDigits );
1684 }
1685 }
1686 textStream << textLine.join( ',' ) << u"\n"_s;
1687 }
1688 }
1689 else
1690 {
1691 QStringList textLine;
1692 textLine << QString::number( fid ) << QString::number( point.x(), 'f', coordDigits ) << QString::number( point.y(), 'f', coordDigits ) << QObject::tr( "static dataset" );
1693 const QMap<int, int> &groupToData = mRelativeTimeToData.value( 0 );
1694 for ( int groupIndex : std::as_const( mGroupIndexes ) )
1695 {
1696 if ( !groupToData.contains( groupIndex ) )
1697 continue;
1698 int dataIndex = groupToData.value( groupIndex );
1699 if ( dataIndex < 0 || dataIndex > mDatasets.count() - 1 )
1700 continue;
1701 const DataGroup &dataGroup = mDatasets.at( dataIndex );
1702 QgsMeshDatasetValue value = extractDatasetValue( point, nativeFaceIndex, triangularFaceIndex, mTriangularMesh, dataGroup.activeFaces, dataGroup.datasetValues, dataGroup.metadata );
1703 if ( abs( value.x() ) == std::numeric_limits<double>::quiet_NaN() )
1704 textLine << QString( ' ' );
1705 else
1706 textLine << QString::number( value.scalar(), 'f', datasetDigits );
1707 }
1708 textStream << textLine.join( ',' ) << u"\n"_s;
1709 }
1710 }
1711 featCounter++;
1712 if ( feedback )
1713 {
1714 feedback->setProgress( 100.0 * featCounter / featCount );
1715 if ( feedback->isCanceled() )
1716 return QVariantMap();
1717 }
1718 }
1719
1720 file.close();
1721
1722 QVariantMap ret;
1723 ret[u"OUTPUT"_s] = outputFileName;
1724 return ret;
1725}
1726
@ VectorPoint
Vector point layers.
Definition qgis.h:3729
@ VectorPolygon
Vector polygon layers.
Definition qgis.h:3731
@ VectorLine
Vector line layers.
Definition qgis.h:3730
@ 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:3962
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
Definition qgis.h:3961
@ Double
Double/float values.
Definition qgis.h:4002
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:7401
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.