sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-04-19 00:04:52 -04:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #7194 from wonder-sk/snap-geometries-alg

Browse Source 
[FEATURE] Snap geometries algorithm
This commit is contained in:
Martin Dobias 2018-09-14 12:36:41 +02:00 committed by GitHub
commit 9bcd21fb19
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
11 changed files with 547 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
python/analysis/analysis_auto.sip
View File

@ -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

52
python/analysis/auto_generated/vector/qgsgeometrysnappersinglesource.sip.in Normal file
View File

@ -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 *
************************************************************************/

11
python/plugins/processing/algs/qgis/SnapGeometries.py
View File

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

BIN
python/plugins/processing/tests/testdata/custom/circular_strings.gpkg vendored
View File

Binary file not shown.

9
python/plugins/processing/tests/testdata/custom/snap_geometries.geojson vendored Normal file
View File

@ -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 ] ] ] } }
]
}

24
python/plugins/processing/tests/testdata/expected/snap_geometries.gml vendored Normal file
View File

@ -0,0 +1,24 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ snap_geometries.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>20</gml:X><gml:Y>5</gml:Y></gml:coord>
<gml:coord><gml:X>31</gml:X><gml:Y>20.5</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>
<gml:featureMember>
<ogr:snap_geometries fid="snap_geometries.0">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:3857"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>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</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
</ogr:snap_geometries>
</gml:featureMember>
<gml:featureMember>
<ogr:snap_geometries fid="snap_geometries.1">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:3857"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>31,10 30.0,10.123 21.0,10.123 20.123,10.123 20,5 31,10</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
</ogr:snap_geometries>
</gml:featureMember>
</ogr:FeatureCollection>

23
python/plugins/processing/tests/testdata/expected/snap_geometries.xsd vendored Normal file
View File

@ -0,0 +1,23 @@
<?xml version="1.0" encoding="UTF-8"?>
<xs:schema targetNamespace="http://ogr.maptools.org/" xmlns:ogr="http://ogr.maptools.org/" xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:gml="http://www.opengis.net/gml" elementFormDefault="qualified" version="1.0">
<xs:import namespace="http://www.opengis.net/gml" schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:_FeatureCollection"/>
<xs:complexType name="FeatureCollectionType">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureCollectionType">
<xs:attribute name="lockId" type="xs:string" use="optional"/>
<xs:attribute name="scope" type="xs:string" use="optional"/>
</xs:extension>
</xs:complexContent>
</xs:complexType>
<xs:element name="snap_geometries" type="ogr:snap_geometries_Type" substitutionGroup="gml:_Feature"/>
<xs:complexType name="snap_geometries_Type">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureType">
<xs:sequence>
<xs:element name="geometryProperty" type="gml:PolygonPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
</xs:sequence>
</xs:extension>
</xs:complexContent>
</xs:complexType>
</xs:schema>

16
python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml vendored
View File

@ -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:

2
src/analysis/CMakeLists.txt
View File

@ -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

356
src/analysis/vector/qgsgeometrysnappersinglesource.cpp Normal file
View File

@ -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<AnchorPoint> &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<QgsFeatureId> 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<AnchorPoint> &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<QgsFeatureId> 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<AnchorPoint> &pnts )
{
// Find point ( should always find one point )
QList<QgsFeatureId> 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<AnchorPoint> &pnts, double thresh )
{
QVector<QgsPoint> newPoints;
QVector<int> 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<QgsFeatureId> 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<QgsFeatureId> fids = index.intersects( rect );
QVector<AnchorAlongSegment> 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<AnchorPoint> &pnts, double thresh )
{
bool changed = false;
if ( QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( g ) )
{
changed |= snapLineString( linestring, index, pnts, thresh );
}
else if ( QgsPolygon *polygon = qgsgeometry_cast<QgsPolygon *>( g ) )
{
if ( QgsLineString *exteriorRing = qgsgeometry_cast<QgsLineString *>( polygon->exteriorRing() ) )
changed |= snapLineString( exteriorRing, index, pnts, thresh );
for ( int i = 0; i < polygon->numInteriorRings(); ++i )
{
if ( QgsLineString *interiorRing = qgsgeometry_cast<QgsLineString *>( polygon->interiorRing( i ) ) )
changed |= snapLineString( interiorRing, index, pnts, thresh );
}
}
else if ( QgsGeometryCollection *collection = qgsgeometry_cast<QgsGeometryCollection *>( g ) )
{
for ( int i = 0; i < collection->numGeometries(); ++i )
changed |= snapGeometry( collection->geometryN( i ), index, pnts, thresh );
}
else if ( QgsPoint *pt = qgsgeometry_cast<QgsPoint *>( 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<AnchorPoint> 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;
}

54
src/analysis/vector/qgsgeometrysnappersinglesource.h Normal file
View File

@ -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