[QGIS Commit] r13035 - trunk/qgis/python/plugins/fTools/tools

svn_qgis at osgeo.org svn_qgis at osgeo.org
Tue Mar 9 19:16:05 EST 2010


Author: cfarmer
Date: 2010-03-09 19:16:03 -0500 (Tue, 09 Mar 2010)
New Revision: 13035

Modified:
   trunk/qgis/python/plugins/fTools/tools/doIntersectLines.py
   trunk/qgis/python/plugins/fTools/tools/doPointsInPolygon.py
Log:
Uses spatial index to select intersecting features (used selections before). Tools should both run faster and with fewer map canvas glitches.

Modified: trunk/qgis/python/plugins/fTools/tools/doIntersectLines.py
===================================================================
--- trunk/qgis/python/plugins/fTools/tools/doIntersectLines.py	2010-03-09 00:52:39 UTC (rev 13034)
+++ trunk/qgis/python/plugins/fTools/tools/doIntersectLines.py	2010-03-10 00:16:03 UTC (rev 13035)
@@ -39,157 +39,133 @@
 
 class Dialog(QDialog, Ui_Dialog):
 
-	def __init__(self, iface):
-		QDialog.__init__(self)
-		self.iface = iface
-		# Set up the user interface from Designer.
-		self.setupUi(self)
-		QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
-		QObject.connect(self.inLine1, SIGNAL("currentIndexChanged(QString)"), self.update1)
-		QObject.connect(self.inLine2, SIGNAL("currentIndexChanged(QString)"), self.update2)
-		self.setWindowTitle( self.tr("Line intersections") )
-		# populate layer list
-		self.progressBar.setValue(0)
-		mapCanvas = self.iface.mapCanvas()
-		for i in range(mapCanvas.layerCount()):
-			layer = mapCanvas.layer(i)
-			if layer.type() == layer.VectorLayer:
-				if layer.geometryType() == QGis.Line:
-					self.inLine1.addItem(layer.name())
-					self.inLine2.addItem(layer.name())
+    def __init__(self, iface):
+        QDialog.__init__(self)
+        self.iface = iface
+        # Set up the user interface from Designer.
+        self.setupUi(self)
+        QObject.connect(self.toolOut, SIGNAL("clicked()"), self.outFile)
+        QObject.connect(self.inLine1, SIGNAL("currentIndexChanged(QString)"), self.update1)
+        QObject.connect(self.inLine2, SIGNAL("currentIndexChanged(QString)"), self.update2)
+        self.setWindowTitle( self.tr("Line intersections") )
+        # populate layer list
+        self.progressBar.setValue(0)
+        mapCanvas = self.iface.mapCanvas()
+        layers = ftools_utils.getLayerNames([QGis.Line])
+        self.inLine1.addItems(layers)
+        self.inLine2.addItems(layers)
 
-	def update1(self, inputLayer):
-		self.inField1.clear()
-		changedLayer = self.getVectorLayerByName(unicode(inputLayer))
-		changedField = self.getFieldList(changedLayer)
-		for i in changedField:
-			self.inField1.addItem(unicode(changedField[i].name()))
+    def update1(self, inputLayer):
+        self.inField1.clear()
+        changedLayer = ftools_utils.getVectorLayerByName(unicode(inputLayer))
+        changedField = ftools_utils.getFieldList(changedLayer)
+        for i in changedField:
+            self.inField1.addItem(unicode(changedField[i].name()))
 
-	def update2(self, inputLayer):
-		self.inField2.clear()
-		changedLayer = self.getVectorLayerByName(unicode(inputLayer))
-		changedField = self.getFieldList(changedLayer)
-		for i in changedField:
-			self.inField2.addItem(unicode(changedField[i].name()))
+    def update2(self, inputLayer):
+        self.inField2.clear()
+        changedLayer = ftools_utils.getVectorLayerByName(unicode(inputLayer))
+        changedField = ftools_utils.getFieldList(changedLayer)
+        for i in changedField:
+            self.inField2.addItem(unicode(changedField[i].name()))
 
