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). Tool should run faster and with fewer map canvas glitches.

git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@13036 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
cfarmer committed Mar 10, 2010
1 parent e899240 commit 781a4d1
Showing 1 changed file with 171 additions and 215 deletions.
386 changes: 171 additions & 215 deletions python/plugins/fTools/tools/doSpatialJoin.py
View file
@@ -1,3 +1,4 @@
# -*- coding: utf-8 -*-
#----------------------------------------------------------- #-----------------------------------------------------------
# #
# Spatial Join # Spatial Join
Expand Down Expand Up @@ -40,222 +41,177 @@


class Dialog(QDialog, Ui_Dialog): class Dialog(QDialog, Ui_Dialog):


def __init__(self, iface): def __init__(self, iface):
QDialog.__init__(self) QDialog.__init__(self)
self.iface = iface self.iface = iface
# Set up the user interface from Designer. # Set up the user interface from Designer.
self.setupUi(self) self.setupUi(self)
QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile) QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
self.setWindowTitle( self.tr("Join attributes by location") ) self.setWindowTitle( self.tr("Join attributes by location") )
# populate layer list # populate layer list
self.progressBar.setValue(0) self.progressBar.setValue(0)
mapCanvas = self.iface.mapCanvas() mapCanvas = self.iface.mapCanvas()
for i in range(mapCanvas.layerCount()): layers = ftools_utils.getLayerNames([QGis.Point, QGis.Line, QGis.Polygon])
layer = mapCanvas.layer(i) self.inShape.addItems(layers)
if layer.type() == layer.VectorLayer: self.joinShape.addItems(layers)
self.inShape.addItem(layer.name())
self.joinShape.addItem(layer.name()) def accept(self):

