Skip to content

Commit

Permalink
Uses spatial index to select intersecting features (used selections b…
Browse files Browse the repository at this point in the history 
…efore). Tools should both run faster and with fewer map canvas glitches.

git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@13035 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
cfarmer committed Mar 10, 2010
1 parent b8855d6 commit e899240
Show file tree
Hide file tree
Showing 2 changed files with 231 additions and 272 deletions.
274 changes: 125 additions & 149 deletions python/plugins/fTools/tools/doIntersectLines.py
View file
Expand Up @@ -39,157 +39,133 @@

class Dialog(QDialog, Ui_Dialog):

def __init__(self, iface):
QDialog.__init__(self)
self.iface = iface
# Set up the user interface from Designer.
self.setupUi(self)
QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
QObject.connect(self.inLine1, SIGNAL("currentIndexChanged(QString)"), self.update1)
QObject.connect(self.inLine2, SIGNAL("currentIndexChanged(QString)"), self.update2)
self.setWindowTitle( self.tr("Line intersections") )
# populate layer list
self.progressBar.setValue(0)
mapCanvas = self.iface.mapCanvas()
for i in range(mapCanvas.layerCount()):
layer = mapCanvas.layer(i)
if layer.type() == layer.VectorLayer:
if layer.geometryType() == QGis.Line:
self.inLine1.addItem(layer.name())
self.inLine2.addItem(layer.name())
def __init__(self, iface):
QDialog.__init__(self)
self.iface = iface
# Set up the user interface from Designer.
self.setupUi(self)
QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
QObject.connect(self.inLine1, SIGNAL("currentIndexChanged(QString)"), self.update1)
QObject.connect(self.inLine2, SIGNAL("currentIndexChanged(QString)"), self.update2)
self.setWindowTitle( self.tr("Line intersections") )
# populate layer list
self.progressBar.setValue(0)
mapCanvas = self.iface.mapCanvas()
layers = ftools_utils.getLayerNames([QGis.Line])
self.inLine1.addItems(layers)
self.inLine2.addItems(layers)

def update1(self, inputLayer):
self.inField1.clear()
changedLayer = self.getVectorLayerByName(unicode(inputLayer))
changedField = self.getFieldList(changedLayer)
for i in changedField:
self.inField1.addItem(unicode(changedField[i].name()))
def update1(self, inputLayer):
self.inField1.clear()
changedLayer = ftools_utils.getVectorLayerByName(unicode(inputLayer))
changedField = ftools_utils.getFieldList(changedLayer)
for i in changedField:
self.inField1.addItem(unicode(changedField[i].name()))

def update2(self, inputLayer):
self.inField2.clear()
changedLayer = self.getVectorLayerByName(unicode(inputLayer))
changedField = self.getFieldList(changedLayer)
for i in changedField:
self.inField2.addItem(unicode(changedField[i].name()))
def update2(self, inputLayer):
self.inField2.clear()
changedLayer = ftools_utils.getVectorLayerByName(unicode(inputLayer))
changedField = ftools_utils.getFieldList(changedLayer)
for i in changedField:
self.inField2.addItem(unicode(changedField[i].name()))