-	def accept(self):
-		if self.inLine1.currentText() == "":
-			QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input line layer") )
-		elif self.outShape.text() == "":
-			QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify output shapefile") )
-		elif self.inLine2.currentText() == "":
-			QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify line intersect layer") )
-		elif self.inField1.currentText() == "":
-			QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input unique ID field") )
-		elif self.inField2.currentText() == "":
-			QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify intersect unique ID field") )
-		else:
-			line1 = self.inLine1.currentText()
-			line2 = self.inLine2.currentText()
-			field1 = self.inField1.currentText()
-			field2 = self.inField2.currentText()
-			outPath = self.outShape.text()
-			if outPath.contains("\\"):
-				outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
-			else:
-				outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
-			if outName.endsWith(".shp"):
-				outName = outName.left(outName.length() - 4)
-			self.outShape.clear()
-			self.compute(line1, line2, field1, field2, outPath, self.progressBar)
-			self.progressBar.setValue(100)
-			addToTOC = QMessageBox.question(self, self.tr("Generate Centroids"), self.tr("Created output point shapefile:\n%1\n\nWould you like to add the new layer to the TOC?").arg( outPath ), QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
-			if addToTOC == QMessageBox.Yes:
-				self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
-				QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
-		self.progressBar.setValue(0)
+    def accept(self):
+        if self.inLine1.currentText() == "":
+            QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input line layer") )
+        elif self.outShape.text() == "":
+            QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify output shapefile") )
+        elif self.inLine2.currentText() == "":
+            QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify line intersect layer") )
+        elif self.inField1.currentText() == "":
+            QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify input unique ID field") )
+        elif self.inField2.currentText() == "":
+            QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Please specify intersect unique ID field") )
+        else:
+            line1 = self.inLine1.currentText()
+            line2 = self.inLine2.currentText()
+            field1 = self.inField1.currentText()
+            field2 = self.inField2.currentText()
+            outPath = self.outShape.text()
+            self.outShape.clear()
+            self.compute(line1, line2, field1, field2, outPath, self.progressBar)
+            self.progressBar.setValue(100)
+            addToTOC = QMessageBox.question(self, self.tr("Generate Centroids"), self.tr("Created output point shapefile:\n%1\n\nWould you like to add the new layer to the TOC?").arg( outPath ), QMessageBox.Yes, QMessageBox.No, QMessageBox.NoButton)
+            if addToTOC == QMessageBox.Yes:
+                if not ftools_utils.addShapeToCanvas( unicode( outPath ) ):
+                    QMessageBox.warning( self, self.tr("Geoprocessing"), self.tr( "Error loading output shapefile:\n%1" ).arg( unicode( outPath ) ))
+        self.progressBar.setValue(0)
 
-	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( QString( self.shapefileName ) )
+    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( QString( self.shapefileName ) )
 
-	def compute(self, line1, line2, field1, field2, outPath, progressBar):
-		layer1 = self.getVectorLayerByName(line1)
-		layer2 = self.getVectorLayerByName(line2)
-		provider1 = layer1.dataProvider()
-		provider2 = layer2.dataProvider()
-		allAttrs = provider1.attributeIndexes()
-		provider1.select(allAttrs)
-		allAttrs = provider2.attributeIndexes()
-		provider2.select(allAttrs)
-		fieldList = self.getFieldList(layer1)
-		index1 = provider1.fieldNameIndex(field1)
-		field1 = fieldList[index1]
-		field1.setName(unicode(field1.name()) + "_1")
-		fieldList = self.getFieldList(layer2)
-		index2 = provider2.fieldNameIndex(field2)
-		field2 = fieldList[index2]
-		field2.setName(unicode(field2.name()) + "_2")
-		fieldList = {0:field1, 1:field2}
-		sRs = provider1.crs()
-		check = QFile(self.shapefileName)
-		if check.exists():
-			if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
-				return
-		writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, QGis.WKBPoint, sRs)
-		#writer = QgsVectorFileWriter(outPath, "UTF-8", fieldList, QGis.WKBPoint, sRs)
-		inFeat = QgsFeature()
-		outFeat = QgsFeature()
-		inGeom = QgsGeometry()
-		tempGeom = QgsGeometry()
-		start = 15.00
-		add = 85.00 / layer1.featureCount()
-		while provider1.nextFeature(inFeat):
-			inGeom = inFeat.geometry()
-			lineList = []
-			#(check, lineList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
-			# Below is a work-around for featuresInRectangle
-			# Which appears to have been removed for QGIS version 1.0
-			layer2.select(inGeom.boundingBox(), False)
-			lineList = layer2.selectedFeatures()
-			if len(lineList) > 0: check = 0
-			else: check = 1
-			if check == 0:
-				for i in lineList:
-					tempGeom = i.geometry()
-					tempList = []
-					atMap1 = inFeat.attributeMap()
-					atMap2 = i.attributeMap()
-					if inGeom.intersects(tempGeom):
-						tempGeom = inGeom.intersection(tempGeom)
-						if tempGeom.type() == QGis.Point:
-							if tempGeom.isMultipart():
-								tempList = tempGeom.asMultiPoint()
-							else:
-								tempList.append(tempGeom.asPoint())
-							for j in tempList:
-								outFeat.setGeometry(tempGeom.fromPoint(j))
-								outFeat.addAttribute(0, atMap1[index1])
-								outFeat.addAttribute(1, atMap2[index2])
-								writer.addFeature(outFeat)
-			start = start + add
-			progressBar.setValue(start)
-		del writer
-
-	def getVectorLayerByName(self, myName):
-		mc = self.iface.mapCanvas()
-		nLayers = mc.layerCount()
-		for l in range(nLayers):
-			layer = mc.layer(l)
-			if layer.name() == unicode(myName):
-				vlayer = QgsVectorLayer(unicode(layer.source()),  unicode(myName),  unicode(layer.dataProvider().name()))
-				if vlayer.isValid():
-					return vlayer
-				else:
-					QMessageBox.information(self, self.tr("Locate Line Intersections"), self.tr("Vector layer is not valid") )
-
-	def getFieldList(self, vlayer):
-		fProvider = vlayer.dataProvider()
-		feat = QgsFeature()
-		allAttrs = fProvider.attributeIndexes()
-		fProvider.select(allAttrs)
-		myFields = fProvider.fields()
-		return myFields
+    def compute(self, line1, line2, field1, field2, outPath, progressBar):
+        layer1 = ftools_utils.getVectorLayerByName(line1)
+        layer2 = ftools_utils.getVectorLayerByName(line2)
+        provider1 = layer1.dataProvider()
+        provider2 = layer2.dataProvider()
+        allAttrs = provider1.attributeIndexes()
+        provider1.select(allAttrs)
+        allAttrs = provider2.attributeIndexes()
+        provider2.select(allAttrs)
+        fieldList = ftools_utils.getFieldList(layer1)
+        index1 = provider1.fieldNameIndex(field1)
+        field1 = fieldList[index1]
+        field1.setName(unicode(field1.name()) + "_1")
+        fieldList = ftools_utils.getFieldList(layer2)
+        index2 = provider2.fieldNameIndex(field2)
+        field2 = fieldList[index2]
+        field2.setName(unicode(field2.name()) + "_2")
+        fieldList = {0:field1, 1:field2}
+        sRs = provider1.crs()
+        check = QFile(self.shapefileName)
+        if check.exists():
+            if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
+                return
+        writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, QGis.WKBPoint, sRs)
+        #writer = QgsVectorFileWriter(outPath, "UTF-8", fieldList, QGis.WKBPoint, sRs)
+        inFeat = QgsFeature()
+        inFeatB = QgsFeature()
+        outFeat = QgsFeature()
+        inGeom = QgsGeometry()
+        tempGeom = QgsGeometry()
+        start = 15.00
+        add = 85.00 / layer1.featureCount()
+        index = ftools_utils.createIndex( provider2 )
+        while provider1.nextFeature(inFeat):
+            inGeom = inFeat.geometry()
+            lineList = []
+            #(check, lineList) = layer2.featuresInRectangle(inGeom.boundingBox(), True, True)
+            # Below is a work-around for featuresInRectangle
+            # Which appears to have been removed for QGIS version 1.0
+            #layer2.select(inGeom.boundingBox(), False)
+            #lineList = layer2.selectedFeatures()
+            lineList = index.intersects( inGeom.boundingBox() )
+            if len(lineList) > 0: check = 0
+            else: check = 1
+            if check == 0:
+                for i in lineList:
+                    provider2.featureAtId( int( i ), inFeatB , True, allAttrs )
+                    tmpGeom = QgsGeometry( inFeatB.geometry() )
+                    #tempGeom = i.geometry()
+                    tempList = []
+                    atMap1 = inFeat.attributeMap()
+                    atMap2 = inFeatB.attributeMap()
+                    if inGeom.intersects(tmpGeom):
+                        tempGeom = inGeom.intersection(tmpGeom)
+                        if tempGeom.type() == QGis.Point:
+                            if tempGeom.isMultipart():
+                                tempList = tempGeom.asMultiPoint()
+                            else:
+                                tempList.append(tempGeom.asPoint())
+                            for j in tempList:
+                                outFeat.setGeometry(tempGeom.fromPoint(j))
+                                outFeat.addAttribute(0, atMap1[index1])
+                                outFeat.addAttribute(1, atMap2[index2])
+                                writer.addFeature(outFeat)
+            start = start + add
+            progressBar.setValue(start)
+        del writer

