Skip to content

Commit 63c3671

Browse files
author
cfarmer
committedMar 10, 2010
Uses spatial index to select intersecting features (used selections before). Tool should run faster and with fewer map canvas glitches.
git-svn-id: http://svn.osgeo.org/qgis/trunk@13036 c8812cc2-4d05-0410-92ff-de0c093fc19c
·
ltr-3_40Before-merge-new_vector_api
1 parent cd23b74 commit 63c3671

File tree

1 file changed

+171
-215
lines changed

1 file changed

+171
-215
lines changed
 

‎python/plugins/fTools/tools/doSpatialJoin.py

Lines changed: 171 additions & 215 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
# -*- coding: utf-8 -*-
12
#-----------------------------------------------------------
23
#
34
# Spatial Join
@@ -40,222 +41,177 @@
4041

4142
class Dialog(QDialog, Ui_Dialog):
4243

43-
def __init__(self, iface):
44-
QDialog.__init__(self)
45-
self.iface = iface
46-
# Set up the user interface from Designer.
47-
self.setupUi(self)
48-
QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
49-
self.setWindowTitle( self.tr("Join attributes by location") )
50-
# populate layer list
51-
self.progressBar.setValue(0)
52-
mapCanvas = self.iface.mapCanvas()
53-
for i in range(mapCanvas.layerCount()):
54-
layer = mapCanvas.layer(i)
55-
if layer.type() == layer.VectorLayer:
56-
self.inShape.addItem(layer.name())
57-
self.joinShape.addItem(layer.name())
58-
59-
def accept(self):
60-
if self.inShape.currentText() == "":
61-
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify target vector layer") )
62-
elif self.outShape.text() == "":
63-
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify output shapefile") )
64-
elif self.joinShape.currentText() == "":
65-
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify join vector layer") )
66-
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()):
67-
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify at least one summary statistic") )
68-
else:
69-
inName = self.inShape.currentText()
70-
joinName = self.joinShape.currentText()
71-
outPath = self.outShape.text()
72-
if self.rdoSummary.isChecked():
73-
summary = True
74-
sumList = []
75-
if self.chkSum.isChecked(): sumList.append("SUM")
76-
if self.chkMean.isChecked(): sumList.append("MEAN")
77-
if self.chkMin.isChecked(): sumList.append("MIN")
78-
if self.chkMax.isChecked(): sumList.append("MAX")
79-
else:
80-
summary = False
81-
sumList = ["all"]
82-
if self.rdoKeep.isChecked(): keep = True
83-
else: keep = False
84-
if outPath.contains("\\"):
85-
outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
86-
else:
87-
outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
88-
if outName.endsWith(".shp"):
89-
outName = outName.left(outName.length() - 4)
90-
self.compute(inName, joinName, outPath, summary, sumList, keep, self.progressBar)
91-
self.outShape.clear()
92-
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)
93-
if addToTOC == QMessageBox.Yes:
94-
self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
95-
QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
96-
self.progressBar.setValue(0)
44+
def __init__(self, iface):
45+
QDialog.__init__(self)
46+
self.iface = iface
47+
# Set up the user interface from Designer.
48+
self.setupUi(self)
49+
QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
50+
self.setWindowTitle( self.tr("Join attributes by location") )
51+
# populate layer list
52+
self.progressBar.setValue(0)
53+
mapCanvas = self.iface.mapCanvas()
54+
layers = ftools_utils.getLayerNames([QGis.Point, QGis.Line, QGis.Polygon])
55+
self.inShape.addItems(layers)
56+
self.joinShape.addItems(layers)
57+
58+
def accept(self):
59+
if self.inShape.currentText() == "":
60+
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify target vector layer") )
61+
elif self.outShape.text() == "":
62+
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify output shapefile") )
63+
elif self.joinShape.currentText() == "":
64+
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify join vector layer") )
65+
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()):
66+
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify at least one summary statistic") )
67+
else:
68+
inName = self.inShape.currentText()
69+
joinName = self.joinShape.currentText()
70+
outPath = self.outShape.text()
71+
if self.rdoSummary.isChecked():
72+
summary = True
73+
sumList = []
74+
if self.chkSum.isChecked(): sumList.append("SUM")
75+
if self.chkMean.isChecked(): sumList.append("MEAN")
76+
if self.chkMin.isChecked(): sumList.append("MIN")
77+
if self.chkMax.isChecked(): sumList.append("MAX")
78+
else:
79+
summary = False
80+
sumList = ["all"]
81+
if self.rdoKeep.isChecked(): keep = True
82+
else: keep = False
83+
if outPath.contains("\\"):
84+
outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
85+
else:
86+
outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
87+
if outName.endsWith(".shp"):
88+
outName = outName.left(outName.length() - 4)
89+
self.compute(inName, joinName, outPath, summary, sumList, keep, self.progressBar)
90+
self.outShape.clear()
91+
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)
92+
if addToTOC == QMessageBox.Yes:
93+
self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
94+
QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
95+
self.progressBar.setValue(0)
9796

