diff --git a/python/analysis/analysis_auto.sip b/python/analysis/analysis_auto.sip index 5d92e3cbaf8..00ec4069e26 100644 --- a/python/analysis/analysis_auto.sip +++ b/python/analysis/analysis_auto.sip @@ -13,6 +13,7 @@ %Include auto_generated/raster/qgsrastercalcnode.sip %Include auto_generated/raster/qgstotalcurvaturefilter.sip %Include auto_generated/vector/qgsgeometrysnapper.sip +%Include auto_generated/vector/qgsgeometrysnappersinglesource.sip %Include auto_generated/vector/qgszonalstatistics.sip %Include auto_generated/interpolation/qgsinterpolator.sip %Include auto_generated/interpolation/qgsgridfilewriter.sip diff --git a/python/analysis/auto_generated/vector/qgsgeometrysnappersinglesource.sip.in b/python/analysis/auto_generated/vector/qgsgeometrysnappersinglesource.sip.in new file mode 100644 index 00000000000..c3391561b94 --- /dev/null +++ b/python/analysis/auto_generated/vector/qgsgeometrysnappersinglesource.sip.in @@ -0,0 +1,52 @@ +/************************************************************************ + * This file has been generated automatically from * + * * + * src/analysis/vector/qgsgeometrysnappersinglesource.h * + * * + * Do not edit manually ! Edit header and run scripts/sipify.pl again * + ************************************************************************/ + + + + +class QgsGeometrySnapperSingleSource +{ +%Docstring + +Makes sure that any two vertices of the vector layer are at least at distance given by the threshold value. +The algorithm moves nearby vertices to one location and adds vertices to segments that are passing around other +vertices within the threshold. It does not remove any vertices. Also, it does not modify geometries unless +needed (it does not snap coordinates to a grid). + +This algorithm comes handy when doing vector overlay operations such as intersection, union or difference +to prevent possible topological errors caused by numerical errors if coordinates are very close to each other. + +After running the algorithm some previously valid geometries may become invalid and therefore it may be useful +to run Fix geometries algorithm afterwards. + +.. note:: + + Originally ported from GRASS implementation of Vect_snap_lines_list() + +.. versionadded:: 3.4 +%End + +%TypeHeaderCode +#include "qgsgeometrysnappersinglesource.h" +%End + public: + + static int run( const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback ); +%Docstring +Run the algorithm on given source and output results to the sink, using threshold value in the source's map units. +Returns number of modified geometries. +%End +}; + +/************************************************************************ + * This file has been generated automatically from * + * * + * src/analysis/vector/qgsgeometrysnappersinglesource.h * + * * + * Do not edit manually ! Edit header and run scripts/sipify.pl again * + ************************************************************************/ diff --git a/python/plugins/processing/algs/qgis/SnapGeometries.py b/python/plugins/processing/algs/qgis/SnapGeometries.py index 9bb85c4c24d..80c270697ca 100644 --- a/python/plugins/processing/algs/qgis/SnapGeometries.py +++ b/python/plugins/processing/algs/qgis/SnapGeometries.py @@ -26,6 +26,7 @@ __copyright__ = '(C) 2016, Nyall Dawson' __revision__ = '$Format:%H$' from qgis.analysis import (QgsGeometrySnapper, + QgsGeometrySnapperSingleSource, QgsInternalGeometrySnapper) from qgis.core import (QgsFeatureSink, QgsProcessing, @@ -71,7 +72,8 @@ class SnapGeometriesToLayer(QgisAlgorithm): self.tr('Prefer closest point, don\'t insert new vertices'), self.tr('Move end points only, prefer aligning nodes'), self.tr('Move end points only, prefer closest point'), - self.tr('Snap end points to end points only')] + self.tr('Snap end points to end points only'), + self.tr('Snap to anchor nodes (single layer only)')] self.addParameter(QgsProcessingParameterEnum( self.BEHAVIOR, self.tr('Behavior'), @@ -106,6 +108,9 @@ class SnapGeometriesToLayer(QgisAlgorithm): total = 100.0 / source.featureCount() if source.featureCount() else 0 if parameters[self.INPUT] != parameters[self.REFERENCE_LAYER]: + if mode == 7: + raise QgsProcessingException(self.tr('This mode applies when the input and reference layer are the same.')) + snapper = QgsGeometrySnapper(reference_source) processed = 0 for f in features: @@ -119,6 +124,10 @@ class SnapGeometriesToLayer(QgisAlgorithm): sink.addFeature(f) processed += 1 feedback.setProgress(processed * total) + elif mode == 7: + # input layer == ref layer + modified_count = QgsGeometrySnapperSingleSource.run(source, sink, tolerance, feedback) + feedback.pushInfo(self.tr('Snapped {} geometries.').format(modified_count)) else: # snapping internally snapper = QgsInternalGeometrySnapper(tolerance, mode) diff --git a/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg b/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg index fd342b87dd8..ff47737efa6 100644 Binary files a/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg and b/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg differ diff --git a/python/plugins/processing/tests/testdata/custom/snap_geometries.geojson b/python/plugins/processing/tests/testdata/custom/snap_geometries.geojson new file mode 100644 index 00000000000..c661eda464d --- /dev/null +++ b/python/plugins/processing/tests/testdata/custom/snap_geometries.geojson @@ -0,0 +1,9 @@ +{ +"type": "FeatureCollection", +"name": "snap_geometries", +"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3857" } }, +"features": [ +{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 20.123, 10.123 ], [ 20.123, 20.456 ], [ 20.2, 20.456 ], [ 20.2, 20.5 ], [ 30.0, 20.5 ], [ 30.0, 10.123 ], [ 29.8, 10.123 ], [ 21.0, 10.123 ], [ 20.123, 10.123 ] ] ] } }, +{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 31.0, 10.0 ], [ 20.0, 10.0 ], [ 20.0, 5.0 ], [ 31.0, 10.0 ] ] ] } } +] +} diff --git a/python/plugins/processing/tests/testdata/expected/snap_geometries.gml b/python/plugins/processing/tests/testdata/expected/snap_geometries.gml new file mode 100644 index 00000000000..610846d3811 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/snap_geometries.gml @@ -0,0 +1,24 @@ + + + + + 205 + 3120.5 + + + + + + 20.123,10.123 20.123,20.456 20.123,20.456 20.123,20.456 30.0,20.5 30.0,10.123 30.0,10.123 21.0,10.123 20.123,10.123 + + + + + 31,10 30.0,10.123 21.0,10.123 20.123,10.123 20,5 31,10 + + + diff --git a/python/plugins/processing/tests/testdata/expected/snap_geometries.xsd b/python/plugins/processing/tests/testdata/expected/snap_geometries.xsd new file mode 100644 index 00000000000..c27240ec294 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/snap_geometries.xsd @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml index 749d524d777..b8f1ea60a6e 100755 --- a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml +++ b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml @@ -2045,6 +2045,22 @@ tests: name: expected/snap_internal.gml type: vector + - algorithm: qgis:snapgeometries + name: Test Snap Geometries (to each other) + params: + BEHAVIOR: '7' + INPUT: + name: custom/snap_geometries.geojson + type: vector + REFERENCE_LAYER: + name: custom/snap_geometries.geojson + type: vector + TOLERANCE: 0.5 + results: + OUTPUT: + name: expected/snap_geometries.gml + type: vector + - algorithm: qgis:poleofinaccessibility name: Pole of inaccessibility (polygons) params: diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index f8b89ad2bfd..e6cb0ebd3c6 100755 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -123,6 +123,7 @@ SET(QGIS_ANALYSIS_SRCS raster/qgsrastermatrix.cpp vector/mersenne-twister.cpp vector/qgsgeometrysnapper.cpp + vector/qgsgeometrysnappersinglesource.cpp vector/qgszonalstatistics.cpp network/qgsgraph.cpp @@ -231,6 +232,7 @@ SET(QGIS_ANALYSIS_HDRS vector/mersenne-twister.h vector/qgsgeometrysnapper.h + vector/qgsgeometrysnappersinglesource.h vector/qgszonalstatistics.h vector/geometry_checker/qgsgeometrycheckerutils.h vector/geometry_checker/qgsfeaturepool.h diff --git a/src/analysis/vector/qgsgeometrysnappersinglesource.cpp b/src/analysis/vector/qgsgeometrysnappersinglesource.cpp new file mode 100644 index 00000000000..a08eaddb24c --- /dev/null +++ b/src/analysis/vector/qgsgeometrysnappersinglesource.cpp @@ -0,0 +1,356 @@ +/*************************************************************************** + qgsgeometrysnappersinglesource.cpp + --------------------- + Date : May 2018 + Copyright : (C) 2018 by Martin Dobias + Email : wonder dot sk at gmail dot com + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include "qgsgeometrysnappersinglesource.h" + +#include "qgsfeatureiterator.h" +#include "qgsfeaturesink.h" +#include "qgsfeaturesource.h" +#include "qgsfeedback.h" +#include "qgsgeometrycollection.h" +#include "qgsgeometryutils.h" +#include "qgslinestring.h" +#include "qgspolygon.h" +#include "qgsspatialindex.h" + +//! record about vertex coordinates and index of anchor to which it is snapped +struct AnchorPoint +{ + //! coordinates of the point + double x, y; + + /** + * Anchor information: + * 0+ - index of anchor to which this point should be snapped + * -1 - initial value (undefined) + * -2 - this point is an anchor, i.e. do not snap this point (snap others to this point) + */ + int anchor; +}; + + +//! record about anchor being along a segment +struct AnchorAlongSegment +{ + int anchor; //!< Index of the anchor point + double along; //!< Distance of the anchor point along the segment +}; + + +static void buildSnapIndex( QgsFeatureIterator &fi, QgsSpatialIndex &index, QVector &pnts, QgsFeedback *feedback, int &count, int totalCount ) +{ + QgsFeature f; + int pntId = 0; + + while ( fi.nextFeature( f ) ) + { + if ( feedback->isCanceled() ) + break; + + QgsGeometry g = f.geometry(); + + for ( auto it = g.vertices_begin(); it != g.vertices_end(); ++it ) + { + QgsPoint pt = *it; + QgsRectangle rect( pt.x(), pt.y(), pt.x(), pt.y() ); + + QList ids = index.intersects( rect ); + if ( ids.isEmpty() ) + { + // add to tree and to structure + index.insertFeature( pntId, pt.boundingBox() ); + + AnchorPoint xp; + xp.x = pt.x(); + xp.y = pt.y(); + xp.anchor = -1; + pnts.append( xp ); + pntId++; + } + } + + ++count; + feedback->setProgress( 100. * count / totalCount ); + } +} + + +static void assignAnchors( QgsSpatialIndex &index, QVector &pnts, double thresh ) +{ + double thresh2 = thresh * thresh; + int nanchors = 0, ntosnap = 0; + for ( int point = 0; point < pnts.count(); ++point ) + { + if ( pnts[point].anchor >= 0 ) + continue; + + pnts[point].anchor = -2; // make it anchor + nanchors++; + + // Find points in threshold + double x = pnts[point].x, y = pnts[point].y; + QgsRectangle rect( x - thresh, y - thresh, x + thresh, y + thresh ); + + const QList ids = index.intersects( rect ); + for ( QgsFeatureId pointb : ids ) + { + if ( pointb == point ) + continue; + + double dx = pnts[pointb].x - pnts[point].x; + double dy = pnts[pointb].y - pnts[point].y; + double dist2 = dx * dx + dy * dy; + if ( dist2 > thresh2 ) + continue; // outside threshold + + if ( pnts[pointb].anchor == -1 ) + { + // doesn't have an anchor yet + pnts[pointb].anchor = point; + ntosnap++; + } + else if ( pnts[pointb].anchor >= 0 ) + { + // check distance to previously assigned anchor + double dx2 = pnts[pnts[pointb].anchor].x - pnts[pointb].x; + double dy2 = pnts[pnts[pointb].anchor].y - pnts[pointb].y; + double dist2_a = dx2 * dx2 + dy2 * dy2; + if ( dist2 < dist2_a ) + pnts[pointb].anchor = point; // replace old anchor + } + } + } +} + + +static bool snapPoint( QgsPoint *pt, QgsSpatialIndex &index, QVector &pnts ) +{ + // Find point ( should always find one point ) + QList fids = index.intersects( QgsRectangle( pt->x(), pt->y(), pt->x(), pt->y() ) ); + Q_ASSERT( fids.count() == 1 ); + + int spoint = fids[0]; + int anchor = pnts[spoint].anchor; + + if ( anchor >= 0 ) + { + // to be snapped + pt->setX( pnts[anchor].x ); + pt->setY( pnts[anchor].y ); + return true; + } + + return false; +} + + +static bool snapLineString( QgsLineString *linestring, QgsSpatialIndex &index, QVector &pnts, double thresh ) +{ + QVector newPoints; + QVector anchors; // indexes of anchors for vertices + double thresh2 = thresh * thresh; + double minDistX, minDistY; // coordinates of the closest point on the segment line + bool changed = false; + + // snap vertices + for ( int v = 0; v < linestring->numPoints(); v++ ) + { + double x = linestring->xAt( v ); + double y = linestring->yAt( v ); + QgsRectangle rect( x, y, x, y ); + + // Find point ( should always find one point ) + QList fids = index.intersects( rect ); + Q_ASSERT( fids.count() == 1 ); + + int spoint = fids.first(); + int anchor = pnts[spoint].anchor; + if ( anchor >= 0 ) + { + // to be snapped + linestring->setXAt( v, pnts[anchor].x ); + linestring->setYAt( v, pnts[anchor].y ); + anchors.append( anchor ); // point on new location + changed = true; + } + else + { + anchors.append( spoint ); // old point + } + } + + // Snap all segments to anchors in threshold + for ( int v = 0; v < linestring->numPoints() - 1; v++ ) + { + double x1 = linestring->xAt( v ); + double x2 = linestring->xAt( v + 1 ); + double y1 = linestring->yAt( v ); + double y2 = linestring->yAt( v + 1 ); + + newPoints << linestring->pointN( v ); + + // Box + double xmin = x1, xmax = x2, ymin = y1, ymax = y2; + if ( xmin > xmax ) + std::swap( xmin, xmax ); + if ( ymin > ymax ) + std::swap( ymin, ymax ); + + QgsRectangle rect( xmin - thresh, ymin - thresh, xmax + thresh, ymax + thresh ); + + // Find points + const QList fids = index.intersects( rect ); + + QVector newVerticesAlongSegment; + + // Snap to anchor in threshold different from end points + for ( QgsFeatureId fid : fids ) + { + int spoint = fid; + + if ( spoint == anchors[v] || spoint == anchors[v + 1] ) + continue; // end point + if ( pnts[spoint].anchor >= 0 ) + continue; // point is not anchor + + // Check the distance + double dist2 = QgsGeometryUtils::sqrDistToLine( pnts[spoint].x, pnts[spoint].y, x1, y1, x2, y2, minDistX, minDistY, 0 ); + // skip points that are behind segment's endpoints or extremely close to them + double dx1 = minDistX - x1, dx2 = minDistX - x2; + double dy1 = minDistY - y1, dy2 = minDistY - y2; + bool isOnSegment = !qgsDoubleNear( dx1 * dx1 + dy1 * dy1, 0 ) && !qgsDoubleNear( dx2 * dx2 + dy2 * dy2, 0 ); + if ( isOnSegment && dist2 <= thresh2 ) + { + // an anchor is in the threshold + AnchorAlongSegment item; + item.anchor = spoint; + item.along = QgsPointXY( x1, y1 ).distance( minDistX, minDistY ); + newVerticesAlongSegment << item; + } + } + + if ( !newVerticesAlongSegment.isEmpty() ) + { + // sort by distance along the segment + std::sort( newVerticesAlongSegment.begin(), newVerticesAlongSegment.end(), []( const AnchorAlongSegment & p1, const AnchorAlongSegment & p2 ) + { + return ( p1.along < p2.along ? -1 : ( p1.along > p2.along ) ); + } ); + + // insert new vertices + for ( int i = 0; i < newVerticesAlongSegment.count(); i++ ) + { + int anchor = newVerticesAlongSegment[i].anchor; + newPoints << QgsPoint( pnts[anchor].x, pnts[anchor].y, 0 ); + } + changed = true; + } + } + + // append end point + newPoints << linestring->pointN( linestring->numPoints() - 1 ); + + // replace linestring's points + if ( changed ) + linestring->setPoints( newPoints ); + + return changed; +} + + +static bool snapGeometry( QgsAbstractGeometry *g, QgsSpatialIndex &index, QVector &pnts, double thresh ) +{ + bool changed = false; + if ( QgsLineString *linestring = qgsgeometry_cast( g ) ) + { + changed |= snapLineString( linestring, index, pnts, thresh ); + } + else if ( QgsPolygon *polygon = qgsgeometry_cast( g ) ) + { + if ( QgsLineString *exteriorRing = qgsgeometry_cast( polygon->exteriorRing() ) ) + changed |= snapLineString( exteriorRing, index, pnts, thresh ); + for ( int i = 0; i < polygon->numInteriorRings(); ++i ) + { + if ( QgsLineString *interiorRing = qgsgeometry_cast( polygon->interiorRing( i ) ) ) + changed |= snapLineString( interiorRing, index, pnts, thresh ); + } + } + else if ( QgsGeometryCollection *collection = qgsgeometry_cast( g ) ) + { + for ( int i = 0; i < collection->numGeometries(); ++i ) + changed |= snapGeometry( collection->geometryN( i ), index, pnts, thresh ); + } + else if ( QgsPoint *pt = qgsgeometry_cast( g ) ) + { + changed |= snapPoint( pt, index, pnts ); + } + + return changed; +} + + +int QgsGeometrySnapperSingleSource::run( const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback ) +{ + // the logic here comes from GRASS implementation of Vect_snap_lines_list() + + int count = 0; + int totalCount = source.featureCount() * 2; + + // step 1: record all point locations in a spatial index + extra data structure to keep + // reference to which other point they have been snapped to (in the next phase). + + QgsSpatialIndex index; + QVector pnts; + QgsFeatureRequest request; + request.setSubsetOfAttributes( QgsAttributeList() ); + QgsFeatureIterator fi = source.getFeatures( request ); + buildSnapIndex( fi, index, pnts, feedback, count, totalCount ); + + if ( feedback->isCanceled() ) + return 0; + + // step 2: go through all registered points and if not yet marked mark it as anchor and + // assign this anchor to all not yet marked points in threshold + + assignAnchors( index, pnts, thresh ); + + // step 3: alignment of vertices and segments to the anchors + // Go through all lines and: + // 1) for all vertices: if not anchor snap it to its anchor + // 2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) + + int modified = 0; + QgsFeature f; + fi = source.getFeatures(); + while ( fi.nextFeature( f ) ) + { + if ( feedback->isCanceled() ) + break; + + QgsGeometry geom = f.geometry(); + if ( snapGeometry( geom.get(), index, pnts, thresh ) ) + { + f.setGeometry( geom ); + ++modified; + } + + sink.addFeature( f, QgsFeatureSink::FastInsert ); + + ++count; + feedback->setProgress( 100. * count / totalCount ); + } + + return modified; +} diff --git a/src/analysis/vector/qgsgeometrysnappersinglesource.h b/src/analysis/vector/qgsgeometrysnappersinglesource.h new file mode 100644 index 00000000000..f413d02cfba --- /dev/null +++ b/src/analysis/vector/qgsgeometrysnappersinglesource.h @@ -0,0 +1,54 @@ +/*************************************************************************** + qgsgeometrysnappersinglesource.h + --------------------- + Date : May 2018 + Copyright : (C) 2018 by Martin Dobias + Email : wonder dot sk at gmail dot com + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#ifndef QGSGEOMETRYSNAPPERSINGLESOURCE_H +#define QGSGEOMETRYSNAPPERSINGLESOURCE_H + +#include "qgis_analysis.h" + +class QgsFeatureSink; +class QgsFeatureSource; +class QgsFeedback; + +/** + * \ingroup analysis + * + * Makes sure that any two vertices of the vector layer are at least at distance given by the threshold value. + * The algorithm moves nearby vertices to one location and adds vertices to segments that are passing around other + * vertices within the threshold. It does not remove any vertices. Also, it does not modify geometries unless + * needed (it does not snap coordinates to a grid). + * + * This algorithm comes handy when doing vector overlay operations such as intersection, union or difference + * to prevent possible topological errors caused by numerical errors if coordinates are very close to each other. + * + * After running the algorithm some previously valid geometries may become invalid and therefore it may be useful + * to run Fix geometries algorithm afterwards. + * + * \note Originally ported from GRASS implementation of Vect_snap_lines_list() + * + * \since QGIS 3.4 + */ +class ANALYSIS_EXPORT QgsGeometrySnapperSingleSource +{ + public: + + /** + * Run the algorithm on given source and output results to the sink, using threshold value in the source's map units. + * Returns number of modified geometries. + */ + static int run( const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback ); +}; + +#endif // QGSGEOMETRYSNAPPERSINGLESOURCE_H