|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | + |
| 3 | +""" |
| 4 | +*************************************************************************** |
| 5 | + SpatialJoin.py |
| 6 | + --------------------- |
| 7 | + Date : October 2013 |
| 8 | + Copyright : (C) 2013 by Joshua Arnott |
| 9 | + Email : josh at snorfalorpagus dot net |
| 10 | +*************************************************************************** |
| 11 | +* * |
| 12 | +* This program is free software; you can redistribute it and/or modify * |
| 13 | +* it under the terms of the GNU General Public License as published by * |
| 14 | +* the Free Software Foundation; either version 2 of the License, or * |
| 15 | +* (at your option) any later version. * |
| 16 | +* * |
| 17 | +*************************************************************************** |
| 18 | +""" |
| 19 | +from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException |
| 20 | +__author__ = 'Joshua Arnott' |
| 21 | +__date__ = 'October 2013' |
| 22 | +__copyright__ = '(C) 2013, Joshua Arnott' |
| 23 | +# This will get replaced with a git SHA1 when you do a git archive |
| 24 | +__revision__ = '$Format:%H$' |
| 25 | + |
| 26 | +from processing.core.GeoAlgorithm import GeoAlgorithm |
| 27 | +from PyQt4.QtCore import * |
| 28 | +from PyQt4.QtGui import * |
| 29 | +from qgis.core import * |
| 30 | +from processing.parameters.ParameterVector import ParameterVector |
| 31 | +from processing.parameters.ParameterSelection import ParameterSelection |
| 32 | +from processing.parameters.ParameterString import ParameterString |
| 33 | +from processing.tools import dataobjects |
| 34 | +from processing.outputs.OutputVector import OutputVector |
| 35 | +from processing.tools import vector |
| 36 | +from processing.core.ProcessingLog import ProcessingLog |
| 37 | +import os |
| 38 | + |
| 39 | +def myself(L): |
| 40 | + #median computation |
| 41 | + nVal = len(L) |
| 42 | + if nVal == 1: |
| 43 | + return L[0] |
| 44 | + L.sort() |
| 45 | + #test for list length |
| 46 | + medianVal = 0 |
| 47 | + if nVal > 1: |
| 48 | + if ( nVal % 2 ) == 0: |
| 49 | + #index begin at 0 |
| 50 | + #remove 1 to index in standard median computation |
| 51 | + medianVal = 0.5 * ( (L[ (nVal) / 2 - 1]) + (L[ (nVal) / 2 ] )) |
| 52 | + else: |
| 53 | + medianVal = L[ (nVal + 1) / 2 - 1] |
| 54 | + return medianVal |
| 55 | + |
| 56 | +class SpatialJoin(GeoAlgorithm): |
| 57 | + ''' |
| 58 | + Join by location |
| 59 | + |
| 60 | + Port of the spatial join algorithm from fTools to the Processing Toolbox. |
| 61 | + ''' |
| 62 | + |
| 63 | + INPUT1 = "INPUT1" |
| 64 | + INPUT2 = "INPUT2" |
| 65 | + SUMMARY = "SUMMARY" |
| 66 | + STATS = "STATS" |
| 67 | + GEOMETRY = "GEOMETRY" |
| 68 | + KEEP = "KEEP" |
| 69 | + OUTPUT = "OUTPUT" |
| 70 | + |
| 71 | + SUMMARYS = [ |
| 72 | + 'Take attributes of the first located feature', |
| 73 | + 'Take summary of intersecting features' |
| 74 | + ] |
| 75 | + |
| 76 | + GEOMETRYS = [ |
| 77 | + 'Use geometry from target layer', |
| 78 | + 'Use geometry from joined layer (multipart if summary)' |
| 79 | + ] |
| 80 | + |
| 81 | + KEEPS = [ |
| 82 | + 'Only keep matching records', |
| 83 | + 'Keep all records (including non-matching target records)' |
| 84 | + ] |
| 85 | + |
| 86 | + #=========================================================================== |
| 87 | + # def getIcon(self): |
| 88 | + # return QIcon(os.path.dirname(__file__) + "/icons/join_location.png") |
| 89 | + #=========================================================================== |
| 90 | + |
| 91 | + def defineCharacteristics(self): |
| 92 | + self.name = "Join by location" |
| 93 | + self.group = "Vector general tools" |
| 94 | + self.addParameter(ParameterVector(SpatialJoin.INPUT1, "Target vector layer", [ParameterVector.VECTOR_TYPE_ANY])) |
| 95 | + self.addParameter(ParameterVector(SpatialJoin.INPUT2, "Join vector layer", [ParameterVector.VECTOR_TYPE_ANY])) |
| 96 | + self.addParameter(ParameterSelection(self.SUMMARY, "Attribute summary", self.SUMMARYS, 0)) |
| 97 | + self.addParameter(ParameterString(self.STATS, "Statistics for summary (comma separated)", "sum,mean,min,max,median")) |
| 98 | + self.addParameter(ParameterSelection(self.GEOMETRY, "Output geometry", self.GEOMETRYS, 0)) |
| 99 | + self.addParameter(ParameterSelection(self.KEEP, "Output table", self.KEEPS, 0)) |
| 100 | + self.addOutput(OutputVector(SpatialJoin.OUTPUT, "Output layer")) |
| 101 | + |
| 102 | + def processAlgorithm(self, progress): |
| 103 | + progress.setText('Analysing inputs...') |
| 104 | + |
| 105 | + summary = self.getParameterValue(self.SUMMARY) == 1 |
| 106 | + sumList = self.getParameterValue(self.STATS).upper().replace(' ','').split(',') |
| 107 | + use_geom = self.getParameterValue(self.GEOMETRY) |
| 108 | + keep = self.getParameterValue(self.KEEP) == 1 |
| 109 | + |
| 110 | + input1 = self.getParameterValue(self.INPUT1) |
| 111 | + layer1 = dataobjects.getObjectFromUri(input1) |
| 112 | + provider1 = layer1.dataProvider() |
| 113 | + fieldList1 = provider1.fields() |
| 114 | + |
| 115 | + input2 = self.getParameterValue(self.INPUT2) |
| 116 | + layer2 = dataobjects.getObjectFromUri(input2) |
| 117 | + provider2 = layer2.dataProvider() |
| 118 | + fieldList2 = provider2.fields() |
| 119 | + |
| 120 | + fieldList = QgsFields() |
| 121 | + if not summary: |
| 122 | + fieldList2 = vector.testForUniqueness(fieldList1, fieldList2) |
| 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 xrange(len(fieldList2)): |
| 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, 0, "Summary field" ) |
| 135 | + fieldList.append(field) |
| 136 | + fieldList2 = vector.testForUniqueness(fieldList1, fieldList) |
| 137 | + fieldList1.extend(fieldList) |
| 138 | + seq = range(0, len(fieldList1)) |
| 139 | + fieldList1 = dict(zip(seq, fieldList1)) |
| 140 | + |
| 141 | + progress.setPercentage(13) |
| 142 | + fields = QgsFields() |
| 143 | + for f in fieldList1.values(): |
| 144 | + fields.append(f) |
| 145 | + output = self.getOutputFromName(self.OUTPUT) |
| 146 | + |
| 147 | + if use_geom == 0: |
| 148 | + # from target layer |
| 149 | + crs = provider1.crs() |
| 150 | + geometry_type = provider1.geometryType() |
| 151 | + else: |
| 152 | + # from joined layer |
| 153 | + crs = provider2.crs() |
| 154 | + if summary: |
| 155 | + geometry_type = self.singleToMultiGeom(provider2.geometryType()) |
| 156 | + else: |
| 157 | + geometry_type = provider2.geometryType() |
| 158 | + |
| 159 | + writer = output.getVectorWriter(fields, geometry_type, crs) |
| 160 | + |
| 161 | + inFeat = QgsFeature() |
| 162 | + outFeat = QgsFeature() |
| 163 | + inFeatB = QgsFeature() |
| 164 | + inGeom = QgsGeometry() |
| 165 | + |
| 166 | + progress.setPercentage(15) |
| 167 | + start = 15.00 |
| 168 | + add = 85.00 / provider1.featureCount() |
| 169 | + |
| 170 | + progress.setText('Creating spatial index...') |
| 171 | + index = vector.spatialindex(layer2) |
| 172 | + progress.setText('Processing spatial join...') |
| 173 | + fit1 = provider1.getFeatures() |
| 174 | + while fit1.nextFeature(inFeat): |
| 175 | + inGeom = inFeat.geometry() |
| 176 | + atMap1 = inFeat.attributes() |
| 177 | + if use_geom == 0: |
| 178 | + outFeat.setGeometry(inGeom) |
| 179 | + none = True |
| 180 | + joinList = [] |
| 181 | + if inGeom.type() == QGis.Point: |
| 182 | + joinList = index.intersects( inGeom.buffer(10,2).boundingBox() ) |
| 183 | + if len(joinList) > 0: check = 0 |
| 184 | + else: check = 1 |
| 185 | + else: |
| 186 | + joinList = index.intersects( inGeom.boundingBox() ) |
| 187 | + if len(joinList) > 0: check = 0 |
| 188 | + else: check = 1 |
| 189 | + if check == 0: |
| 190 | + count = 0 |
| 191 | + multi_feature = [] |
| 192 | + for i in joinList: |
| 193 | + provider2.getFeatures( QgsFeatureRequest().setFilterFid( int(i) ) ).nextFeature( inFeatB ) |
| 194 | + if inGeom.intersects(inFeatB.geometry()): |
| 195 | + count = count + 1 |
| 196 | + atMap2 = inFeatB.attributes() |
| 197 | + if not summary: |
| 198 | + # first located feature |
| 199 | + atMap = atMap1 |
| 200 | + atMap2 = atMap2 |
| 201 | + atMap.extend(atMap2) |
| 202 | + atMap = dict(zip(seq, atMap)) |
| 203 | + if use_geom == 1: |
| 204 | + outFeat.setGeometry(inFeatB.geometry()) |
| 205 | + none = False |
| 206 | + break |
| 207 | + else: |
| 208 | + for j in numFields.keys(): |
| 209 | + numFields[j].append(atMap2[j]) |
| 210 | + if use_geom == 0: |
| 211 | + if none: |
| 212 | + outFeat.setGeometry(inGeom) |
| 213 | + else: |
| 214 | + feature_list = self.extractAsMulti(inFeatB.geometry()) |
| 215 | + multi_feature.extend(feature_list) |
| 216 | + none = False |
| 217 | + if summary and not none: |
| 218 | + atMap = atMap1 |
| 219 | + for j in numFields.keys(): |
| 220 | + for k in sumList: |
| 221 | + if k == "SUM": atMap.append(sum(numFields[j])) |
| 222 | + elif k == "MEAN": atMap.append(sum(numFields[j]) / count) |
| 223 | + elif k == "MIN": atMap.append(min(numFields[j])) |
| 224 | + elif k == "MEDIAN": atMap.append(myself(numFields[j])) |
| 225 | + else: atMap.append(max(numFields[j])) |
| 226 | + numFields[j] = [] |
| 227 | + atMap.append(count) |
| 228 | + atMap = dict(zip(seq, atMap)) |
| 229 | + if use_geom == 1: |
| 230 | + outGeom = QgsGeometry(self.convertGeometry(multi_feature, geometry_type)) |
| 231 | + outFeat.setGeometry(outGeom) |
| 232 | + if none: |
| 233 | + outFeat.setAttributes(atMap1) |
| 234 | + else: |
| 235 | + outFeat.setAttributes(atMap.values()) |
| 236 | + if keep: # keep all records |
| 237 | + writer.addFeature(outFeat) |
| 238 | + else: # keep only matching records |
| 239 | + if not none: |
| 240 | + writer.addFeature(outFeat) |
| 241 | + start = start + add |
| 242 | + progress.setPercentage(start) |
| 243 | + |
| 244 | + del writer |
| 245 | + |
| 246 | + def singleToMultiGeom(self, wkbType): |
| 247 | + try: |
| 248 | + if wkbType in (QGis.WKBPoint, QGis.WKBMultiPoint, |
| 249 | + QGis.WKBPoint25D, QGis.WKBMultiPoint25D): |
| 250 | + return QGis.WKBMultiPoint |
| 251 | + elif wkbType in (QGis.WKBLineString, QGis.WKBMultiLineString, |
| 252 | + QGis.WKBMultiLineString25D, |
| 253 | + QGis.WKBLineString25D): |
| 254 | + |
| 255 | + return QGis.WKBMultiLineString |
| 256 | + elif wkbType in (QGis.WKBPolygon, QGis.WKBMultiPolygon, |
| 257 | + QGis.WKBMultiPolygon25D, QGis.WKBPolygon25D): |
| 258 | + |
| 259 | + return QGis.WKBMultiPolygon |
| 260 | + else: |
| 261 | + return QGis.WKBUnknown |
| 262 | + except Exception, err: |
| 263 | + print unicode(err) |
| 264 | + |
| 265 | + def extractAsMulti(self, geom): |
| 266 | + if geom.type() == QGis.Point: |
| 267 | + if geom.isMultipart(): |
| 268 | + return geom.asMultiPoint() |
| 269 | + else: |
| 270 | + return [geom.asPoint()] |
| 271 | + elif geom.type() == QGis.Line: |
| 272 | + if geom.isMultipart(): |
| 273 | + return geom.asMultiPolyline() |
| 274 | + else: |
| 275 | + return [geom.asPolyline()] |
| 276 | + else: |
| 277 | + if geom.isMultipart(): |
| 278 | + return geom.asMultiPolygon() |
| 279 | + else: |
| 280 | + return [geom.asPolygon()] |
| 281 | + |
| 282 | + def convertGeometry(self, geom_list, vType): |
| 283 | + if vType == QGis.Point: |
| 284 | + return QgsGeometry().fromMultiPoint(geom_list) |
| 285 | + elif vType == QGis.Line: |
| 286 | + return QgsGeometry().fromMultiPolyline(geom_list) |
| 287 | + else: |
| 288 | + return QgsGeometry().fromMultiPolygon(geom_list) |
0 commit comments