98-
def outFile(self):
99-
self.outShape.clear()
100-
( self.shapefileName, self.encoding ) = ftools_utils.saveDialog( self )
101-
if self.shapefileName is None or self.encoding is None:
102-
return
103-
self.outShape.setText( QString( self.shapefileName ) )
97+
def outFile(self):
98+
self.outShape.clear()
99+
( self.shapefileName, self.encoding ) = ftools_utils.saveDialog( self )
100+
if self.shapefileName is None or self.encoding is None:
101+
return
102+
self.outShape.setText( QString( self.shapefileName ) )
104103

105-
def compute(self, inName, joinName, outName, summary, sumList, keep, progressBar):
106-
layer1 = self.getVectorLayerByName(inName)
107-
provider1 = layer1.dataProvider()
108-
allAttrs = provider1.attributeIndexes()
109-
provider1.select(allAttrs)
110-
fieldList1 = self.getFieldList(layer1).values()
104+
def compute(self, inName, joinName, outName, summary, sumList, keep, progressBar):
105+
layer1 = ftools_utils.getVectorLayerByName(inName)
106+
provider1 = layer1.dataProvider()
107+
allAttrs = provider1.attributeIndexes()
108+
provider1.select(allAttrs)
109+
fieldList1 = ftools_utils.getFieldList(layer1).values()
111110

112-
layer2 = self.getVectorLayerByName(joinName)
113-
provider2 = layer2.dataProvider()
114-
allAttrs = provider2.attributeIndexes()
115-
provider2.select(allAttrs)
116-
fieldList2 = self.getFieldList(layer2)
117-
fieldList = []
118-
if provider1.crs() <> provider2.crs():
119-
QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results."))
120-
if not summary:
121-
fieldList2 = self.testForUniqueness(fieldList1, fieldList2.values())
122-
seq = range(0, len(fieldList1) + len(fieldList2))
123-
fieldList1.extend(fieldList2)
124-
fieldList1 = dict(zip(seq, fieldList1))
125-
else:
126-
numFields = {}
127-
for j in fieldList2.keys():
128-
if fieldList2[j].type() == QVariant.Int or fieldList2[j].type() == QVariant.Double:
129-
numFields[j] = []
130-
for i in sumList:
131-
field = QgsField(i + unicode(fieldList2[j].name()), QVariant.Double, "real", 24, 16, self.tr("Summary field") )
132-
fieldList.append(field)
133-
field = QgsField("COUNT", QVariant.Double, "real", 24, 16, self.tr("Summary field") )
134-
fieldList.append(field)
135-
fieldList2 = self.testForUniqueness(fieldList1, fieldList)
136-
fieldList1.extend(fieldList)
137-
seq = range(0, len(fieldList1))
138-
fieldList1 = dict(zip(seq, fieldList1))
111+
layer2 = ftools_utils.getVectorLayerByName(joinName)
112+
provider2 = layer2.dataProvider()
113+
allAttrs = provider2.attributeIndexes()
114+
provider2.select(allAttrs)
115+
fieldList2 = ftools_utils.getFieldList(layer2)
116+
fieldList = []
117+
if provider1.crs() <> provider2.crs():
118+
QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results."))
119+
if not summary:
120+
fieldList2 = ftools_utils.testForUniqueness(fieldList1, fieldList2.values())
121+
seq = range(0, len(fieldList1) + len(fieldList2))
122+
fieldList1.extend(fieldList2)
123+
fieldList1 = dict(zip(seq, fieldList1))
124+
else:
125+
numFields = {}
126+
for j in fieldList2.keys():
127+
if fieldList2[j].type() == QVariant.Int or fieldList2[j].type() == QVariant.Double:
128+
numFields[j] = []
129+
for i in sumList:
130+
field = QgsField(i + unicode(fieldList2[j].name()), QVariant.Double, "real", 24, 16, self.tr("Summary field") )
131+
fieldList.append(field)
132+
field = QgsField("COUNT", QVariant.Double, "real", 24, 16, self.tr("Summary field") )
133+
fieldList.append(field)
134+
fieldList2 = ftools_utils.testForUniqueness(fieldList1, fieldList)
135+
fieldList1.extend(fieldList)
136+
seq = range(0, len(fieldList1))
137+
fieldList1 = dict(zip(seq, fieldList1))
139138

