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

Uses spatial index to select intersecting features (used selections before). Tool should run faster and with fewer map canvas glitches.

Browse Source 
git-svn-id: http://svn.osgeo.org/qgis/trunk@13036 c8812cc2-4d05-0410-92ff-de0c093fc19c
This commit is contained in:
cfarmer 2010-03-10 00:49:23 +00:00
parent cd23b747a9
commit 63c36711a0
1 changed files with 171 additions and 215 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

386
python/plugins/fTools/tools/doSpatialJoin.py
View File

@ -1,3 +1,4 @@
# -*- coding: utf-8 -*-
#----------------------------------------------------------- #-----------------------------------------------------------
# #
# Spatial Join # Spatial Join
@ -40,222 +41,177 @@ from ui_frmSpatialJoin import Ui_Dialog
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