diff --git a/python/plugins/processing/tests/AlgorithmsTestBase.py b/python/plugins/processing/tests/AlgorithmsTestBase.py index e2255a0e32aa..08637ac8186d 100644 --- a/python/plugins/processing/tests/AlgorithmsTestBase.py +++ b/python/plugins/processing/tests/AlgorithmsTestBase.py @@ -42,6 +42,7 @@ QgsCoordinateReferenceSystem, QgsFeatureRequest, QgsMapLayer, + QgsMeshLayer, QgsProject, QgsApplication, QgsProcessingContext, @@ -267,7 +268,7 @@ def load_param(self, param, id=None): parameter based on its key `type` and return the appropriate parameter to pass to the algorithm. """ try: - if param["type"] in ("vector", "raster", "table"): + if param["type"] in ("vector", "raster", "table", "mesh"): return self.load_layer(id, param).id() elif param["type"] == "vrtlayers": vals = [] @@ -377,6 +378,8 @@ def load_layer(self, id, param): options = QgsRasterLayer.LayerOptions() options.loadDefaultStyle = False lyr = QgsRasterLayer(filepath, param["name"], "gdal", options) + elif param["type"] == "mesh": + lyr = QgsMeshLayer(filepath, param["name"], "mdal") self.assertTrue( lyr.isValid(), f'Could not load layer "{filepath}" from param {param}' diff --git a/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon.gml b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon.gml new file mode 100644 index 000000000000..b15b731e4923 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon.gml @@ -0,0 +1,16 @@ + + + 1000 20003000 3000 + + + + 1000 20003000 3000 + 1000 2000 2000 2000 3000 2000 2000 3000 1000 3000 1000 2000 + + + diff --git a/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon.xsd b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon.xsd new file mode 100644 index 000000000000..21f8960b3e6f --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon.xsd @@ -0,0 +1,47 @@ + + + + + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon_complex_mesh_hole.gml b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon_complex_mesh_hole.gml new file mode 100644 index 000000000000..4653ded25f57 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon_complex_mesh_hole.gml @@ -0,0 +1,16 @@ + + + -0.28285357 -10.3031914941.54679257 15.32259074 + + + + -0.28285357 -10.3031914941.54679257 15.32259074 + -0.28285357 -5.15175219 0 0 0.0 7.5 0 10 5.0 9.5 5.10513141 15.27002503 22.95118899 15.32259074 28.07634543 15.24374218 26.65707134 9.80319149 25.57947434 -1.10419274 21.66332916 -5.91395494 15.80225282 -10.30319149 10.25657071 -10.25062578 -0.28285357 -5.1517521910 0 10 10 15.53942428 9.96088861 16.42244099 7.57792771 15.01376721 -0.18429287 10 032.85840157 5.894131 37.37097674 6.83705716 41.54679257 12.62931782 34.20543893 11.48433606 32.85840157 5.894131 + + + diff --git a/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon_complex_mesh_hole.xsd b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon_complex_mesh_hole.xsd new file mode 100644 index 000000000000..58c91b2bda13 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/mesh_surface_to_polygon_complex_mesh_hole.xsd @@ -0,0 +1,47 @@ + + + + + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/python/plugins/processing/tests/testdata/mesh_two_part_with_hole.2dm b/python/plugins/processing/tests/testdata/mesh_two_part_with_hole.2dm new file mode 100644 index 000000000000..238373d02c20 --- /dev/null +++ b/python/plugins/processing/tests/testdata/mesh_two_part_with_hole.2dm @@ -0,0 +1,68 @@ +MESH2D +ND 1 0 0 0 +ND 2 10 0 0 +ND 3 10 10 0 +ND 4 5 4 0 +ND 5 5 2 0 +ND 6 0 7.5 0 +ND 7 0 10 0 +ND 8 5 9.5 0 +ND 9 5 6.75 0 +ND 10 15.53942428 9.96088861 0 +ND 11 15.01376721 -0.18429287 0 +ND 12 12.38548185 14.61295369 0 +ND 13 22.29411765 9.8557572 0 +ND 14 5.10513141 15.27002503 0 +ND 15 22.95118899 15.32259074 0 +ND 16 10.17772215 -5.23060075 0 +ND 17 -0.28285357 -5.15175219 0 +ND 18 14.90863579 -5.41458073 0 +ND 19 10.25657071 -10.25062578 0 +ND 20 15.80225282 -10.30319149 0 +ND 21 21.47934919 -1.07790989 0 +ND 22 21.66332916 -5.91395494 0 +ND 23 26.65707134 9.80319149 0 +ND 24 25.57947434 -1.10419274 0 +ND 25 28.07634543 15.24374218 0 +ND 26 16.42244099 7.57792771 0 +ND 27 15.93402616 -0.31148313 0 +ND 28 32.85840157 5.894131 0 +ND 29 34.20543893 11.48433606 0 +ND 30 37.37097674 6.83705716 0 +ND 31 41.54679257 12.62931782 0 +E3T 1 3 8 2 +E3T 2 8 9 2 +E3T 3 9 4 2 +E3T 4 4 5 2 +E3T 5 5 1 2 +E3T 6 5 6 1 +E3T 7 6 5 7 +E3T 8 10 12 3 +E3T 9 12 8 3 +E3T 10 13 12 10 +E3T 11 12 14 8 +E3T 12 13 15 12 +E3T 13 12 15 14 +E3T 14 7 5 4 +E3T 15 9 7 4 +E3T 16 9 8 7 +E3T 17 1 16 2 +E3T 18 1 17 16 +E4Q 19 16 18 11 2 +E3T 20 17 19 16 +E4Q 21 19 20 18 16 +E3T 22 18 22 21 +E3T 23 20 22 18 +E3T 24 13 23 15 +E3T 25 13 21 23 +E3T 26 21 24 23 +E3T 27 21 22 24 +E3T 28 23 25 15 +E3T 29 26 11 13 +E3T 30 26 13 10 +E3T 31 11 18 27 +E3T 32 27 18 21 +E3T 33 11 27 13 +E3T 34 27 21 13 +E3T 35 28 30 29 +E3T 36 30 31 29 diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests5.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests5.yaml index 00b79ac80025..0bd255e05a92 100644 --- a/python/plugins/processing/tests/testdata/qgis_algorithm_tests5.yaml +++ b/python/plugins/processing/tests/testdata/qgis_algorithm_tests5.yaml @@ -247,6 +247,28 @@ tests: name: expected/create_grid_negative_overlay.gml type: vector + - algorithm: native:surfacetopolygon + name: Create polygon from surface of Mesh layer + params: + INPUT: + name: quad_and_triangle.2dm + type: mesh + results: + OUTPUT: + name: expected/mesh_surface_to_polygon.gml + type: vector + + - algorithm: native:surfacetopolygon + name: Create polygon from surface of Mesh layer - multiple parts with hole + params: + INPUT: + name: mesh_two_part_with_hole.2dm + type: mesh + results: + OUTPUT: + name: expected/mesh_surface_to_polygon_complex_mesh_hole.gml + type: vector + - algorithm: native:rasterminmax name: Raster min max params: diff --git a/python/plugins/processing/tests/testdata/quad_and_triangle.2dm b/python/plugins/processing/tests/testdata/quad_and_triangle.2dm new file mode 100644 index 000000000000..01a1aeaa95b0 --- /dev/null +++ b/python/plugins/processing/tests/testdata/quad_and_triangle.2dm @@ -0,0 +1,8 @@ +MESH2D 1000.000 2000.000 0.000000 200 300 1.000 1.000 +ND 1 1000.000 2000.000 20.000 +ND 2 2000.000 2000.000 30.000 +ND 3 3000.000 2000.000 40.000 +ND 4 2000.000 3000.000 50.000 +ND 5 1000.000 3000.000 10.000 +E4Q 1 1 2 4 5 1 +E3T 2 2 3 4 1 diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index b4739fe4cb65..cb560fb36d50 100644 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -144,6 +144,7 @@ set(QGIS_ANALYSIS_SRCS processing/qgsalgorithmloadlayer.cpp processing/qgsalgorithmmeancoordinates.cpp processing/qgsalgorithmmergelines.cpp + processing/qgsalgorithmmeshsurfacetopolygon.cpp processing/qgsalgorithmmergevector.cpp processing/qgsalgorithmminimumenclosingcircle.cpp processing/qgsalgorithmmultidifference.cpp diff --git a/src/analysis/processing/qgsalgorithmmeshsurfacetopolygon.cpp b/src/analysis/processing/qgsalgorithmmeshsurfacetopolygon.cpp new file mode 100644 index 000000000000..7f1d0a1b8589 --- /dev/null +++ b/src/analysis/processing/qgsalgorithmmeshsurfacetopolygon.cpp @@ -0,0 +1,312 @@ +/*************************************************************************** + qgsalgorithmmeshsurfacetopolygon.cpp + --------------------------- + begin : September 2024 + copyright : (C) 2024 by Jan Caha + email : jan.caha at outlook dot com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include "qgsalgorithmmeshsurfacetopolygon.h" +#include "qgsprocessingparametermeshdataset.h" +#include "qgsmeshlayer.h" +#include "qgspolygon.h" +#include "qgslinestring.h" +#include "qgsmultilinestring.h" +#include "qgsmultipolygon.h" +#include "qgsgeometryengine.h" + +#include + +///@cond PRIVATE + + +QString QgsMeshSurfaceToPolygonAlgorithm::shortHelpString() const +{ + return QObject::tr( "This algorithm exports a polygon file containing mesh layer boundary. It may contain holes and it may be a multi-part polygon." ); +} + +QString QgsMeshSurfaceToPolygonAlgorithm::name() const +{ + return QStringLiteral( "surfacetopolygon" ); +} + +QString QgsMeshSurfaceToPolygonAlgorithm::displayName() const +{ + return QObject::tr( "Surface to polygon" ); +} + +QString QgsMeshSurfaceToPolygonAlgorithm::group() const +{ + return QObject::tr( "Mesh" ); +} + +QString QgsMeshSurfaceToPolygonAlgorithm::groupId() const +{ + return QStringLiteral( "mesh" ); +} + +QgsProcessingAlgorithm *QgsMeshSurfaceToPolygonAlgorithm::createInstance() const +{ + return new QgsMeshSurfaceToPolygonAlgorithm(); +} + +void QgsMeshSurfaceToPolygonAlgorithm::initAlgorithm( const QVariantMap &configuration ) +{ + Q_UNUSED( configuration ); + + addParameter( new QgsProcessingParameterMeshLayer( QStringLiteral( "INPUT" ), QObject::tr( "Input mesh layer" ) ) ); + + addParameter( new QgsProcessingParameterCrs( QStringLiteral( "CRS_OUTPUT" ), QObject::tr( "Output coordinate system" ), QVariant(), true ) ); + + addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output vector layer" ), Qgis::ProcessingSourceType::VectorPolygon ) ); +} + +bool QgsMeshSurfaceToPolygonAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback * ) +{ + QgsMeshLayer *meshLayer = parameterAsMeshLayer( parameters, QStringLiteral( "INPUT" ), context ); + + if ( !meshLayer || !meshLayer->isValid() ) + return false; + + if ( meshLayer->isEditable() ) + throw QgsProcessingException( QObject::tr( "Input mesh layer in edit mode is not supported" ) ); + + QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, QStringLiteral( "CRS_OUTPUT" ), context ); + if ( !outputCrs.isValid() ) + outputCrs = meshLayer->crs(); + mTransform = QgsCoordinateTransform( meshLayer->crs(), outputCrs, context.transformContext() ); + if ( !meshLayer->nativeMesh() ) + meshLayer->updateTriangularMesh( mTransform ); //necessary to load the native mesh + + mNativeMesh = *meshLayer->nativeMesh(); + + return true; +} + + +QVariantMap QgsMeshSurfaceToPolygonAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) +{ + if ( feedback ) + { + if ( feedback->isCanceled() ) + return QVariantMap(); + feedback->setProgress( 0 ); + feedback->pushInfo( QObject::tr( "Creating output vector layer" ) ); + } + + QgsCoordinateReferenceSystem outputCrs = parameterAsCrs( parameters, QStringLiteral( "CRS_OUTPUT" ), context ); + QString identifier; + std::unique_ptr sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, identifier, QgsFields(), Qgis::WkbType::MultiPolygon, outputCrs ) ); + if ( !sink ) + return QVariantMap(); + + if ( feedback ) + { + if ( feedback->isCanceled() ) + return QVariantMap(); + feedback->setProgress( 0 ); + } + + QgsGeometry lines; + QgsMeshFace face; + QMap, int> edges; // edge as key and count of edge usage as value + std::pair edge; + + if ( feedback ) + feedback->setProgressText( QObject::tr( "Parsing mesh faces to extract edges." ) ); + + for ( int i = 0; i < mNativeMesh.faceCount(); i++ ) + { + if ( feedback ) + { + if ( feedback->isCanceled() ) + return QVariantMap(); + } + + face = mNativeMesh.face( i ); + + for ( int j = 0; j < face.size(); j++ ) + { + int indexEnd; + if ( j == face.size() - 1 ) + indexEnd = 0; + else + indexEnd = j + 1; + int edgeFirstVertex = face.at( j ); + int edgeSecondVertex = face.at( indexEnd ); + + // make vertex sorted to avoid have 1,2 and 2,1 as different keys + if ( edgeSecondVertex < edgeFirstVertex ) + edge = std::make_pair( edgeSecondVertex, edgeFirstVertex ); + else + edge = std::make_pair( edgeFirstVertex, edgeSecondVertex ); + + // if edge exist in map increase its count otherwise set count to 1 + auto it = edges.find( edge ); + if ( it != edges.end() ) + { + int count = edges.take( edge ) + 1; + edges.insert( edge, count ); + } + else + { + edges.insert( edge, 1 ); + } + } + + if ( feedback ) + feedback->setProgress( 100.0 * static_cast( i ) / mNativeMesh.faceCount() ); + } + + if ( feedback ) + { + feedback->setProgress( 0 ); + feedback->setProgressText( QObject::tr( "Parsing mesh edges." ) ); + } + + std::unique_ptr multiLineString( new QgsMultiLineString() ); + + int i = 0; + for ( auto it = edges.begin(); it != edges.end(); it++ ) + { + if ( feedback ) + { + if ( feedback->isCanceled() ) + return QVariantMap(); + } + + // only consider edges with count 1 which are on the edge of mesh surface + if ( it.value() == 1 ) + { + std::unique_ptr line( new QgsLineString( mNativeMesh.vertex( it.key().first ), mNativeMesh.vertex( it.key().second ) ) ); + multiLineString->addGeometry( line.release() ); + } + + if ( feedback ) + feedback->setProgress( 100.0 * static_cast( i ) / edges.size() ); + + i++; + } + + if ( feedback ) + { + feedback->setProgressText( QObject::tr( "Creating final geometry." ) ); + if ( feedback->isCanceled() ) + return QVariantMap(); + } + + // merge lines + QgsGeometry mergedLines = QgsGeometry( multiLineString.release() ); + mergedLines = mergedLines.mergeLines(); + QgsAbstractGeometry *multiLinesAbstract = mergedLines.get(); + + // set of polygons to construct result + QVector polygons; + + // for every part create polygon and add to resulting multipolygon + for ( auto pit = multiLinesAbstract->const_parts_begin(); pit != multiLinesAbstract->const_parts_end(); ++pit ) + { + if ( feedback ) + { + if ( feedback->isCanceled() ) + return QVariantMap(); + } + + // individula polygon - can be either polygon or hole in polygon + QgsPolygon *polygon = new QgsPolygon(); + polygon->setExteriorRing( qgsgeometry_cast( *pit )->clone() ); + + // add first polygon, no need to check anything + if ( polygons.empty() ) + { + polygons.push_back( polygon ); + continue; + } + + // engine for spatial relations + std::unique_ptr engine( QgsGeometry::createGeometryEngine( polygon ) ); + + // need to check if polygon is not either contained (hole) or covering (main polygon) with another + // this solves meshes with holes + bool isHole = false; + + for ( int i = 0; i < polygons.count(); i++ ) + { + QgsPolygon *p = qgsgeometry_cast( polygons.at( i ) ); + + // polygon covers another, turn contained polygon into interior ring + if ( engine->contains( p ) ) + { + polygons.removeAt( i ); + polygon->addInteriorRing( p->exteriorRing()->clone() ); + break; + } + // polygon is within another, make it interior rind and do not add it + else if ( engine->within( p ) ) + { + p->addInteriorRing( polygon->exteriorRing()->clone() ); + isHole = true; + break; + } + } + + // if is not a hole polygon add it to the vector of polygons + if ( !isHole ) + polygons.append( polygon ); + else + // polygon was used as a hole, it is not needed anymore, delete it to avoid memory leak + delete polygon; + } + + // create resulting multipolygon + std::unique_ptr multiPolygon = std::make_unique(); + multiPolygon->addGeometries( polygons ); + + if ( feedback ) + { + if ( feedback->isCanceled() ) + return QVariantMap(); + } + + // create final geom and transform it + QgsGeometry resultGeom = QgsGeometry( multiPolygon.release() ); + + try + { + resultGeom.transform( mTransform ); + } + catch ( QgsCsException & ) + { + if ( feedback ) + feedback->reportError( QObject::tr( "Could not transform point to destination CRS" ) ); + } + + QgsFeature feat; + feat.setGeometry( resultGeom ); + + if ( !sink->addFeature( feat, QgsFeatureSink::FastInsert ) ) + throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) ); + + if ( feedback ) + { + feedback->pushInfo( QObject::tr( "Output vector layer created" ) ); + if ( feedback->isCanceled() ) + return QVariantMap(); + } + + QVariantMap ret; + ret[QStringLiteral( "OUTPUT" )] = identifier; + + return ret; +} + +///@endcond PRIVATE diff --git a/src/analysis/processing/qgsalgorithmmeshsurfacetopolygon.h b/src/analysis/processing/qgsalgorithmmeshsurfacetopolygon.h new file mode 100644 index 000000000000..482d3c2b8dff --- /dev/null +++ b/src/analysis/processing/qgsalgorithmmeshsurfacetopolygon.h @@ -0,0 +1,50 @@ +/*************************************************************************** + qgsalgorithmmeshsurfacetopolygon.h + --------------------------- + begin : September 2024 + copyright : (C) 2024 by Jan Caha + email : jan.caha at outlook dot com + ***************************************************************************/ + +/*************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#ifndef QGSALGORITHMMESHCONVEXHULL_H +#define QGSALGORITHMMESHCONVEXHULL_H + +#define SIP_NO_FILE + +#include "qgsprocessingalgorithm.h" +#include "qgsmeshdataprovider.h" +#include "qgstriangularmesh.h" + +///@cond PRIVATE + +class QgsMeshSurfaceToPolygonAlgorithm : public QgsProcessingAlgorithm +{ + public: + QgsMeshSurfaceToPolygonAlgorithm() = default; + QString shortHelpString() const override; + QString name() const override; + QString displayName() const override; + QString group() const override; + QString groupId() const override; + void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override; + bool prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + QVariantMap processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + QgsProcessingAlgorithm *createInstance() const override SIP_FACTORY; + + private: + QgsMesh mNativeMesh; + QgsCoordinateTransform mTransform; +}; + +///@endcond PRIVATE + +#endif // QGSALGORITHMMESHCONVEXHULL_H diff --git a/src/analysis/processing/qgsnativealgorithms.cpp b/src/analysis/processing/qgsnativealgorithms.cpp index 933a0ed22ed6..00c65aebfdb0 100644 --- a/src/analysis/processing/qgsnativealgorithms.cpp +++ b/src/analysis/processing/qgsnativealgorithms.cpp @@ -129,6 +129,7 @@ #include "qgsalgorithmmeancoordinates.h" #include "qgsalgorithmmergelines.h" #include "qgsalgorithmmergevector.h" +#include "qgsalgorithmmeshsurfacetopolygon.h" #include "qgsalgorithmminimumenclosingcircle.h" #include "qgsalgorithmmultidifference.h" #include "qgsalgorithmmultiintersection.h" @@ -430,6 +431,7 @@ void QgsNativeAlgorithms::loadAlgorithms() addAlgorithm( new QgsMeshContoursAlgorithm ); addAlgorithm( new QgsMeshExportCrossSection ); addAlgorithm( new QgsMeshExportTimeSeries ); + addAlgorithm( new QgsMeshSurfaceToPolygonAlgorithm() ); addAlgorithm( new QgsMinimumEnclosingCircleAlgorithm() ); addAlgorithm( new QgsMultiDifferenceAlgorithm() ); addAlgorithm( new QgsMultiIntersectionAlgorithm() );