Modified: trunk/qgis/python/plugins/fTools/tools/doPointsInPolygon.py
===================================================================
--- trunk/qgis/python/plugins/fTools/tools/doPointsInPolygon.py	2010-03-09 00:52:39 UTC (rev 13034)
+++ trunk/qgis/python/plugins/fTools/tools/doPointsInPolygon.py	2010-03-10 00:16:03 UTC (rev 13035)
@@ -1,4 +1,5 @@
-#-----------------------------------------------------------
+# -*- coding: utf-8 -*-
+#-----------------------------------------------------------
 # 
 # Points in Polygon
 #
@@ -38,127 +39,109 @@
 
 class Dialog(QDialog, Ui_Dialog):
 
-	def __init__(self, iface):
-		QDialog.__init__(self)
-		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("Count Points in Polygon"))
-		# populate layer list
-		self.progressBar.setValue(0)
-		mapCanvas = self.iface.mapCanvas()
-		for i in range(mapCanvas.layerCount()):
-			layer = mapCanvas.layer(i)
-			if layer.type() == layer.VectorLayer:
-				if layer.geometryType() == QGis.Polygon:
-					self.inPolygon.addItem(layer.name())
-				elif layer.geometryType() == QGis.Point:
-					self.inPoint.addItem(layer.name())
-		
-	def accept(self):
-		if self.inPolygon.currentText() == "":
-			QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify input polygon vector layer"))
-		elif self.outShape.text() == "":
-			QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify output shapefile"))
-		elif self.inPoint.currentText() == "":
-			QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify input point vector layer"))
-		elif self.lnField.text() == "":
-			QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify output count field"))
-		else:
-			inPoly = self.inPolygon.currentText()
-			inPts = self.inPoint.currentText()
-			inField = self.lnField.text()
-			outPath = self.outShape.text()
-			if outPath.contains("\\"):
-				outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
-			else:
-				outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
-			if outName.endsWith(".shp"):
-				outName = outName.left(outName.length() - 4)
-			self.compute(inPoly, inPts, inField, outPath, self.progressBar)
-			self.outShape.clear()
-			addToTOC = QMessageBox.question(self, self.tr("Count Points in Polygon"), 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)
-			if addToTOC == QMessageBox.Yes:
-				self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
-				QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
-		self.progressBar.setValue(0)
+    def __init__(self, iface):
+        QDialog.__init__(self)
+        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("Count Points in Polygon"))
+        # populate layer list
+        self.progressBar.setValue(0)
+        mapCanvas = self.iface.mapCanvas()
+        layers = ftools_utils.getLayerNames([QGis.Polygon])
+        self.inPolygon.addItems(layers)
+        layers = ftools_utils.getLayerNames([QGis.Point])
+        self.inPoint.addItems(layers)
+    
+    def accept(self):
+        if self.inPolygon.currentText() == "":
+            QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify input polygon vector layer"))
+        elif self.outShape.text() == "":
+            QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify output shapefile"))
+        elif self.inPoint.currentText() == "":
+            QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify input point vector layer"))
+        elif self.lnField.text() == "":
+            QMessageBox.information(self, self.tr("Count Points In Polygon"), self.tr("Please specify output count field"))
+        else:
+            inPoly = self.inPolygon.currentText()
+            inPts = self.inPoint.currentText()
+            inField = self.lnField.text()
+            outPath = self.outShape.text()
+            if outPath.contains("\\"):
+                outName = outPath.right((outPath.length() - outPath.lastIndexOf("\\")) - 1)
+            else:
+                outName = outPath.right((outPath.length() - outPath.lastIndexOf("/")) - 1)
+            if outName.endsWith(".shp"):
+                outName = outName.left(outName.length() - 4)
+            self.compute(inPoly, inPts, inField, outPath, self.progressBar)
+            self.outShape.clear()
+            addToTOC = QMessageBox.question(self, self.tr("Count Points in Polygon"), 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)
+            if addToTOC == QMessageBox.Yes:
+                self.vlayer = QgsVectorLayer(outPath, unicode(outName), "ogr")
+                QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)
+        self.progressBar.setValue(0)
 