140-
sRs = provider1.crs()
141-
progressBar.setValue(13)
142-
check = QFile(self.shapefileName)
143-
if check.exists():
144-
if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
145-
return
146-
writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList1, provider1.geometryType(), sRs)
147-
#writer = QgsVectorFileWriter(outName, "UTF-8", fieldList1, provider1.geometryType(), sRs)
148-
inFeat = QgsFeature()
149-
outFeat = QgsFeature()
150-
joinFeat = QgsFeature()
151-
inGeom = QgsGeometry()
152-
progressBar.setValue(15)
153-
start = 15.00
154-
add = 85.00 / provider1.featureCount()
155-
provider1.rewind()
156-
157-
while provider1.nextFeature(inFeat):
158-
inGeom = inFeat.geometry()
159-
atMap1 = inFeat.attributeMap()
160-
outFeat.setGeometry(inGeom)
161-
none = True
162-
joinList = []
163-
if inGeom.type() == QGis.Point:
164-
#(check, joinList) = layer2.featuresInRectangle(inGeom.buffer(10,2).boundingBox(), True, True)
165-
layer2.select(inGeom.buffer(10,2).boundingBox(), False)
166-
joinList = layer2.selectedFeatures()
167-
if len(joinList) > 0: check = 0
168-
else: check = 1
169-
else:
170-
#(check, joinList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
171-
layer2.select(inGeom.boundingBox(), False)
172-
joinList = layer2.selectedFeatures()
173-
if len(joinList) > 0: check = 0
174-
else: check = 1
175-
if check == 0:
176-
count = 0
177-
for i in joinList:
178-
tempGeom = i.geometry()
179-
if inGeom.intersects(tempGeom):
180-
count = count + 1
181-
none = False
182-
atMap2 = i.attributeMap()
183-
if not summary:
184-
atMap = atMap1.values()
185-
atMap2 = atMap2.values()
186-
atMap.extend(atMap2)
187-
atMap = dict(zip(seq, atMap))
188-
break
189-
else:
190-
for j in numFields.keys():
191-
numFields[j].append(atMap2[j].toDouble()[0])
192-
if summary and not none:
193-
atMap = atMap1.values()
194-
for j in numFields.keys():
195-
for k in sumList:
196-
if k == "SUM": atMap.append(QVariant(sum(numFields[j])))
197-
elif k == "MEAN": atMap.append(QVariant(sum(numFields[j]) / count))
198-
elif k == "MIN": atMap.append(QVariant(min(numFields[j])))
199-
else: atMap.append(QVariant(max(numFields[j])))
200-
numFields[j] = []
201-
atMap.append(QVariant(count))
202-
atMap = dict(zip(seq, atMap))
203-
if none:
204-
outFeat.setAttributeMap(atMap1)
205-
else:
206-
outFeat.setAttributeMap(atMap)
207-
if keep: # keep all records
208-
writer.addFeature(outFeat)
209-
else: # keep only matching records
210-
if not none:
211-
writer.addFeature(outFeat)
212-
start = start + add
213-
progressBar.setValue(start)
214-
del writer
215-
216-
def testForUniqueness(self, fieldList1, fieldList2):
217-
changed = True
218-
while changed:
219-
changed = False
220-
for i in fieldList1:
221-
for j in fieldList2:
222-
if j.name() == i.name():
223-
j = self.createUniqueFieldName(j)
224-
changed = True
225-
return fieldList2
226-
227-
def createUniqueFieldName(self, field):
228-
check = field.name().right(2)
229-
if check.startsWith("_"):
230-
(val,test) = check.right(1).toInt()
231-
if test:
232-
if val < 2:
233-
val = 2
234-
else:
235-
val = val + 1
236-
field.setName(field.name().left(len(field.name())-1) + unicode(val))
237-
else:
238-
field.setName(field.name() + "_2")
239-
else:
240-
field.setName(field.name() + "_2")
241-
return field
242-
243-
def getVectorLayerByName(self, myName):
244-
mc = self.iface.mapCanvas()
245-
nLayers = mc.layerCount()
246-
for l in range(nLayers):
247-
layer = mc.layer(l)
248-
if layer.name() == unicode(myName):
249-
vlayer = QgsVectorLayer(unicode(layer.source()), unicode(myName), unicode(layer.dataProvider().name()))
250-
if vlayer.isValid():
251-
return vlayer
252-
else:
253-
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Vector layer is not valid"))
254-
255-
def getFieldList(self, vlayer):
256-
fProvider = vlayer.dataProvider()
257-
feat = QgsFeature()
258-
allAttrs = fProvider.attributeIndexes()
259-
fProvider.select(allAttrs)
260-
myFields = fProvider.fields()
261-
return myFields
139+
sRs = provider1.crs()
140+
progressBar.setValue(13)
141+
check = QFile(self.shapefileName)
142+
if check.exists():
143+
if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
144+
return
145+
writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList1, provider1.geometryType(), sRs)
146+
#writer = QgsVectorFileWriter(outName, "UTF-8", fieldList1, provider1.geometryType(), sRs)
147+
inFeat = QgsFeature()
148+
outFeat = QgsFeature()
149+
inFeatB = QgsFeature()
150+
inGeom = QgsGeometry()
151+
progressBar.setValue(15)
152+
start = 15.00
153+
add = 85.00 / provider1.featureCount()
154+
provider1.rewind()
155+
index = ftools_utils.createIndex(provider2)
156+
while provider1.nextFeature(inFeat):
157+
inGeom = inFeat.geometry()
158+
atMap1 = inFeat.attributeMap()
159+
outFeat.setGeometry(inGeom)
160+
none = True
161+
joinList = []
162+
if inGeom.type() == QGis.Point:
163+
#(check, joinList) = layer2.featuresInRectangle(inGeom.buffer(10,2).boundingBox(), True, True)
164+
#layer2.select(inGeom.buffer(10,2).boundingBox(), False)
165+
#joinList = layer2.selectedFeatures()
166+
joinList = index.intersects( inGeom.buffer(10,2).boundingBox() )
167+
if len(joinList) > 0: check = 0
168+
else: check = 1
169+
else:
170+
#(check, joinList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
171+
#layer2.select(inGeom.boundingBox(), False)
172+
#joinList = layer2.selectedFeatures()
173+
joinList = index.intersects( inGeom.boundingBox() )
174+
if len(joinList) > 0: check = 0
175+
else: check = 1
176+
if check == 0:
177+
count = 0
178+
for i in joinList:
179+
#tempGeom = i.geometry()
180+
provider2.featureAtId(int(i), inFeatB , True, allAttrs)
181+
tmpGeom = QgsGeometry( inFeatB.geometry() )
182+
if inGeom.intersects(tmpGeom):
183+
count = count + 1
184+
none = False
185+
atMap2 = inFeatB.attributeMap()
186+
if not summary:
187+
atMap = atMap1.values()
188+
atMap2 = atMap2.values()
189+
atMap.extend(atMap2)
190+
atMap = dict(zip(seq, atMap))
191+
break
192+
else:
193+
for j in numFields.keys():
194+
numFields[j].append(atMap2[j].toDouble()[0])
195+
if summary and not none:
196+
atMap = atMap1.values()
197+
for j in numFields.keys():
198+
for k in sumList:
199+
if k == "SUM": atMap.append(QVariant(sum(numFields[j])))
200+
elif k == "MEAN": atMap.append(QVariant(sum(numFields[j]) / count))
201+
elif k == "MIN": atMap.append(QVariant(min(numFields[j])))
202+
else: atMap.append(QVariant(max(numFields[j])))
203+
numFields[j] = []
204+
atMap.append(QVariant(count))
205+
atMap = dict(zip(seq, atMap))
206+
if none:
207+
outFeat.setAttributeMap(atMap1)
208+
else:
209+
outFeat.setAttributeMap(atMap)
210+
if keep: # keep all records
211+
writer.addFeature(outFeat)
212+
else: # keep only matching records
213+
if not none:
214+
writer.addFeature(outFeat)
215+
start = start + add
216+
progressBar.setValue(start)
217+
del writer

0 commit comments

Comments
 (0)
Please sign in to comment.