Skip to content
Permalink
Browse files

[processing] Use prepared geometries for intersects tests

Wherever possible use prepared geometry engines for intersects
type tests, as it's much faster
  • Loading branch information
nyalldawson committed Oct 17, 2016
1 parent 86368f3 commit a05b610a8bda6a845af8403d012408b4baf41e6f
@@ -253,10 +253,14 @@ def processAlgorithm(self, progress):
min = -1
selFeat = QgsFeature()

# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(geom2Eliminate.geometry())
engine.prepareGeometry()

while fit.nextFeature(selFeat):
selGeom = selFeat.geometry()

if geom2Eliminate.intersects(selGeom):
if engine.intersects(selGeom.geometry()):
# We have a candidate
iGeom = geom2Eliminate.intersection(selGeom)

@@ -89,9 +89,16 @@ def processAlgorithm(self, progress):
atMapA = inFeatA.attributes()
intersects = index.intersects(geom.boundingBox())
request = QgsFeatureRequest().setFilterFids(intersects)

engine = None
if len(intersects) > 0:
# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(geom.geometry())
engine.prepareGeometry()

for inFeatB in vlayerB.getFeatures(request):
tmpGeom = inFeatB.geometry()
if geom.intersects(tmpGeom):
if engine.intersects(tmpGeom.geometry()):
atMapB = inFeatB.attributes()
int_geom = QgsGeometry(geom.intersection(tmpGeom))
if int_geom.wkbType() == QgsWkbTypes.Unknown or QgsWkbTypes.flatType(int_geom.geometry().wkbType()) == QgsWkbTypes.GeometryCollection:
@@ -101,8 +101,12 @@ def processAlgorithm(self, progress):
hasIntersections = False
lines = spatialIndex.intersects(inGeom.boundingBox())

engine = None
if len(lines) > 0:
hasIntersections = True
# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(inGeom.geometry())
engine.prepareGeometry()

if hasIntersections:
request = QgsFeatureRequest().setFilterFids(lines)
@@ -113,7 +117,7 @@ def processAlgorithm(self, progress):
attrsA = inFeatA.attributes()
attrsB = inFeatB.attributes()

if inGeom.intersects(tmpGeom):
if engine.intersects(tmpGeom.geometry()):
tempGeom = inGeom.intersection(tmpGeom)
if tempGeom.type() == QgsWkbTypes.PointGeometry:
if tempGeom.isMultipart():
@@ -28,7 +28,7 @@
__revision__ = '$Format:%H$'

from osgeo import gdal
from qgis.core import Qgis, QgsFields, QgsField, QgsFeature, QgsPoint, QgsGeometry, QgsWkbTypes
from qgis.core import Qgis, QgsFields, QgsField, QgsFeature, QgsPoint, QgsGeometry, QgsWkbTypes, QgsPointV2
from qgis.PyQt.QtCore import QVariant
from processing.core.GeoAlgorithm import GeoAlgorithm
from processing.core.parameters import ParameterRaster
@@ -73,7 +73,6 @@ def processAlgorithm(self, progress):

outFeature = QgsFeature()
outFeature.setFields(fields)
point = QgsPoint()

fid = 0
polyId = 0
@@ -93,14 +92,19 @@ def processAlgorithm(self, progress):
(startRow, startColumn) = raster.mapToPixel(xMin, yMax, geoTransform)
(endRow, endColumn) = raster.mapToPixel(xMax, yMin, geoTransform)

# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(geom.geometry())
engine.prepareGeometry()

for row in range(startRow, endRow + 1):
for col in range(startColumn, endColumn + 1):
(x, y) = raster.pixelToMap(row, col, geoTransform)
point = QgsPointV2()
point.setX(x)
point.setY(y)

if geom.contains(point):
outFeature.setGeometry(QgsGeometry.fromPoint(point))
if engine.contains(point):
outFeature.setGeometry(QgsGeometry(point))
outFeature['id'] = fid
outFeature['poly_id'] = polyId
outFeature['point_id'] = pointId
@@ -108,6 +108,11 @@ def processAlgorithm(self, progress):
count = 0
total = 100.0 / (area / pSpacing)
y = extent.yMaximum() - inset

extent_geom = QgsGeometry.fromRect(extent)
extent_engine = QgsGeometry.createGeometryEngine(extent_geom.geometry())
extent_engine.prepareGeometry()

