← Back to branch summary

~ubuntu-branches/ubuntu/wily/qgis/wily

« back to all changes in this revision

Viewing changes to python/plugins/fTools/tools/doSpatialJoin.py

  • Committer: Bazaar Package Importer
  • Author(s): Johan Van de Wauw
  • Date: 2010-07-11 20:23:24 UTC
  • mfrom: (3.1.4 squeeze)
  • Revision ID: james.westby@ubuntu.com-20100711202324-5ktghxa7hracohmr
Tags: 1.4.0+12730-3ubuntu1
https://launchpad.net/bugs/540941
* Merge from Debian unstable (LP: #540941).
* Fix compilation issues with QT 4.7
* Add build-depends on libqt4-webkit-dev 

Show diffs side-by-side

added added

removed removed

Lines of Context:
python/plugins/fTools/tools/doSpatialJoin.py
 
1
#-----------------------------------------------------------
 
2
# 
 
3
# Spatial Join
 
4
#
 
5
# Performing an attribute join between vector layers based on spatial
 
6
# relationship. Attributes from one vector layer are appended to the
 
7
# attribute table of another layer based on the relative locations of
 
8
# the features in the two layers.
 
9
#
 
10
# Copyright (C) 2008  Carson Farmer
 
11
#
 
12
# EMAIL: carson.farmer (at) gmail.com
 
13
# WEB  : www.geog.uvic.ca/spar/carson
 
14
#
 
15
#-----------------------------------------------------------
 
16
# 
 
17
# licensed under the terms of GNU GPL 2
 
18
# 
 
19
# This program is free software; you can redistribute it and/or modify
 
20
# it under the terms of the GNU General Public License as published by
 
21
# the Free Software Foundation; either version 2 of the License, or
 
22
# (at your option) any later version.
 
23
# 
 
24
# This program is distributed in the hope that it will be useful,
 
25
# but WITHOUT ANY WARRANTY; without even the implied warranty of
 
26
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
27
# GNU General Public License for more details.
 
28
# 
 
29
# You should have received a copy of the GNU General Public License along
 
30
# with this program; if not, write to the Free Software Foundation, Inc.,
 
31
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
32
# 
 
33
#---------------------------------------------------------------------
 
34
 
 
35
from PyQt4.QtCore import *
 
36
from PyQt4.QtGui import *
 
37
import ftools_utils
 
38
from qgis.core import *
 
39
from ui_frmSpatialJoin import Ui_Dialog
 
40
 
 
41
class Dialog(QDialog, Ui_Dialog):
 
42
 
 
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("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, "Spatial Join", "Please specify target vector layer")
 
62
                elif self.outShape.text() == "":
 
63
                        QMessageBox.information(self, "Spatial Join", "Please specify output shapefile")
 
64
                elif self.joinShape.currentText() == "":
 
65
                        QMessageBox.information(self, "Spatial Join", "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, "Spatial Join", "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, "Spatial Join", "Created output Shapefile:\n" + outPath 
 
93
                        + "\n\nWould you like to add the new layer to the TOC?", QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
 
94
                        if addToTOC == QMessageBox.Yes:
 
95
                                self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
 
96
                                QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
 
97
                self.progressBar.setValue(0)
 
98
 
 
99
        def outFile(self):
 
100
                self.outShape.clear()
 
101
                ( self.shapefileName, self.encoding ) = ftools_utils.saveDialog( self )
 
102
                if self.shapefileName is None or self.encoding is None:
 
103
                        return
 
104
                self.outShape.setText( QString( self.shapefileName ) )
 
105
 
 
106
        def compute(self, inName, joinName, outName, summary, sumList, keep, progressBar):
 
107
                layer1 = self.getVectorLayerByName(inName)
 
108
                provider1 = layer1.dataProvider()
 
109
                allAttrs = provider1.attributeIndexes()
 
110
                provider1.select(allAttrs)
 
111
                fieldList1 = self.getFieldList(layer1).values()
 
112
 
 
113
                layer2 = self.getVectorLayerByName(joinName)
 
114
                provider2 = layer2.dataProvider()
 
115
                allAttrs = provider2.attributeIndexes()
 
116
                provider2.select(allAttrs)
 
117
                fieldList2 = self.getFieldList(layer2)
 
118
                fieldList = []
 
119
                if provider1.crs() <> provider2.crs():
 
120
                        QMessageBox.warning(self, "CRS warning!", "Warning: Input layers have non-matching CRS.\nThis may cause unexpected results.")
 
121
                if not summary:
 
122
                        fieldList2 = self.testForUniqueness(fieldList1, fieldList2.values())
 
123
                        seq = range(0, len(fieldList1) + len(fieldList2))
 
124
                        fieldList1.extend(fieldList2)
 
125
                        fieldList1 = dict(zip(seq, fieldList1))
 
126
                else:
 
127
                        numFields = {}
 
128
                        for j in fieldList2.keys():
 
129
                                if fieldList2[j].type() == QVariant.Int or fieldList2[j].type() == QVariant.Double:
 
130
                                        numFields[j] = []
 
131
                                        for i in sumList:
 
132
                                                field = QgsField(i + unicode(fieldList2[j].name()), QVariant.Double, "real", 24, 16, "Summary field")
 
133
                                                fieldList.append(field)
 
134
                        field = QgsField("COUNT", QVariant.Double, "real", 24, 16, "Summary field")
 
135
                        fieldList.append(field)
 
136
                        fieldList2 = self.testForUniqueness(fieldList1, fieldList)
 
137
                        fieldList1.extend(fieldList)
 
138
                        seq = range(0, len(fieldList1))
 
139
                        fieldList1 = dict(zip(seq, fieldList1))
 
140
 
 
141
                sRs = provider1.crs()
 
142
                progressBar.setValue(13)
 
143
                check = QFile(self.shapefileName)
 
144
                if check.exists():
 
145
                        if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
 
146
                                return
 
147
                writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList1, provider1.geometryType(), sRs)
 
148
                #writer = QgsVectorFileWriter(outName, "UTF-8", fieldList1, provider1.geometryType(), sRs)
 
149
                inFeat = QgsFeature()
 
150
                outFeat = QgsFeature()
 
151
                joinFeat = QgsFeature()
 
152
                inGeom = QgsGeometry()
 
153
                progressBar.setValue(15)
 
154
                start = 15.00
 
155
                add = 85.00 / provider1.featureCount()
 
156
                provider1.rewind()
 
157
 
 
158
                while provider1.nextFeature(inFeat):
 
159
                        inGeom = inFeat.geometry()
 
160
                        atMap1 = inFeat.attributeMap()
 
161
                        outFeat.setGeometry(inGeom)
 
162
                        none = True
 
163
                        joinList = []
 
164
                        if inGeom.type() == QGis.Point:
 
165
                                #(check, joinList) = layer2.featuresInRectangle(inGeom.buffer(10,2).boundingBox(), True, True)
 
166
                                layer2.select(inGeom.buffer(10,2).boundingBox(), False)
 
167
                                joinList = layer2.selectedFeatures()
 
168
                                if len(joinList) > 0: check = 0
 
169
                                else: check = 1
 
170
                        else:
 
171
                                #(check, joinList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
 
172
                                layer2.select(inGeom.boundingBox(), False)
 
173
                                joinList = layer2.selectedFeatures()
 
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
                                        if inGeom.intersects(tempGeom):
 
181
                                                count = count + 1
 
182
                                                none = False
 
183
                                                atMap2 = i.attributeMap()
 
184
                                                if not summary: 
 
185
                                                        atMap = atMap1.values()
 
186
                                                        atMap2 = atMap2.values()
 
187
                                                        atMap.extend(atMap2)
 
188
                                                        atMap = dict(zip(seq, atMap))                                   
 
189
                                                        break
 
190
                                                else:
 
191
                                                        for j in numFields.keys():
 
192
                                                                numFields[j].append(atMap2[j].toDouble()[0])
 
193
                                if summary and not none:
 
194
                                        atMap = atMap1.values()
 
195
                                        for j in numFields.keys():
 
196
                                                for k in sumList:
 
197
                                                        if k == "SUM": atMap.append(QVariant(sum(numFields[j])))
 
198
                                                        elif k == "MEAN": atMap.append(QVariant(sum(numFields[j]) / count))
 
199
                                                        elif k == "MIN": atMap.append(QVariant(min(numFields[j])))
 
200
                                                        else: atMap.append(QVariant(max(numFields[j])))
 
201
                                                numFields[j] = []
 
202
                                        atMap.append(QVariant(count))
 
203
                                        atMap = dict(zip(seq, atMap))
 
204
                        if none:
 
205
                                outFeat.setAttributeMap(atMap1)
 
206
                        else:
 
207
                                outFeat.setAttributeMap(atMap)
 
208
                        if keep: # keep all records
 
209
                                writer.addFeature(outFeat)
 
210
                        else: # keep only matching records
 
211
                                if not none:
 
212
                                        writer.addFeature(outFeat)
 
213
                        start = start + add
 
214
                        progressBar.setValue(start)
 
215
                del writer
 
216
 
 
217
        def testForUniqueness(self, fieldList1, fieldList2):
 
218
                changed = True
 
219
                while changed:
 
220
                        changed = False
 
221
                        for i in fieldList1:
 
222
                                for j in fieldList2:
 
223
                                        if j.name() == i.name():
 
224
                                                j = self.createUniqueFieldName(j)
 
225
                                                changed = True
 
226
                return fieldList2
 
227
                                
 
228
        def createUniqueFieldName(self, field):
 
229
                check = field.name().right(2)
 
230
                if check.startsWith("_"):
 
231
                        (val,test) = check.right(1).toInt()
 
232
                        if test:
 
233
                                if val < 2:
 
234
                                        val = 2
 
235
                                else:
 
236
                                        val = val + 1
 
237
                                field.setName(field.name().left(len(field.name())-1) + unicode(val))
 
238
                        else:
 
239
                                field.setName(field.name() + "_2")
 
240
                else:
 
241
                        field.setName(field.name() + "_2")
 
242
                return field
 
243
                                
 
244
        def getVectorLayerByName(self, myName):
 
245
                mc = self.iface.mapCanvas()
 
246
                nLayers = mc.layerCount()
 
247
                for l in range(nLayers):
 
248
                        layer = mc.layer(l)
 
249
                        if layer.name() == unicode(myName):
 
250
                                vlayer = QgsVectorLayer(unicode(layer.source()),  unicode(myName),  unicode(layer.dataProvider().name()))
 
251
                                if vlayer.isValid():
 
252
                                        return vlayer
 
253
                                else:
 
254
                                        QMessageBox.information(self, "Spatial Join", "Vector layer is not valid")
 
255
 
 
256
        def getFieldList(self, vlayer):
 
257
                fProvider = vlayer.dataProvider()
 
258
                feat = QgsFeature()
 
259
                allAttrs = fProvider.attributeIndexes()
 
260
                fProvider.select(allAttrs)
 
261
                myFields = fProvider.fields()
 
262
                return myFields