if self.inShape.currentText() == "":
def accept(self): QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify target vector layer") )
if self.inShape.currentText() == "": elif self.outShape.text() == "":
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify target vector layer") ) QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify output shapefile") )
elif self.outShape.text() == "": elif self.joinShape.currentText() == "":
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify output shapefile") ) QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify join vector layer") )
elif self.joinShape.currentText() == "": elif self.rdoSummary.isChecked() and not (self.chkMean.isChecked() or self.chkSum.isChecked() or self.chkMin.isChecked() or self.chkMax.isChecked() or self.chkMean.isChecked()):
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify join vector layer") ) QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify at least one summary statistic") )
elif self.rdoSummary.isChecked() and not (self.chkMean.isChecked() or self.chkSum.isChecked() or self.chkMin.isChecked() or self.chkMax.isChecked() or self.chkMean.isChecked()): else:
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify at least one summary statistic") ) inName = self.inShape.currentText()
else: joinName = self.joinShape.currentText()
inName = self.inShape.currentText() outPath = self.outShape.text()
joinName = self.joinShape.currentText() if self.rdoSummary.isChecked():
outPath = self.outShape.text() summary = True
if self.rdoSummary.isChecked(): sumList = []
summary = True if self.chkSum.isChecked(): sumList.append("SUM")
sumList = [] if self.chkMean.isChecked(): sumList.append("MEAN")
if self.chkSum.isChecked(): sumList.append("SUM") if self.chkMin.isChecked(): sumList.append("MIN")
if self.chkMean.isChecked(): sumList.append("MEAN") if self.chkMax.isChecked(): sumList.append("MAX")
if self.chkMin.isChecked(): sumList.append("MIN") else:
if self.chkMax.isChecked(): sumList.append("MAX") summary = False
else: sumList = ["all"]
summary = False if self.rdoKeep.isChecked(): keep = True
sumList = ["all"] else: keep = False
if self.rdoKeep.isChecked(): keep = True if outPath.contains("\\"):
else: keep = False outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
if outPath.contains("\\"): else:
outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1) outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
else: if outName.endsWith(".shp"):
outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1) outName = outName.left(outName.length() - 4)
if outName.endsWith(".shp"): self.compute(inName, joinName, outPath, summary, sumList, keep, self.progressBar)
outName = outName.left(outName.length() - 4) self.outShape.clear()
self.compute(inName, joinName, outPath, summary, sumList, keep, self.progressBar) addToTOC = QMessageBox.question(self, self.tr("Spatial Join"), self.tr("Created output shapefile:\n%1\n\nWould you like to add the new layer to the TOC?").arg(unicode(outPath)), QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
self.outShape.clear() if addToTOC == QMessageBox.Yes:
addToTOC = QMessageBox.question(self, self.tr("Spatial Join"), self.tr("Created output shapefile:\n%1\n\nWould you like to add the new layer to the TOC?").arg(unicode(outPath)), QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton) self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
if addToTOC == QMessageBox.Yes: QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr") self.progressBar.setValue(0)
QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
self.progressBar.setValue(0)


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


def compute(self, inName, joinName, outName, summary, sumList, keep, progressBar): def compute(self, inName, joinName, outName, summary, sumList, keep, progressBar):
layer1 = self.getVectorLayerByName(inName) layer1 = ftools_utils.getVectorLayerByName(inName)
provider1 = layer1.dataProvider() provider1 = layer1.dataProvider()
allAttrs = provider1.attributeIndexes() allAttrs = provider1.attributeIndexes()
provider1.select(allAttrs) provider1.select(allAttrs)
fieldList1 = self.getFieldList(layer1).values() fieldList1 = ftools_utils.getFieldList(layer1).values()


layer2 = self.getVectorLayerByName(joinName) layer2 = ftools_utils.getVectorLayerByName(joinName)
provider2 = layer2.dataProvider() provider2 = layer2.dataProvider()
allAttrs = provider2.attributeIndexes() allAttrs = provider2.attributeIndexes()
provider2.select(allAttrs) provider2.select(allAttrs)
fieldList2 = self.getFieldList(layer2) fieldList2 = ftools_utils.getFieldList(layer2)
fieldList = [] fieldList = []
if provider1.crs() <> provider2.crs(): if provider1.crs() <> provider2.crs():
QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results.")) QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results."))
if not summary: if not summary:
fieldList2 = self.testForUniqueness(fieldList1, fieldList2.values()) fieldList2 = ftools_utils.testForUniqueness(fieldList1, fieldList2.values())
seq = range(0, len(fieldList1) + len(fieldList2)) seq = range(0, len(fieldList1) + len(fieldList2))
fieldList1.extend(fieldList2) fieldList1.extend(fieldList2)
fieldList1 = dict(zip(seq, fieldList1)) fieldList1 = dict(zip(seq, fieldList1))
else: else:
numFields = {} numFields = {}
for j in fieldList2.keys(): for j in fieldList2.keys():
if fieldList2[j].type() == QVariant.Int or fieldList2[j].type() == QVariant.Double: if fieldList2[j].type() == QVariant.Int or fieldList2[j].type() == QVariant.Double:
numFields[j] = [] numFields[j] = []
for i in sumList: for i in sumList:
field = QgsField(i + unicode(fieldList2[j].name()), QVariant.Double, "real", 24, 16, self.tr("Summary field") ) field = QgsField(i + unicode(fieldList2[j].name()), QVariant.Double, "real", 24, 16, self.tr("Summary field") )
fieldList.append(field) fieldList.append(field)
field = QgsField("COUNT", QVariant.Double, "real", 24, 16, self.tr("Summary field") ) field = QgsField("COUNT", QVariant.Double, "real", 24, 16, self.tr("Summary field") )
fieldList.append(field) fieldList.append(field)
fieldList2 = self.testForUniqueness(fieldList1, fieldList) fieldList2 = ftools_utils.testForUniqueness(fieldList1, fieldList)
fieldList1.extend(fieldList) fieldList1.extend(fieldList)
seq = range(0, len(fieldList1)) seq = range(0, len(fieldList1))
fieldList1 = dict(zip(seq, fieldList1)) fieldList1 = dict(zip(seq, fieldList1))


sRs = provider1.crs() sRs = provider1.crs()
progressBar.setValue(13) progressBar.setValue(13)
check = QFile(self.shapefileName) check = QFile(self.shapefileName)
if check.exists(): if check.exists():
if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName): if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
return return
writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList1, provider1.geometryType(), sRs) writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList1, provider1.geometryType(), sRs)
#writer = QgsVectorFileWriter(outName, "UTF-8", fieldList1, provider1.geometryType(), sRs) #writer = QgsVectorFileWriter(outName, "UTF-8", fieldList1, provider1.geometryType(), sRs)
inFeat = QgsFeature() inFeat = QgsFeature()
outFeat = QgsFeature() outFeat = QgsFeature()
joinFeat = QgsFeature() inFeatB = QgsFeature()
inGeom = QgsGeometry() inGeom = QgsGeometry()
progressBar.setValue(15) progressBar.setValue(15)
start = 15.00 start = 15.00
add = 85.00 / provider1.featureCount() add = 85.00 / provider1.featureCount()
provider1.rewind() provider1.rewind()