while y >= extent.yMinimum():
x = extent.xMinimum() + inset
while x <= extent.xMaximum():
@@ -118,7 +123,7 @@ def processAlgorithm(self, progress):
else:
geom = QgsGeometry().fromPoint(QgsPoint(x, y))

if geom.intersects(extent):
if extent_engine.intersects(geom.geometry()):
f.setAttribute('id', count)
f.setGeometry(geom)
writer.addFeature(f)
@@ -76,6 +76,12 @@ def processAlgorithm(self, progress):
inLines = [inGeom]
lines = spatialIndex.intersects(inGeom.boundingBox())

engine = None
if len(lines) > 0:
# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(inGeom.geometry())
engine.prepareGeometry()

if len(lines) > 0: # hasIntersections
splittingLines = []

@@ -88,19 +94,22 @@ def processAlgorithm(self, progress):

splitGeom = inFeatB.geometry()

if inGeom.intersects(splitGeom):
if engine.intersects(splitGeom.geometry()):
splittingLines.append(splitGeom)

if len(splittingLines) > 0:
for splitGeom in splittingLines:
splitterPList = vector.extractPoints(splitGeom)
outLines = []

split_geom_engine = QgsGeometry.createGeometryEngine(splitGeom.geometry())
split_geom_engine.prepareGeometry()

while len(inLines) > 0:
inGeom = inLines.pop()
inPoints = vector.extractPoints(inGeom)

if inGeom.intersects(splitGeom):
if split_geom_engine.intersects(inGeom.geometry()):
try:
result, newGeometries, topoTestPoints = inGeom.splitGeometry(splitterPList, False)
except:
@@ -100,14 +100,18 @@ def processAlgorithm(self, progress):
length = 0
hasIntersections = False
lines = spatialIndex.intersects(inGeom.boundingBox())
engine = None
if len(lines) > 0:
hasIntersections = True
# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(inGeom.geometry())
engine.prepareGeometry()

if hasIntersections:
request = QgsFeatureRequest().setFilterFids(lines).setSubsetOfAttributes([])
for ftLine in lineLayer.getFeatures(request):
tmpGeom = ftLine.geometry()
if inGeom.intersects(tmpGeom):
if engine.intersects(tmpGeom.geometry()):
outGeom = inGeom.intersection(tmpGeom)
length += distArea.measureLength(outGeom)
count += 1
@@ -106,13 +106,17 @@ def processAlgorithm(self, progress):
self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
else:
request = QgsFeatureRequest().setFilterFids(intersects)

engine = QgsGeometry.createGeometryEngine(geom.geometry())
engine.prepareGeometry()

for inFeatB in vlayerB.getFeatures(request):
count += 1

atMapB = inFeatB.attributes()
tmpGeom = inFeatB.geometry()

if geom.intersects(tmpGeom):
if engine.intersects(tmpGeom.geometry()):
int_geom = geom.intersection(tmpGeom)
lstIntersectingB.append(tmpGeom)

@@ -198,11 +202,16 @@ def processAlgorithm(self, progress):
self.tr('Feature geometry error: One or more output features ignored due to invalid geometry.'))
else:
request = QgsFeatureRequest().setFilterFids(intersects)

# use prepared geometries for faster intersection tests
engine = QgsGeometry.createGeometryEngine(diff_geom.geometry())
engine.prepareGeometry()

for inFeatB in vlayerA.getFeatures(request):
atMapB = inFeatB.attributes()
tmpGeom = inFeatB.geometry()

