New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Improve split with lines #3539
Improve split with lines #3539
Changes from 2 commits
63e3fd3
6cf47be
c282e26
5093ec6
164a85a
45711d3
6b2b4c5
6ea0049
be6672d
51a5657
73b283c
b0644ea
88587fd
1cd1158
10648df
dc800db
cc0b2e6
a3ae0b2
1ef7ed5
72118f9
132e76a
5625d6e
156fce9
be2223f
35d106b
f3d7e39
a4cc572
28d7cea
054d430
45617fb
4981bfc
4e96912
fb124ba
b069e95
c180cf3
08231b8
159fda6
5991ecc
a6bd9f0
31a6189
23de13c
52e29b9
d08c02d
397df65
02ed0aa
32094e9
d239a97
a8feec8
256efcc
e1ee477
652addb
85d1fd7
b6d5f35
dff239c
426c5be
2287230
75bd622
fcc96de
af16bbf
0589566
d81533e
fbc12a8
3ef7b3b
f24cda4
dbf6107
8b1adc5
d729951
9872b48
514d443
1d3f1f0
d5c307e
d6f09c0
b440939
880647e
1f81a7c
798bc09
6e9631e
7ae6269
7299e6b
968e02d
b13b4f9
9ad365e
9afd6b0
9950b08
645514a
8898c94
dd68d81
85c8c97
d3f5314
4a5faa0
5d78d60
29d33b4
d237e27
8908eea
c87839f
4a4ffa4
38a4aac
afd5d1e
5e3d8fe
08eaeaa
b0801f3
5346023
b5c1d0f
8f8624a
9ddf78e
a373f95
6bfd0fd
adcb772
6c4eeda
e9c4090
1d538a6
87d2ac3
6727ea7
5e487cf
5b4a88f
e7d8b2c
08d350c
a6eb7b6
f18d1c3
f5f2850
b32a719
f2e3d53
582a56d
ff691a6
3242321
4682eaf
747097d
95271c8
b5864cd
091f488
631bd48
249c8dc
60bbd09
af016cf
33ee514
6638569
959f97f
ec49341
9ee7873
5c3aea3
f78f2a0
fc65334
1da400a
089b663
964ecfd
e624518
fb8b931
24ffa15
a8a05ba
6b6896b
1ffdadd
aea5770
bf43d3d
10db745
ebcfa72
5f8fd2b
263ba81
5e1a69f
13a0e48
034e4f3
72b4e72
5f3ba72
8fcf834
14cfa26
f2ed214
fc18fd9
f70a3b9
9679b6a
cac8de5
b4533cd
a958c8a
08f8ca7
19368cf
cc985c2
2260780
4c2725b
e9fa3e0
74dfd1f
9bb3235
a510516
5df9cbc
64b6a44
c123d3b
2270603
eca83e3
6d0e8e6
eda412d
ee71077
e27e533
3c3e17a
14c930a
5cbf9d5
d0f8863
137b198
5992f74
c1b6edc
f77ab4d
00eb261
f9be179
a61e8bb
9dffe64
d657c77
b421a53
fb5cdd8
82082b4
22dc096
8c8be00
3dcf03c
271e67e
188033a
e5f62e4
377cba0
86ab302
20dc7fb
33caf66
0c58555
d716640
efd2992
bc130be
0a0b3a7
e226529
8dbe6ea
2894236
e3e70f7
b7fa540
dc82ad4
ee80be4
311f482
d712b49
3c51a93
489e00d
2652aa5
eeafe67
d559d7f
f71430e
a972c6d
ae02f51
621e366
2f445a3
a4f4595
7eabb35
81d7d18
1838a71
3dd57e8
c20d2cb
44a5f07
09e9e51
783553f
7f5a011
c710f3c
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
This file was deleted.
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,173 @@ | ||
# -*- coding: utf-8 -*- | ||
|
||
""" | ||
*************************************************************************** | ||
SplitWithLines.py | ||
--------------------- | ||
Date : November 2014 | ||
Revised : September 2016 | ||
Copyright : (C) 2014 by Bernhard Ströbl | ||
Email : bernhard dot stroebl at jena dot de | ||
*************************************************************************** | ||
* * | ||
* 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. * | ||
* * | ||
*************************************************************************** | ||
""" | ||
|
||
__author__ = 'Bernhard Ströbl' | ||
__date__ = 'November 2014' | ||
__copyright__ = '(C) 2014, Bernhard Ströbl' | ||
|
||
# This will get replaced with a git SHA1 when you do a git archive | ||
|
||
__revision__ = '$Format:%H$' | ||
|
||
from qgis.core import Qgis, QgsFeatureRequest, QgsFeature, QgsGeometry, QgsWkbTypes | ||
from processing.core.GeoAlgorithm import GeoAlgorithm | ||
from processing.core.parameters import ParameterVector | ||
from processing.core.outputs import OutputVector | ||
from processing.core.ProcessingLog import ProcessingLog | ||
from processing.tools import dataobjects | ||
from processing.tools import vector | ||
|
||
|
||
class SplitWithLines(GeoAlgorithm): | ||
|
||
INPUT_A = 'INPUT_A' | ||
INPUT_B = 'INPUT_B' | ||
|
||
OUTPUT = 'OUTPUT' | ||
|
||
def defineCharacteristics(self): | ||
self.name, self.i18n_name = self.trAlgorithm('Split with lines') | ||
self.group, self.i18n_group = self.trAlgorithm('Vector overlay tools') | ||
self.addParameter(ParameterVector(self.INPUT_A, | ||
self.tr('Input layer, single geometries only'), [dataobjects.TYPE_VECTOR_POLYGON, | ||
dataobjects.TYPE_VECTOR_LINE])) | ||
self.addParameter(ParameterVector(self.INPUT_B, | ||
self.tr('Split layer'), [dataobjects.TYPE_VECTOR_LINE])) | ||
|
||
self.addOutput(OutputVector(self.OUTPUT, self.tr('Splitted'))) | ||
|
||
def processAlgorithm(self, progress): | ||
layerA = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT_A)) | ||
splitLayer = dataobjects.getObjectFromUri(self.getParameterValue(self.INPUT_B)) | ||
|
||
sameLayer = self.getParameterValue(self.INPUT_A) == self.getParameterValue(self.INPUT_B) | ||
fieldList = layerA.fields() | ||
|
||
writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fieldList, | ||
layerA.wkbType(), layerA.crs()) | ||
|
||
spatialIndex = vector.spatialindex(splitLayer) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Unfortunately this + line 71 means that all features are fetched from the provider twice - this could be a performance killer on large layers. I'd suggest you could create the spatialIndex here, and then add features to the index from within the loop at line 72. The index creation will be slower then the current approach, but the boost from only fetching features once should more than offset this. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I was under the impression that it would be wise to use the tools processing provides for certain tasks and thus did not think about implementing it myself. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thanks! The inbuilt processing methods are fine for simple use cases (eg looping once through features and applying an operation to each), but to really optimise complex algs you need to use the qgis api itself. |
||
splitGeoms = {} | ||
|
||
for aSplitFeature in vector.features(splitLayer): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. You could speed this up a bit by passing a feature request to vector.features which uses |
||
splitGeoms[aSplitFeature.id()] = aSplitFeature.geometry() | ||
# honor the case that user has selection on split layer and has setting "use selection" | ||
|
||
outFeat = QgsFeature() | ||
features = vector.features(layerA) | ||
total = 100.0 / float(len(features)) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. you'll need to handle the case that len(features)==0 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. You're right, although this would only happen with an empty layer. This was basically a copy-paste operation from some other algo, so I expect to see this line in others, too. But will fix There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yeah, that's how I know to look out for it... I hit it in some other algorithm with the same block There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. :-) |
||
allowedWkbTypes = [2, #WkbLineString | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This looks fragile, and doesn't handle Z/M geometries - I'd suggest instead checking the QgsGeometry.type() for QgsWkbTypes.LineGeometry and PolygonGeometry instead, and then using QgsGeometry.isMultipart() to avoid the multipart types. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Wouldn't it be either Line (type = 1) or Polygon (type = 2) anyways because INPUT_A is set to be either There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. basically I would prefer to disable multi-layers to be used as input in the first place, but there is no method to do so AFAIK. |
||
3, # WkbPolygon | ||
-2147483646, #WkbLineString25D | ||
-2147483645] # WkbPolygon25D | ||
# MultiGeometries are not allowed because the result of a splitted ring is not clearly defined: | ||
# 1) add both parts as new features | ||
# 2) store one part as a new feature and the other one as ring of the multi geometry | ||
# 2a) which part is which, seems arbitrary | ||
|
||
multiGeoms = 0 # how many multi geometries were encountered | ||
|
||
for current, inFeatA in enumerate(features): | ||
inGeom = inFeatA.geometry() | ||
|
||
if allowedWkbTypes.count(inGeom.wkbType()) == 0: | ||
multiGeoms += 1 | ||
else: | ||
attrsA = inFeatA.attributes() | ||
outFeat.setAttributes(attrsA) | ||
inGeoms = [inGeom] | ||
lines = spatialIndex.intersects(inGeom.boundingBox()) | ||
|
||
if len(lines) > 0: # has intersection of bounding boxes | ||
splittingLines = [] | ||
|
||
for i in lines: | ||
try: | ||
splitGeom = splitGeoms[i] | ||
except: | ||
continue | ||
|
||
# check if trying to self-intersect | ||
if sameLayer: | ||
if inFeatA.id() == i: | ||
continue | ||
|
||
if inGeom.intersects(splitGeom): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If you want to make this algorithm magnitudes faster you could use prepared geometries here - have a look at https://github.com/qgis/QGIS/blob/master/python/plugins/processing/algs/qgis/Union.py#L110 and https://github.com/qgis/QGIS/blob/master/python/plugins/processing/algs/qgis/Union.py#L119 for how this is done There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @nyalldawson I tried to implement this but get a runtime error There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. You'll need to call There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. That did the trick! |
||
splittingLines.append(splitGeom) | ||
|
||
if len(splittingLines) > 0: | ||
for splitGeom in splittingLines: | ||
splitterPList = vector.extractPoints(splitGeom) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If you want to optimise further - you could set splitterPList to None here, and then only extract the points when splitterPList is required (ie, if it's still None). Extracting points is slow and they may not be used if none of the intersects checks on line 126 pass. |
||
outGeoms = [] | ||
|
||
while len(inGeoms) > 0: | ||
inGeom = inGeoms.pop() | ||
inPoints = vector.extractPoints(inGeom) | ||
|
||
if inGeom.intersects(splitGeom): | ||
try: | ||
result, newGeometries, topoTestPoints = inGeom.splitGeometry(splitterPList, False) | ||
except: | ||
ProcessingLog.addToLog(ProcessingLog.LOG_WARNING, | ||
self.tr('Geometry exception while splitting')) | ||
result = 1 | ||
|
||
# splitGeometry: If there are several intersections | ||
# between geometry and splitLine, only the first one is considered. | ||
if result == 0: # split occurred | ||
|
||
if inPoints == vector.extractPoints(inGeom): | ||
# bug in splitGeometry: sometimes it returns 0 but | ||
# the geometry is unchanged | ||
outGeoms.append(inGeom) | ||
else: | ||
inGeoms.append(inGeom) | ||
|
||
for aNewGeom in newGeometries: | ||
inGeoms.append(aNewGeom) | ||
else: | ||
outGeoms.append(inGeom) | ||
else: | ||
outGeoms.append(inGeom) | ||
|
||
inGeoms = outGeoms | ||
|
||
for aGeom in inGeoms: | ||
passed = True | ||
|
||
if aGeom.wkbType == 2 or aGeom.wkbType == -2147483646: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Would be better to use the proper enum values from QgsWkbTypes here - and better to check the flatType so that Z/M/ZM geometry types will be properly handled. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @nyalldawson You mean to only check for WKBLineString? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. At the moment this is only passing LineString or LineString25D. I'd suggest changing it to
Then it'll support LineStringZ, LineStringM, ...ZM, curved types, etc, and be future proof in case more are added in future. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. ok, although I do not fully understand the implemented it now. |
||
passed = len(aGeom.asPolyline()) > 2 | ||
|
||
if not passed: | ||
passed = (len(aGeom.asPolyline()) == 2 and | ||
aGeom.asPolyline()[0] != aGeom.asPolyline()[1]) | ||
# sometimes splitting results in lines of zero length | ||
|
||
if passed: | ||
outFeat.setGeometry(aGeom) | ||
writer.addFeature(outFeat) | ||
|
||
progress.setPercentage(int(current * total)) | ||
|
||
if multiGeoms > 0: | ||
ProcessingLog.addToLog(ProcessingLog.LOG_INFO, | ||
self.tr('Feature geometry error: %s input features ignored due to multi-geometry.') % str(multiGeoms)) | ||
|
||
del writer |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,32 @@ | ||
<GMLFeatureClassList> | ||
<GMLFeatureClass> | ||
<Name>polygons_split_with_lines</Name> | ||
<ElementPath>polygons_split_with_lines</ElementPath> | ||
<!--POLYGON--> | ||
<GeometryType>3</GeometryType> | ||
<SRSName>EPSG:4326</SRSName> | ||
<DatasetSpecificInfo> | ||
<FeatureCount>6</FeatureCount> | ||
<ExtentXMin>-1.00000</ExtentXMin> | ||
<ExtentXMax>10.00000</ExtentXMax> | ||
<ExtentYMin>-3.00000</ExtentYMin> | ||
<ExtentYMax>6.00000</ExtentYMax> | ||
</DatasetSpecificInfo> | ||
<PropertyDefn> | ||
<Name>name</Name> | ||
<ElementPath>name</ElementPath> | ||
<Type>String</Type> | ||
<Width>5</Width> | ||
</PropertyDefn> | ||
<PropertyDefn> | ||
<Name>intval</Name> | ||
<ElementPath>intval</ElementPath> | ||
<Type>Integer</Type> | ||
</PropertyDefn> | ||
<PropertyDefn> | ||
<Name>floatval</Name> | ||
<ElementPath>floatval</ElementPath> | ||
<Type>Real</Type> | ||
</PropertyDefn> | ||
</GMLFeatureClass> | ||
</GMLFeatureClassList> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"Splitted" -> "split"
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fixed