From 8daf8c371aade774051981615d53c1275deebd58 Mon Sep 17 00:00:00 2001 From: Joshua Arnott Date: Tue, 22 Oct 2013 18:42:52 +0100 Subject: [PATCH 1/3] Ported 'Join by location' from fTools to Processing --- .../processing/algs/QGISAlgorithmProvider.py | 3 +- .../processing/algs/ftools/SpatialJoin.py | 209 ++++++++++++++++++ python/plugins/processing/tools/vector.py | 14 ++ 3 files changed, 225 insertions(+), 1 deletion(-) create mode 100644 python/plugins/processing/algs/ftools/SpatialJoin.py diff --git a/python/plugins/processing/algs/QGISAlgorithmProvider.py b/python/plugins/processing/algs/QGISAlgorithmProvider.py index c6fc541280af..f631496f0bf5 100644 --- a/python/plugins/processing/algs/QGISAlgorithmProvider.py +++ b/python/plugins/processing/algs/QGISAlgorithmProvider.py @@ -78,6 +78,7 @@ from processing.algs.ftools.DensifyGeometriesInterval import \ DensifyGeometriesInterval from processing.algs.ftools.Eliminate import Eliminate +from processing.algs.ftools.SpatialJoin import SpatialJoin from processing.algs.mmqgisx.MMQGISXAlgorithms import \ mmqgisx_delete_columns_algorithm, \ @@ -135,7 +136,7 @@ def __init__(self): Intersection(), Union(), Clip(), ExtentFromLayer(), RandomSelection(), RandomSelectionWithinSubsets(), SelectByLocation(), RandomExtract(), RandomExtractWithinSubsets(), - ExtractByLocation(), + ExtractByLocation(), SpatialJoin(), # ------ mmqgisx ------ mmqgisx_delete_columns_algorithm(), mmqgisx_delete_duplicate_geometries_algorithm(), diff --git a/python/plugins/processing/algs/ftools/SpatialJoin.py b/python/plugins/processing/algs/ftools/SpatialJoin.py new file mode 100644 index 000000000000..661ae11ac9c9 --- /dev/null +++ b/python/plugins/processing/algs/ftools/SpatialJoin.py @@ -0,0 +1,209 @@ +# -*- coding: utf-8 -*- + +""" +*************************************************************************** + SpatialJoin.py + --------------------- + Date : October 2013 + Copyright : (C) 2013 by Joshua Arnott + Email : josh at snorfalorpagus dot net +*************************************************************************** +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +*************************************************************************** +""" +from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException +__author__ = 'Joshua Arnott' +__date__ = 'October 2013' +__copyright__ = '(C) 2013, Joshua Arnott' +# This will get replaced with a git SHA1 when you do a git archive +__revision__ = '$Format:%H$' + +from processing.core.GeoAlgorithm import GeoAlgorithm +from PyQt4.QtCore import * +from PyQt4.QtGui import * +from qgis.core import * +from processing.parameters.ParameterVector import ParameterVector +from processing.parameters.ParameterSelection import ParameterSelection +from processing.parameters.ParameterString import ParameterString +from processing.tools import dataobjects +from processing.outputs.OutputVector import OutputVector +from processing.tools import vector +from processing.core.ProcessingLog import ProcessingLog +import os + +def myself(L): + #median computation + nVal = len(L) + if nVal == 1: + return L[0] + L.sort() + #test for list length + medianVal = 0 + if nVal > 1: + if ( nVal % 2 ) == 0: + #index begin at 0 + #remove 1 to index in standard median computation + medianVal = 0.5 * ( (L[ (nVal) / 2 - 1]) + (L[ (nVal) / 2 ] )) + else: + medianVal = L[ (nVal + 1) / 2 - 1] + return medianVal + +class SpatialJoin(GeoAlgorithm): + ''' + Join attributes by location + + Port of the spatial join algorithm from fTools to the Processing Toolbox. + ''' + + INPUT1 = "INPUT1" + INPUT2 = "INPUT2" + SUMMARY = "SUMMARY" + STATS = "STATS" + KEEP = "KEEP" + OUTPUT = "OUTPUT" + + SUMMARYS = [ + 'Take attributes of the first located feature', + 'Take summary of intersecting features' + ] + + KEEPS = [ + 'Only keep matching records', + 'Keep all records (including non-matching target records)' + ] + + #=========================================================================== + # def getIcon(self): + # return QIcon(os.path.dirname(__file__) + "/icons/join_location.png") + #=========================================================================== + + def defineCharacteristics(self): + self.name = "Join attributes by location" + self.group = "Vector general tools" + self.addParameter(ParameterVector(SpatialJoin.INPUT1, "Target vector layer", [ParameterVector.VECTOR_TYPE_ANY])) + self.addParameter(ParameterVector(SpatialJoin.INPUT2, "Join vector layer", [ParameterVector.VECTOR_TYPE_ANY])) + self.addParameter(ParameterSelection(self.SUMMARY, "Attribute summary", self.SUMMARYS, 0)) + self.addParameter(ParameterString(self.STATS, "Statistics for summary (comma separated)", "sum,mean,min,max,med")) + self.addParameter(ParameterSelection(self.KEEP, "Output table", self.KEEPS, 0)) + self.addOutput(OutputVector(SpatialJoin.OUTPUT, "Output layer")) + + def processAlgorithm(self, progress): + progress.setText('Analysing inputs...') + + summary = self.getParameterValue(self.SUMMARY) == 1 + sumList = self.getParameterValue(self.STATS).upper().replace(' ','').split(',') + keep = self.getParameterValue(self.KEEP) == 1 + + input1 = self.getParameterValue(self.INPUT1) + layer1 = dataobjects.getObjectFromUri(input1) + provider1 = layer1.dataProvider() + fieldList1 = provider1.fields() + + input2 = self.getParameterValue(self.INPUT2) + layer2 = dataobjects.getObjectFromUri(input2) + provider2 = layer2.dataProvider() + fieldList2 = provider2.fields() + + fieldList = QgsFields() + if not summary: + fieldList2 = vector.testForUniqueness(fieldList1, fieldList2) + seq = range(0, len(fieldList1) + len(fieldList2)) + fieldList1.extend(fieldList2) + fieldList1 = dict(zip(seq, fieldList1)) + else: + numFields = {} + for j in xrange(len(fieldList2)): + if fieldList2[j].type() == QVariant.Int or fieldList2[j].type() == QVariant.Double: + numFields[j] = [] + for i in sumList: + field = QgsField(i + unicode(fieldList2[j].name()), QVariant.Double, "real", 24, 16, "Summary field" ) + fieldList.append(field) + field = QgsField("COUNT", QVariant.Double, "real", 24, 0, "Summary field" ) + fieldList.append(field) + fieldList2 = vector.testForUniqueness(fieldList1, fieldList) + fieldList1.extend(fieldList) + seq = range(0, len(fieldList1)) + fieldList1 = dict(zip(seq, fieldList1)) + + sRs = provider1.crs() + progress.setPercentage(13) + fields = QgsFields() + for f in fieldList1.values(): + fields.append(f) + output = self.getOutputFromName(self.OUTPUT) + writer = output.getVectorWriter(fields, provider1.geometryType(), sRs) + + inFeat = QgsFeature() + outFeat = QgsFeature() + inFeatB = QgsFeature() + inGeom = QgsGeometry() + + progress.setPercentage(15) + start = 15.00 + add = 85.00 / provider1.featureCount() + + progress.setText('Creating spatial index...') + index = vector.spatialindex(layer2) + progress.setText('Processing spatial join...') + fit1 = provider1.getFeatures() + while fit1.nextFeature(inFeat): + inGeom = inFeat.geometry() + atMap1 = inFeat.attributes() + outFeat.setGeometry(inGeom) + none = True + joinList = [] + if inGeom.type() == QGis.Point: + joinList = index.intersects( inGeom.buffer(10,2).boundingBox() ) + if len(joinList) > 0: check = 0 + else: check = 1 + else: + joinList = index.intersects( inGeom.boundingBox() ) + if len(joinList) > 0: check = 0 + else: check = 1 + if check == 0: + count = 0 + for i in joinList: + provider2.getFeatures( QgsFeatureRequest().setFilterFid( int(i) ) ).nextFeature( inFeatB ) + if inGeom.intersects(inFeatB.geometry()): + count = count + 1 + none = False + atMap2 = inFeatB.attributes() + if not summary: + atMap = atMap1 + atMap2 = atMap2 + atMap.extend(atMap2) + atMap = dict(zip(seq, atMap)) + break + else: + for j in numFields.keys(): + numFields[j].append(atMap2[j]) + if summary and not none: + atMap = atMap1 + for j in numFields.keys(): + for k in sumList: + if k == "SUM": atMap.append(sum(numFields[j])) + elif k == "MEAN": atMap.append(sum(numFields[j]) / count) + elif k == "MIN": atMap.append(min(numFields[j])) + elif k == "MED": atMap.append(myself(numFields[j])) + else: atMap.append(max(numFields[j])) + numFields[j] = [] + atMap.append(count) + atMap = dict(zip(seq, atMap)) + if none: + outFeat.setAttributes(atMap1) + else: + outFeat.setAttributes(atMap.values()) + if keep: # keep all records + writer.addFeature(outFeat) + else: # keep only matching records + if not none: + writer.addFeature(outFeat) + start = start + add + progress.setPercentage(start) + + del writer diff --git a/python/plugins/processing/tools/vector.py b/python/plugins/processing/tools/vector.py index 4eefb7bfad7e..c83dde99a977 100644 --- a/python/plugins/processing/tools/vector.py +++ b/python/plugins/processing/tools/vector.py @@ -141,6 +141,20 @@ def values(layer, *attributes): ret[attr] = values return ret +def testForUniqueness( fieldList1, fieldList2 ): + '''Returns a modified version of fieldList2, removing naming + collisions with fieldList1.''' + changed = True + while changed: + changed = False + for i in range(0,len(fieldList1)): + for j in range(0,len(fieldList2)): + if fieldList1[i].name() == fieldList2[j].name(): + field = fieldList2[j] + name = createUniqueFieldName( field.name(), fieldList1 ) + fieldList2[j] = QgsField(name, field.type(), len=field.length(), prec=field.precision(), comment=field.comment()) + changed = True + return fieldList2 def spatialindex(layer): """Creates a spatial index for the passed vector layer. From 4e57af6931674b1cf3e0636cbe5f233200060aad Mon Sep 17 00:00:00 2001 From: Joshua Arnott Date: Tue, 22 Oct 2013 23:49:08 +0100 Subject: [PATCH 2/3] 'Join by location' can now use the geometry from the joined layer as output. --- .../processing/algs/ftools/SpatialJoin.py | 56 +++++++++++++++++-- 1 file changed, 50 insertions(+), 6 deletions(-) diff --git a/python/plugins/processing/algs/ftools/SpatialJoin.py b/python/plugins/processing/algs/ftools/SpatialJoin.py index 661ae11ac9c9..d37e0526301e 100644 --- a/python/plugins/processing/algs/ftools/SpatialJoin.py +++ b/python/plugins/processing/algs/ftools/SpatialJoin.py @@ -53,6 +53,16 @@ def myself(L): medianVal = L[ (nVal + 1) / 2 - 1] return medianVal +# translate single part to multi part geometry type +multipart = { + QGis.WKBPoint: QGis.WKBMultiPoint, + QGis.WKBLineString: QGis.WKBMultiLineString, + QGis.WKBPolygon: QGis.WKBMultiPolygon, + QGis.WKBPoint25D: QGis.WKBMultiPoint25D, + QGis.WKBLineString25D: QGis.WKBMultiLineString25D, + QGis.WKBPolygon25D: QGis.WKBMultiPolygon25D, +} + class SpatialJoin(GeoAlgorithm): ''' Join attributes by location @@ -64,6 +74,7 @@ class SpatialJoin(GeoAlgorithm): INPUT2 = "INPUT2" SUMMARY = "SUMMARY" STATS = "STATS" + GEOMETRY = "GEOMETRY" KEEP = "KEEP" OUTPUT = "OUTPUT" @@ -71,6 +82,11 @@ class SpatialJoin(GeoAlgorithm): 'Take attributes of the first located feature', 'Take summary of intersecting features' ] + + GEOMETRYS = [ + 'Use geometry from target layer', + 'Use geometry from joined layer (multipart if summary)' + ] KEEPS = [ 'Only keep matching records', @@ -88,7 +104,8 @@ def defineCharacteristics(self): self.addParameter(ParameterVector(SpatialJoin.INPUT1, "Target vector layer", [ParameterVector.VECTOR_TYPE_ANY])) self.addParameter(ParameterVector(SpatialJoin.INPUT2, "Join vector layer", [ParameterVector.VECTOR_TYPE_ANY])) self.addParameter(ParameterSelection(self.SUMMARY, "Attribute summary", self.SUMMARYS, 0)) - self.addParameter(ParameterString(self.STATS, "Statistics for summary (comma separated)", "sum,mean,min,max,med")) + self.addParameter(ParameterString(self.STATS, "Statistics for summary (comma separated)", "sum,mean,min,max,median")) + self.addParameter(ParameterSelection(self.GEOMETRY, "Output geometry", self.GEOMETRYS, 0)) self.addParameter(ParameterSelection(self.KEEP, "Output table", self.KEEPS, 0)) self.addOutput(OutputVector(SpatialJoin.OUTPUT, "Output layer")) @@ -97,6 +114,7 @@ def processAlgorithm(self, progress): summary = self.getParameterValue(self.SUMMARY) == 1 sumList = self.getParameterValue(self.STATS).upper().replace(' ','').split(',') + use_geom = self.getParameterValue(self.GEOMETRY) keep = self.getParameterValue(self.KEEP) == 1 input1 = self.getParameterValue(self.INPUT1) @@ -130,13 +148,25 @@ def processAlgorithm(self, progress): seq = range(0, len(fieldList1)) fieldList1 = dict(zip(seq, fieldList1)) - sRs = provider1.crs() progress.setPercentage(13) fields = QgsFields() for f in fieldList1.values(): fields.append(f) output = self.getOutputFromName(self.OUTPUT) - writer = output.getVectorWriter(fields, provider1.geometryType(), sRs) + + if use_geom == 0: + # from target layer + crs = provider1.crs() + geometry_type = provider1.geometryType() + else: + # from joined layer + crs = provider2.crs() + if summary: + geometry_type = multipart[provider2.geometryType()] + else: + geometry_type = provider2.geometryType() + + writer = output.getVectorWriter(fields, geometry_type, crs) inFeat = QgsFeature() outFeat = QgsFeature() @@ -154,7 +184,8 @@ def processAlgorithm(self, progress): while fit1.nextFeature(inFeat): inGeom = inFeat.geometry() atMap1 = inFeat.attributes() - outFeat.setGeometry(inGeom) + if use_geom == 0: + outFeat.setGeometry(inGeom) none = True joinList = [] if inGeom.type() == QGis.Point: @@ -171,17 +202,30 @@ def processAlgorithm(self, progress): provider2.getFeatures( QgsFeatureRequest().setFilterFid( int(i) ) ).nextFeature( inFeatB ) if inGeom.intersects(inFeatB.geometry()): count = count + 1 - none = False atMap2 = inFeatB.attributes() if not summary: + # first located feature atMap = atMap1 atMap2 = atMap2 atMap.extend(atMap2) atMap = dict(zip(seq, atMap)) + if use_geom == 1: + outFeat.setGeometry(inFeatB.geometry()) + none = False break else: for j in numFields.keys(): numFields[j].append(atMap2[j]) + if none: + # first located feature in summary + outFeat.setGeometry(inFeatB.geometry()) + else: + # combine (union) with existing geometry + existing = outFeat.geometry() + additional = inFeatB.geometry() + union = existing.combine(additional) + outFeat.setGeometry(union) + none = False if summary and not none: atMap = atMap1 for j in numFields.keys(): @@ -189,7 +233,7 @@ def processAlgorithm(self, progress): if k == "SUM": atMap.append(sum(numFields[j])) elif k == "MEAN": atMap.append(sum(numFields[j]) / count) elif k == "MIN": atMap.append(min(numFields[j])) - elif k == "MED": atMap.append(myself(numFields[j])) + elif k == "MEDIAN": atMap.append(myself(numFields[j])) else: atMap.append(max(numFields[j])) numFields[j] = [] atMap.append(count) From 90fa359fb76c4ac2c602d9f9f4d357b3d95cb890 Mon Sep 17 00:00:00 2001 From: Joshua Arnott Date: Sat, 26 Oct 2013 18:19:52 +0100 Subject: [PATCH 3/3] 'Join by location' now creates multipart features in the same way as SinglePartsToMultiparts. --- .../processing/algs/ftools/SpatialJoin.py | 77 ++++++++++++++----- 1 file changed, 56 insertions(+), 21 deletions(-) diff --git a/python/plugins/processing/algs/ftools/SpatialJoin.py b/python/plugins/processing/algs/ftools/SpatialJoin.py index d37e0526301e..414a804176d3 100644 --- a/python/plugins/processing/algs/ftools/SpatialJoin.py +++ b/python/plugins/processing/algs/ftools/SpatialJoin.py @@ -53,19 +53,9 @@ def myself(L): medianVal = L[ (nVal + 1) / 2 - 1] return medianVal -# translate single part to multi part geometry type -multipart = { - QGis.WKBPoint: QGis.WKBMultiPoint, - QGis.WKBLineString: QGis.WKBMultiLineString, - QGis.WKBPolygon: QGis.WKBMultiPolygon, - QGis.WKBPoint25D: QGis.WKBMultiPoint25D, - QGis.WKBLineString25D: QGis.WKBMultiLineString25D, - QGis.WKBPolygon25D: QGis.WKBMultiPolygon25D, -} - class SpatialJoin(GeoAlgorithm): ''' - Join attributes by location + Join by location Port of the spatial join algorithm from fTools to the Processing Toolbox. ''' @@ -99,7 +89,7 @@ class SpatialJoin(GeoAlgorithm): #=========================================================================== def defineCharacteristics(self): - self.name = "Join attributes by location" + self.name = "Join by location" self.group = "Vector general tools" self.addParameter(ParameterVector(SpatialJoin.INPUT1, "Target vector layer", [ParameterVector.VECTOR_TYPE_ANY])) self.addParameter(ParameterVector(SpatialJoin.INPUT2, "Join vector layer", [ParameterVector.VECTOR_TYPE_ANY])) @@ -162,7 +152,7 @@ def processAlgorithm(self, progress): # from joined layer crs = provider2.crs() if summary: - geometry_type = multipart[provider2.geometryType()] + geometry_type = self.singleToMultiGeom(provider2.geometryType()) else: geometry_type = provider2.geometryType() @@ -198,6 +188,7 @@ def processAlgorithm(self, progress): else: check = 1 if check == 0: count = 0 + multi_feature = [] for i in joinList: provider2.getFeatures( QgsFeatureRequest().setFilterFid( int(i) ) ).nextFeature( inFeatB ) if inGeom.intersects(inFeatB.geometry()): @@ -216,15 +207,12 @@ def processAlgorithm(self, progress): else: for j in numFields.keys(): numFields[j].append(atMap2[j]) - if none: - # first located feature in summary - outFeat.setGeometry(inFeatB.geometry()) + if use_geom == 0: + if none: + outFeat.setGeometry(inGeom) else: - # combine (union) with existing geometry - existing = outFeat.geometry() - additional = inFeatB.geometry() - union = existing.combine(additional) - outFeat.setGeometry(union) + feature_list = self.extractAsMulti(inFeatB.geometry()) + multi_feature.extend(feature_list) none = False if summary and not none: atMap = atMap1 @@ -238,6 +226,9 @@ def processAlgorithm(self, progress): numFields[j] = [] atMap.append(count) atMap = dict(zip(seq, atMap)) + if use_geom == 1: + outGeom = QgsGeometry(self.convertGeometry(multi_feature, geometry_type)) + outFeat.setGeometry(outGeom) if none: outFeat.setAttributes(atMap1) else: @@ -251,3 +242,47 @@ def processAlgorithm(self, progress): progress.setPercentage(start) del writer + + def singleToMultiGeom(self, wkbType): + try: + if wkbType in (QGis.WKBPoint, QGis.WKBMultiPoint, + QGis.WKBPoint25D, QGis.WKBMultiPoint25D): + return QGis.WKBMultiPoint + elif wkbType in (QGis.WKBLineString, QGis.WKBMultiLineString, + QGis.WKBMultiLineString25D, + QGis.WKBLineString25D): + + return QGis.WKBMultiLineString + elif wkbType in (QGis.WKBPolygon, QGis.WKBMultiPolygon, + QGis.WKBMultiPolygon25D, QGis.WKBPolygon25D): + + return QGis.WKBMultiPolygon + else: + return QGis.WKBUnknown + except Exception, err: + print unicode(err) + + def extractAsMulti(self, geom): + if geom.type() == QGis.Point: + if geom.isMultipart(): + return geom.asMultiPoint() + else: + return [geom.asPoint()] + elif geom.type() == QGis.Line: + if geom.isMultipart(): + return geom.asMultiPolyline() + else: + return [geom.asPolyline()] + else: + if geom.isMultipart(): + return geom.asMultiPolygon() + else: + return [geom.asPolygon()] + + def convertGeometry(self, geom_list, vType): + if vType == QGis.Point: + return QgsGeometry().fromMultiPoint(geom_list) + elif vType == QGis.Line: + return QgsGeometry().fromMultiPolyline(geom_list) + else: + return QgsGeometry().fromMultiPolygon(geom_list)