mirror of
https://github.com/qgis/QGIS.git
synced 2025-02-27 00:33:48 -05:00
899c57a91f
filter_null method previously filtered the NULL values using NULL type in qgis.core which is a QPyNullVariant.
f9e0093ce7
This was subsequently changed to 'is not None' comparison which is incorrect. The current version fails to filter out the null values.
277 lines
12 KiB
Python
277 lines
12 KiB
Python
# -*- coding: utf-8 -*-
|
|
#-----------------------------------------------------------
|
|
#
|
|
# fTools
|
|
# Copyright (C) 2008-2011 Carson Farmer
|
|
# EMAIL: carson.farmer (at) gmail.com
|
|
# WEB : http://www.ftools.ca/fTools.html
|
|
#
|
|
# A collection of data management and analysis tools for vector data
|
|
#
|
|
#-----------------------------------------------------------
|
|
#
|
|
# licensed under the terms of GNU GPL 2
|
|
#
|
|
# 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.
|
|
#
|
|
# This program is distributed in the hope that it will be useful,
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
# GNU General Public License for more details.
|
|
#
|
|
# You should have received a copy of the GNU General Public License along
|
|
# with this program; if not, write to the Free Software Foundation, Inc.,
|
|
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
|
#
|
|
#---------------------------------------------------------------------
|
|
|
|
from PyQt4.QtCore import QObject, SIGNAL, QVariant, QFile
|
|
from PyQt4.QtGui import QDialog, QDialogButtonBox, QMessageBox
|
|
from qgis.core import QGis, QgsVectorFileWriter, QgsVectorLayer, QgsMapLayerRegistry, QgsFields, QgsField, QgsFeature, QgsGeometry, NULL
|
|
|
|
import ftools_utils
|
|
from ui_frmSpatialJoin import Ui_Dialog
|
|
|
|
|
|
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
|
|
|
|
|
|
def filter_null(vals):
|
|
"""Takes an iterator of values and returns a new iterator returning the same values but skipping any NULL values"""
|
|
return (v for v in vals if v != NULL)
|
|
|
|
|
|
class Dialog(QDialog, Ui_Dialog):
|
|
|
|
def __init__(self, iface):
|
|
QDialog.__init__(self, iface.mainWindow())
|
|
self.iface = iface
|
|
# Set up the user interface from Designer.
|
|
self.setupUi(self)
|
|
QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
|
|
self.setWindowTitle(self.tr("Join attributes by location"))
|
|
self.buttonOk = self.buttonBox_2.button(QDialogButtonBox.Ok)
|
|
# populate layer list
|
|
self.progressBar.setValue(0)
|
|
layers = ftools_utils.getLayerNames([QGis.Point, QGis.Line, QGis.Polygon])
|
|
self.inShape.addItems(layers)
|
|
self.joinShape.addItems(layers)
|
|
|
|
def accept(self):
|
|
self.buttonOk.setEnabled(False)
|
|
if self.inShape.currentText() == "":
|
|
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify target vector layer"))
|
|
elif self.outShape.text() == "":
|
|
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify output shapefile"))
|
|
elif self.joinShape.currentText() == "":
|
|
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify join vector layer"))
|
|
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() or self.chkMedian.isChecked()):
|
|
QMessageBox.information(self, self.tr("Spatial Join"), self.tr("Please specify at least one summary statistic"))
|
|
else:
|
|
inName = self.inShape.currentText()
|
|
joinName = self.joinShape.currentText()
|
|
outPath = self.outShape.text()
|
|
if self.rdoSummary.isChecked():
|
|
summary = True
|
|
sumList = []
|
|
if self.chkSum.isChecked():
|
|
sumList.append("SUM")
|
|
if self.chkMean.isChecked():
|
|
sumList.append("MEAN")
|
|
if self.chkMin.isChecked():
|
|
sumList.append("MIN")
|
|
if self.chkMax.isChecked():
|
|
sumList.append("MAX")
|
|
if self.chkMedian.isChecked():
|
|
sumList.append("MED")
|
|
else:
|
|
summary = False
|
|
sumList = ["all"]
|
|
if self.rdoKeep.isChecked():
|
|
keep = True
|
|
else:
|
|
keep = False
|
|
outName = ftools_utils.getShapefileName(outPath)
|
|
res = self.compute(inName, joinName, outPath, summary, sumList, keep, self.progressBar)
|
|
self.outShape.clear()
|
|
if res:
|
|
addToTOC = QMessageBox.question(
|
|
self, self.tr("Spatial Join"),
|
|
self.tr("Created output shapefile:\n%s\n\nWould you like to add the new layer to the TOC?") % (unicode(outPath)),
|
|
QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
|
|
if addToTOC == QMessageBox.Yes:
|
|
self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
|
|
QgsMapLayerRegistry.instance().addMapLayers([self.vlayer])
|
|
self.progressBar.setValue(0)
|
|
self.buttonOk.setEnabled(True)
|
|
|
|
def outFile(self):
|
|
self.outShape.clear()
|
|
(self.shapefileName, self.encoding) = ftools_utils.saveDialog(self)
|
|
if self.shapefileName is None or self.encoding is None:
|
|
return
|
|
self.outShape.setText(self.shapefileName)
|
|
|
|
def compute(self, inName, joinName, outName, summary, sumList, keep, progressBar):
|
|
layer1 = ftools_utils.getVectorLayerByName(inName)
|
|
provider1 = layer1.dataProvider()
|
|
fieldList1 = ftools_utils.getFieldList(layer1)
|
|
|
|
layer2 = ftools_utils.getVectorLayerByName(joinName)
|
|
provider2 = layer2.dataProvider()
|
|
|
|
fieldList2 = ftools_utils.getFieldList(layer2)
|
|
fieldList = QgsFields()
|
|
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."))
|
|
if not summary:
|
|
fieldList2 = ftools_utils.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, self.tr("Summary field"))
|
|
fieldList.append(field)
|
|
field = QgsField("COUNT", QVariant.Double, "real", 24, 16, self.tr("Summary field"))
|
|
fieldList.append(field)
|
|
fieldList2 = ftools_utils.testForUniqueness(fieldList1, fieldList)
|
|
fieldList1.extend(fieldList)
|
|
seq = range(0, len(fieldList1))
|
|
fieldList1 = dict(zip(seq, fieldList1))
|
|
|
|
sRs = provider1.crs()
|
|
progressBar.setValue(13)
|
|
check = QFile(self.shapefileName)
|
|
if check.exists():
|
|
if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
|
|
QMessageBox.warning(
|
|
self, self.tr('Error deleting shapefile'),
|
|
self.tr("Can't delete existing shapefile\n%s") % (self.shapefileName))
|
|
return False
|
|
fields = QgsFields()
|
|
for f in fieldList1.values():
|
|
fields.append(f)
|
|
writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fields, provider1.geometryType(), sRs)
|
|
#writer = QgsVectorFileWriter(outName, "UTF-8", fieldList1, provider1.geometryType(), sRs)
|
|
inFeat = QgsFeature()
|
|
outFeat = QgsFeature()
|
|
inFeatB = QgsFeature()
|
|
inGeom = QgsGeometry()
|
|
progressBar.setValue(15)
|
|
start = 15.00
|
|
add = 85.00 / provider1.featureCount()
|
|
|
|
index = ftools_utils.createIndex(provider2)
|
|
|
|
# cache all features from provider2 to avoid huge number of feature requests in the inner loop
|
|
mapP2 = {}
|
|
for f in provider2.getFeatures():
|
|
mapP2[f.id()] = QgsFeature(f)
|
|
|
|
fit1 = provider1.getFeatures()
|
|
while fit1.nextFeature(inFeat):
|
|
inGeom = inFeat.geometry()
|
|
atMap1 = inFeat.attributes()
|
|
outFeat.setGeometry(inGeom)
|
|
none = True
|
|
joinList = []
|
|
if inGeom.type() == QGis.Point:
|
|
#(check, joinList) = layer2.featuresInRectangle(inGeom.buffer(10,2).boundingBox(), True, True)
|
|
#layer2.select(inGeom.buffer(10,2).boundingBox(), False)
|
|
#joinList = layer2.selectedFeatures()
|
|
joinList = index.intersects(inGeom.buffer(10, 2).boundingBox())
|
|
if len(joinList) > 0:
|
|
check = 0
|
|
else:
|
|
check = 1
|
|
else:
|
|
#(check, joinList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
|
|
#layer2.select(inGeom.boundingBox(), False)
|
|
#joinList = layer2.selectedFeatures()
|
|
joinList = index.intersects(inGeom.boundingBox())
|
|
if len(joinList) > 0:
|
|
check = 0
|
|
else:
|
|
check = 1
|
|
if check == 0:
|
|
count = 0
|
|
for i in joinList:
|
|
inFeatB = mapP2[i] # cached feature from provider2
|
|
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(filter_null(numFields[j])))
|
|
elif k == "MEAN":
|
|
try:
|
|
nn_count = sum(1 for _ in filter_null(numFields[j]))
|
|
atMap.append(sum(filter_null(numFields[j])) / nn_count)
|
|
except ZeroDivisionError:
|
|
atMap.append(NULL)
|
|
elif k == "MIN":
|
|
try:
|
|
atMap.append(min(filter_null(numFields[j])))
|
|
except ValueError:
|
|
atMap.append(NULL)
|
|
elif k == "MED":
|
|
atMap.append(myself(numFields[j]))
|
|
else:
|
|
try:
|
|
atMap.append(max(filter_null(numFields[j])))
|
|
except ValueError:
|
|
atMap.append(NULL)
|
|
|
|
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
|
|
progressBar.setValue(start)
|
|
del writer
|
|
return True
|