if diff_geom.intersects(tmpGeom):
if engine.intersects(tmpGeom.geometry()):
add = True
diff_geom = QgsGeometry(diff_geom.difference(tmpGeom))
if diff_geom.isGeosEmpty() or not diff_geom.isGeosValid():
@@ -0,0 +1,15 @@
<GMLFeatureClassList>
<GMLFeatureClass>
<Name>line_intersection</Name>
<ElementPath>line_intersection</ElementPath>
<GeometryType>1</GeometryType>
<SRSName>EPSG:4326</SRSName>
<DatasetSpecificInfo>
<FeatureCount>4</FeatureCount>
<ExtentXMin>3.00000</ExtentXMin>
<ExtentXMax>9.22098</ExtentXMax>
<ExtentYMin>-0.22085</ExtentYMin>
<ExtentYMax>3.22098</ExtentYMax>
</DatasetSpecificInfo>
</GMLFeatureClass>
</GMLFeatureClassList>
@@ -0,0 +1,34 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation=""
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>3</gml:X><gml:Y>-0.2208501518128253</gml:Y></gml:coord>
<gml:coord><gml:X>9.220978076442222</gml:X><gml:Y>3.220978076442222</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:line_intersection fid="lines.0">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>9.220978076442222,3.220978076442222</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:line_intersection>
</gml:featureMember>
<gml:featureMember>
<ogr:line_intersection fid="lines.2">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>3.0,2.8034726769992</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:line_intersection>
</gml:featureMember>
<gml:featureMember>
<ogr:line_intersection fid="lines.3">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>3.155795677799609,1.0</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:line_intersection>
</gml:featureMember>
<gml:featureMember>
<ogr:line_intersection fid="lines.5">
<ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>8.779149848187174,-0.220850151812825</gml:coordinates></gml:Point></ogr:geometryProperty>
</ogr:line_intersection>
</gml:featureMember>
</ogr:FeatureCollection>
@@ -0,0 +1,41 @@
<GMLFeatureClassList>
<GMLFeatureClass>
<Name>sum_line_length</Name>
<ElementPath>sum_line_length</ElementPath>
<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>
<PropertyDefn>
<Name>line_len</Name>
<ElementPath>line_len</ElementPath>
<Type>Real</Type>
</PropertyDefn>
<PropertyDefn>
<Name>line_count</Name>
<ElementPath>line_count</ElementPath>
<Type>Real</Type>
</PropertyDefn>
</GMLFeatureClass>
</GMLFeatureClassList>
@@ -0,0 +1,70 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation=""
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>-1</gml:X><gml:Y>-3</gml:Y></gml:coord>
<gml:coord><gml:X>10</gml:X><gml:Y>6</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>

<gml:featureMember>
<ogr:sum_line_length fid="polys.0">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-1,-1 -1,3 3,3 3,2 2,2 2,-1 -1,-1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>aaaaa</ogr:name>
<ogr:intval>33</ogr:intval>
<ogr:floatval>44.123456</ogr:floatval>
<ogr:line_len>6.000000000000000</ogr:line_len>
<ogr:line_count>2.000000000000000</ogr:line_count>
</ogr:sum_line_length>
</gml:featureMember>
<gml:featureMember>
<ogr:sum_line_length fid="polys.1">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>5,5 6,4 4,4 5,5</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>Aaaaa</ogr:name>
<ogr:intval>-33</ogr:intval>
<ogr:floatval>0</ogr:floatval>
<ogr:line_len>0.000000000000000</ogr:line_len>
<ogr:line_count>0.000000000000000</ogr:line_count>
</ogr:sum_line_length>
</gml:featureMember>
<gml:featureMember>
<ogr:sum_line_length fid="polys.2">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>2,5 2,6 3,6 3,5 2,5</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>bbaaa</ogr:name>
<ogr:floatval>0.123</ogr:floatval>
<ogr:line_len>0.000000000000000</ogr:line_len>
<ogr:line_count>0.000000000000000</ogr:line_count>
</ogr:sum_line_length>
</gml:featureMember>
<gml:featureMember>
<ogr:sum_line_length fid="polys.3">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>6,1 10,1 10,-3 6,-3 6,1</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs><gml:innerBoundaryIs><gml:LinearRing><gml:coordinates>7,0 7,-2 9,-2 9,0 7,0</gml:coordinates></gml:LinearRing></gml:innerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>ASDF</ogr:name>
<ogr:intval>0</ogr:intval>
<ogr:line_len>5.828427124746190</ogr:line_len>
<ogr:line_count>2.000000000000000</ogr:line_count>
</ogr:sum_line_length>
</gml:featureMember>
<gml:featureMember>
<ogr:sum_line_length fid="polys.4">
<ogr:intval>120</ogr:intval>
<ogr:floatval>-100291.43213</ogr:floatval>
<ogr:line_len>0.000000000000000</ogr:line_len>
<ogr:line_count>0.000000000000000</ogr:line_count>
</ogr:sum_line_length>
</gml:featureMember>
<gml:featureMember>
<ogr:sum_line_length fid="polys.5">
<ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>3,2 6,1 6,-3 2,-1 2,2 3,2</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
<ogr:name>elim</ogr:name>
<ogr:intval>2</ogr:intval>
<ogr:floatval>3.33</ogr:floatval>
<ogr:line_len>5.000000000000000</ogr:line_len>
<ogr:line_count>3.000000000000000</ogr:line_count>
</ogr:sum_line_length>
</gml:featureMember>
</ogr:FeatureCollection>

0 comments on commit a05b610

Please sign in to comment.
You can’t perform that action at this time.