-	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( QString( self.shapefileName ) )
+    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( QString( self.shapefileName ) )
 
-	def compute(self, inPoly, inPts, inField, outPath, progressBar):
-		polyLayer = self.getVectorLayerByName(inPoly)
-		pointLayer = self.getVectorLayerByName(inPts)
-		polyProvider = polyLayer.dataProvider()
-		pointProvider = pointLayer.dataProvider()
-		if polyProvider.crs() <> pointProvider.crs():
-			QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results."))
-		allAttrs = polyProvider.attributeIndexes()
-		polyProvider.select(allAttrs)
-		allAttrs = pointProvider.attributeIndexes()
-		pointProvider.select(allAttrs)
-		fieldList = self.getFieldList(polyLayer)
-		index = polyProvider.fieldNameIndex(unicode(inField))
-		if index == -1:
-			index = polyProvider.fieldCount()
-			field = QgsField(unicode(inField), QVariant.Int, "real", 24, 15, self.tr("point count field"))
-			fieldList[index] = field
-		sRs = polyProvider.crs()
-		check = QFile(self.shapefileName)
-		if check.exists():
-			if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
-				return
-		writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, polyProvider.geometryType(), sRs)
-		#writer = QgsVectorFileWriter(outPath, "UTF-8", fieldList, polyProvider.geometryType(), sRs)
-		inFeat = QgsFeature()
-		outFeat = QgsFeature()
-		inGeom = QgsGeometry()
-		start = 15.00
-		add = 85.00 / polyProvider.featureCount()
-		while polyProvider.nextFeature(inFeat):
-			inGeom = inFeat.geometry()
-			atMap = inFeat.attributeMap()
-			outFeat.setAttributeMap(atMap)
-			outFeat.setGeometry(inGeom)
-			pointList = []
-			count = 0
-			#(check, pointList) = pointLayer.featuresInRectangle(inGeom.boundingBox(), True, True)
-			pointLayer.select(inGeom.boundingBox(), False)
-			pointList = pointLayer.selectedFeatures()
-			if len(pointList) > 0: check = 0
-			else: check = 1
-			if check == 0:
-				for i in pointList:
-					if inGeom.contains(i.geometry().asPoint()):
-						count = count + 1
-			outFeat.setAttributeMap(atMap)
-			outFeat.addAttribute(index, QVariant(count))
-			writer.addFeature(outFeat)
-			start = start + 1
-			progressBar.setValue(start)
-		del writer
-				
-	def getVectorLayerByName(self, myName):
-		mc = self.iface.mapCanvas()
-		nLayers = mc.layerCount()
-		for l in range(nLayers):
-			layer = mc.layer(l)
-			if layer.name() == unicode(myName):
-				vlayer = QgsVectorLayer(unicode(layer.source()),  unicode(myName),  unicode(layer.dataProvider().name()))
-				if vlayer.isValid():
-					return vlayer
-				else:
-					QMessageBox.information(self, self.tr("Counts Points In Polygon"), self.tr("Vector layer is not valid"))
-
-	def getFieldList(self, vlayer):
-		fProvider = vlayer.dataProvider()
-		feat = QgsFeature()
-		allAttrs = fProvider.attributeIndexes()
-		fProvider.select(allAttrs)
-		myFields = fProvider.fields()
-		return myFields
+    def compute(self, inPoly, inPts, inField, outPath, progressBar):
+        polyLayer = ftools_utils.getVectorLayerByName(inPoly)
+        pointLayer = ftools_utils.getVectorLayerByName(inPts)
+        polyProvider = polyLayer.dataProvider()
+        pointProvider = pointLayer.dataProvider()
+        if polyProvider.crs() <> pointProvider.crs():
+            QMessageBox.warning(self, self.tr("CRS warning!"), self.tr("Warning: Input layers have non-matching CRS.\nThis may cause unexpected results."))
+        allAttrs = polyProvider.attributeIndexes()
+        polyProvider.select(allAttrs)
+        allAttrs = pointProvider.attributeIndexes()
+        pointProvider.select(allAttrs)
+        fieldList = ftools_utils.getFieldList(polyLayer)
+        index = polyProvider.fieldNameIndex(unicode(inField))
+        if index == -1:
+            index = polyProvider.fieldCount()
+            field = QgsField(unicode(inField), QVariant.Int, "real", 24, 15, self.tr("point count field"))
+            fieldList[index] = field
+        sRs = polyProvider.crs()
+        check = QFile(self.shapefileName)
+        if check.exists():
+            if not QgsVectorFileWriter.deleteShapeFile(self.shapefileName):
+                return
+        writer = QgsVectorFileWriter(self.shapefileName, self.encoding, fieldList, polyProvider.geometryType(), sRs)
+        #writer = QgsVectorFileWriter(outPath, "UTF-8", fieldList, polyProvider.geometryType(), sRs)
+        inFeat = QgsFeature()
+        inFeatB = QgsFeature()
+        outFeat = QgsFeature()
+        inGeom = QgsGeometry()
+        start = 15.00
+        add = 85.00 / polyProvider.featureCount()
+        spatialIndex = ftools_utils.createIndex( pointProvider )
+        while polyProvider.nextFeature(inFeat):
+            inGeom = inFeat.geometry()
+            atMap = inFeat.attributeMap()
+            outFeat.setAttributeMap(atMap)
+            outFeat.setGeometry(inGeom)
+            pointList = []
+            count = 0
+            #(check, pointList) = pointLayer.featuresInRectangle(inGeom.boundingBox(), True, True)
+            #pointLayer.select(inGeom.boundingBox(), False)
+            #pointList = pointLayer.selectedFeatures()
+            pointList = spatialIndex.intersects(inGeom.boundingBox())
+            if len(pointList) > 0: check = 0
+            else: check = 1
+            if check == 0:
+                for i in pointList:
+                    pointProvider.featureAtId( int( i ), inFeatB , True, allAttrs )
+                    tmpGeom = QgsGeometry( inFeatB.geometry() )
+                    if inGeom.contains(tmpGeom.asPoint()):
+                        count = count + 1
+            outFeat.setAttributeMap(atMap)
+            outFeat.addAttribute(index, QVariant(count))
+            writer.addFeature(outFeat)
+            start = start + 1
+            progressBar.setValue(start)
+        del writer



More information about the QGIS-commit mailing list