index = ftools_utils.createIndex(provider2)
while provider1.nextFeature(inFeat): while provider1.nextFeature(inFeat):
inGeom = inFeat.geometry() inGeom = inFeat.geometry()
atMap1 = inFeat.attributeMap() atMap1 = inFeat.attributeMap()
outFeat.setGeometry(inGeom) outFeat.setGeometry(inGeom)
none = True none = True
joinList = [] joinList = []
if inGeom.type() == QGis.Point: if inGeom.type() == QGis.Point:
#(check, joinList) = layer2.featuresInRectangle(inGeom.buffer(10,2).boundingBox(), True, True) #(check, joinList) = layer2.featuresInRectangle(inGeom.buffer(10,2).boundingBox(), True, True)
layer2.select(inGeom.buffer(10,2).boundingBox(), False) #layer2.select(inGeom.buffer(10,2).boundingBox(), False)
joinList = layer2.selectedFeatures() #joinList = layer2.selectedFeatures()
if len(joinList) > 0: check = 0 joinList = index.intersects( inGeom.buffer(10,2).boundingBox() )
else: check = 1 if len(joinList) > 0: check = 0
else: else: check = 1
#(check, joinList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True) else:
layer2.select(inGeom.boundingBox(), False) #(check, joinList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
joinList = layer2.selectedFeatures() #layer2.select(inGeom.boundingBox(), False)
if len(joinList) > 0: check = 0 #joinList = layer2.selectedFeatures()
else: check = 1 joinList = index.intersects( inGeom.boundingBox() )
if check == 0: if len(joinList) > 0: check = 0
count = 0 else: check = 1
for i in joinList: if check == 0:
tempGeom = i.geometry() count = 0
if inGeom.intersects(tempGeom): for i in joinList:
count = count + 1 #tempGeom = i.geometry()
none = False provider2.featureAtId(int(i), inFeatB , True, allAttrs)
atMap2 = i.attributeMap() tmpGeom = QgsGeometry( inFeatB.geometry() )
if not summary: if inGeom.intersects(tmpGeom):
atMap = atMap1.values() count = count + 1
atMap2 = atMap2.values() none = False
atMap.extend(atMap2) atMap2 = inFeatB.attributeMap()
atMap = dict(zip(seq, atMap)) if not summary:
break atMap = atMap1.values()
else: atMap2 = atMap2.values()
for j in numFields.keys(): atMap.extend(atMap2)
numFields[j].append(atMap2[j].toDouble()[0]) atMap = dict(zip(seq, atMap))
if summary and not none: break
atMap = atMap1.values() else:
for j in numFields.keys(): for j in numFields.keys():
for k in sumList: numFields[j].append(atMap2[j].toDouble()[0])
if k == "SUM": atMap.append(QVariant(sum(numFields[j]))) if summary and not none:
elif k == "MEAN": atMap.append(QVariant(sum(numFields[j]) / count)) atMap = atMap1.values()
elif k == "MIN": atMap.append(QVariant(min(numFields[j]))) for j in numFields.keys():
else: atMap.append(QVariant(max(numFields[j]))) for k in sumList:
numFields[j] = [] if k == "SUM": atMap.append(QVariant(sum(numFields[j])))
atMap.append(QVariant(count)) elif k == "MEAN": atMap.append(QVariant(sum(numFields[j]) / count))
atMap = dict(zip(seq, atMap)) elif k == "MIN": atMap.append(QVariant(min(numFields[j])))
if none: else: atMap.append(QVariant(max(numFields[j])))
outFeat.setAttributeMap(atMap1) numFields[j] = []
else: atMap.append(QVariant(count))
outFeat.setAttributeMap(atMap) atMap = dict(zip(seq, atMap))
if keep: # keep all records if none:
writer.addFeature(outFeat) outFeat.setAttributeMap(atMap1)
else: # keep only matching records else:
if not none: outFeat.setAttributeMap(atMap)
writer.addFeature(outFeat) if keep: # keep all records
start = start + add writer.addFeature(outFeat)
progressBar.setValue(start) else: # keep only matching records
del writer if not none:

writer.addFeature(outFeat)
def testForUniqueness(self, fieldList1, fieldList2): start = start + add
changed = True progressBar.setValue(start)
while changed: del writer
changed = False
for i in fieldList1:
for j in fieldList2:
if j.name() == i.name():
j = self.createUniqueFieldName(j)
changed = True
return fieldList2

def createUniqueFieldName(self, field):
check = field.name().right(2)
if check.startsWith("_"):
(val,test) = check.right(1).toInt()
if test:
if val < 2:
val = 2
else:
val = val + 1
field.setName(field.name().left(len(field.name())-1) + unicode(val))
else:
field.setName(field.name() + "_2")
else:
field.setName(field.name() + "_2")
return field

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("Spatial Join"), 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

0 comments on commit 781a4d1

Please sign in to comment.