def accept(self):
if self.inLine1.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input line layer") )
elif self.outShape.text() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify output shapefile") )
elif self.inLine2.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify line intersect layer") )
elif self.inField1.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input unique ID field") )
elif self.inField2.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify intersect unique ID field") )
else:
line1 = self.inLine1.currentText()
line2 = self.inLine2.currentText()
field1 = self.inField1.currentText()
field2 = self.inField2.currentText()
outPath = self.outShape.text()
if outPath.contains("\\"):
outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
else:
outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
if outName.endsWith(".shp"):
outName = outName.left(outName.length() - 4)
self.outShape.clear()
self.compute(line1, line2, field1, field2, outPath, self.progressBar)
self.progressBar.setValue(100)
addToTOC = QMessageBox.question(self, self.tr("Generate Centroids"), self.tr("Created output point shapefile:\n%1\n\nWould you like to add the new layer to the TOC?").arg( outPath ), QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
if addToTOC == QMessageBox.Yes:
self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
self.progressBar.setValue(0)
def accept(self):
if self.inLine1.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input line layer") )
elif self.outShape.text() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify output shapefile") )
elif self.inLine2.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify line intersect layer") )
elif self.inField1.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input unique ID field") )
elif self.inField2.currentText() == "":
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify intersect unique ID field") )
else:
line1 = self.inLine1.currentText()
line2 = self.inLine2.currentText()
field1 = self.inField1.currentText()
field2 = self.inField2.currentText()
outPath = self.outShape.text()
self.outShape.clear()
self.compute(line1, line2, field1, field2, outPath, self.progressBar)
self.progressBar.setValue(100)
addToTOC = QMessageBox.question(self, self.tr("Generate Centroids"), self.tr("Created output point shapefile:\n%1\n\nWould you like to add the new layer to the TOC?").arg( outPath ), QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
if addToTOC == QMessageBox.Yes:
if not ftools_utils.addShapeToCanvas( unicode( outPath ) ):
QMessageBox.warning( self, self.tr("Geoprocessing"), self.tr( "Error loading output shapefile:\n%1" ).arg( unicode( outPath ) ))
self.progressBar.setValue(0)

def outFile(self):
self.outShape.clear()
( self.shapefileName, self.encoding ) = ftools_utils.saveDialog( self )
if self.shapefileName is None or self.encoding is None:
return
self.outShape.setText( QString( self.shapefileName ) )
def outFile(self):
self.outShape.clear()
( self.shapefileName, self.encoding ) = ftools_utils.saveDialog( self )
if self.shapefileName is None or self.encoding is None:
return
self.outShape.setText( QString( self.shapefileName ) )

def compute(self, line1, line2, field1, field2, outPath, progressBar):
layer1 = self.getVectorLayerByName(line1)
layer2 = self.getVectorLayerByName(line2)
provider1 = layer1.dataProvider()
provider2 = layer2.dataProvider()
allAttrs = provider1.attributeIndexes()
provider1.select(allAttrs)
allAttrs = provider2.attributeIndexes()
provider2.select(allAttrs)
fieldList = self.getFieldList(layer1)
index1 = provider1.fieldNameIndex(field1)
field1 = fieldList[index1]
field1.setName(unicode(field1.name()) + "_1")
fieldList = self.getFieldList(layer2)
index2 = provider2.fieldNameIndex(field2)
field2 = fieldList[index2]
field2.setName(unicode(field2.name()) + "_2")
fieldList = {0:field1, 1:field2}
sRs = provider1.crs()
check = QFile(self.shapefileName)
if check.exists():
if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
return
writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, QGis.WKBPoint, sRs)
#writer = QgsVectorFileWriter(outPath, "UTF-8", fieldList, QGis.WKBPoint, sRs)
inFeat = QgsFeature()
outFeat = QgsFeature()
inGeom = QgsGeometry()
tempGeom = QgsGeometry()
start = 15.00
add = 85.00 / layer1.featureCount()
while provider1.nextFeature(inFeat):
inGeom = inFeat.geometry()
lineList = []
#(check, lineList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
# Below is a work-around for featuresInRectangle
# Which appears to have been removed for QGIS version 1.0
layer2.select(inGeom.boundingBox(), False)
lineList = layer2.selectedFeatures()
if len(lineList) > 0: check = 0
else: check = 1
if check == 0:
for i in lineList:
tempGeom = i.geometry()
tempList = []
atMap1 = inFeat.attributeMap()
atMap2 = i.attributeMap()
if inGeom.intersects(tempGeom):
tempGeom = inGeom.intersection(tempGeom)
if tempGeom.type() == QGis.Point:
if tempGeom.isMultipart():
tempList = tempGeom.asMultiPoint()
else:
tempList.append(tempGeom.asPoint())
for j in tempList:
outFeat.setGeometry(tempGeom.fromPoint(j))
outFeat.addAttribute(0, atMap1[index1])
outFeat.addAttribute(1, atMap2[index2])
writer.addFeature(outFeat)
start = start + add
progressBar.setValue(start)
del writer

def getVectorLayerByName(self, myName):
mc = self.iface.mapCanvas()
nLayers = mc.layerCount()
for l in range(nLayers):
layer = mc.layer(l)
if layer.name() == unicode(myName):
vlayer = QgsVectorLayer(unicode(layer.source()), unicode(myName), unicode(layer.dataProvider().name()))
if vlayer.isValid():
return vlayer
else:
QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Vector layer is not valid") )

def getFieldList(self, vlayer):
fProvider = vlayer.dataProvider()
feat = QgsFeature()
allAttrs = fProvider.attributeIndexes()
fProvider.select(allAttrs)
myFields = fProvider.fields()
return myFields
def compute(self, line1, line2, field1, field2, outPath, progressBar):
layer1 = ftools_utils.getVectorLayerByName(line1)
layer2 = ftools_utils.getVectorLayerByName(line2)
provider1 = layer1.dataProvider()
provider2 = layer2.dataProvider()
allAttrs = provider1.attributeIndexes()
provider1.select(allAttrs)
allAttrs = provider2.attributeIndexes()
provider2.select(allAttrs)
fieldList = ftools_utils.getFieldList(layer1)
index1 = provider1.fieldNameIndex(field1)
field1 = fieldList[index1]
field1.setName(unicode(field1.name()) + "_1")
fieldList = ftools_utils.getFieldList(layer2)
index2 = provider2.fieldNameIndex(field2)
field2 = fieldList[index2]
field2.setName(unicode(field2.name()) + "_2")
fieldList = {0:field1, 1:field2}
sRs = provider1.crs()
check = QFile(self.shapefileName)
if check.exists():
if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
return
writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, QGis.WKBPoint, sRs)
#writer = QgsVectorFileWriter(outPath, "UTF-8", fieldList, QGis.WKBPoint, sRs)
inFeat = QgsFeature()
inFeatB = QgsFeature()
outFeat = QgsFeature()
inGeom = QgsGeometry()
tempGeom = QgsGeometry()
start = 15.00
add = 85.00 / layer1.featureCount()
index = ftools_utils.createIndex( provider2 )
while provider1.nextFeature(inFeat):
inGeom = inFeat.geometry()
lineList = []
#(check, lineList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
# Below is a work-around for featuresInRectangle
# Which appears to have been removed for QGIS version 1.0
#layer2.select(inGeom.boundingBox(), False)
#lineList = layer2.selectedFeatures()
lineList = index.intersects( inGeom.boundingBox() )
if len(lineList) > 0: check = 0
else: check = 1
if check == 0:
for i in lineList:
provider2.featureAtId( int( i ), inFeatB , True, allAttrs )
tmpGeom = QgsGeometry( inFeatB.geometry() )
#tempGeom = i.geometry()
tempList = []
atMap1 = inFeat.attributeMap()
atMap2 = inFeatB.attributeMap()
if inGeom.intersects(tmpGeom):
tempGeom = inGeom.intersection(tmpGeom)
if tempGeom.type() == QGis.Point:
if tempGeom.isMultipart():
tempList = tempGeom.asMultiPoint()
else:
tempList.append(tempGeom.asPoint())
for j in tempList:
outFeat.setGeometry(tempGeom.fromPoint(j))
outFeat.addAttribute(0, atMap1[index1])
outFeat.addAttribute(1, atMap2[index2])
writer.addFeature(outFeat)
start = start + add
progressBar.setValue(start)
del writer

0 comments on commit e899240

Please sign in to comment.