Skip to content

Commit

Permalink
review Clip tool
Browse files Browse the repository at this point in the history
  • Loading branch information
alexbruy committed Mar 21, 2013
1 parent 48e9f42 commit aa4f70c
Showing 1 changed file with 71 additions and 76 deletions.
147 changes: 71 additions & 76 deletions python/plugins/sextante/algs/ftools/Clip.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,113 +23,108 @@
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'

from sextante.core.GeoAlgorithm import GeoAlgorithm
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from sextante.parameters.ParameterVector import ParameterVector

from sextante.core.GeoAlgorithm import GeoAlgorithm
from sextante.core.QGisLayers import QGisLayers
from sextante.core.SextanteLog import SextanteLog

from sextante.parameters.ParameterVector import ParameterVector

from sextante.outputs.OutputVector import OutputVector

from sextante.algs.ftools import FToolsUtils as utils
from sextante.core.SextanteLog import SextanteLog

class Clip(GeoAlgorithm):

INPUT = "INPUT"
INPUT2 = "INPUT2"
OVERLAY = "OVERLAY"
OUTPUT = "OUTPUT"

#===========================================================================
# def getIcon(self):
# return QtGui.QIcon(os.path.dirname(__file__) + "/icons/clip.png")
#===========================================================================

def defineCharacteristics(self):
self.name = "Clip"
self.group = "Vector overlay tools"
self.addParameter(ParameterVector(Clip.INPUT, "Input layer", ParameterVector.VECTOR_TYPE_ANY))
self.addParameter(ParameterVector(Clip.OVERLAY, "Clip layer", ParameterVector.VECTOR_TYPE_ANY))
self.addOutput(OutputVector(Clip.OUTPUT, "Clipped"))

def processAlgorithm(self, progress):
settings = QSettings()
systemEncoding = settings.value( "/UI/encoding", "System" ).toString()
output = self.getOutputValue(Clip.OUTPUT)
vlayerA = QGisLayers.getObjectFromUri(self.getParameterValue(Clip.INPUT))
vlayerB = QGisLayers.getObjectFromUri(self.getParameterValue(Clip.INPUT2))
layerA = QGisLayers.getObjectFromUri(self.getParameterValue(Clip.INPUT))
layerB = QGisLayers.getObjectFromUri(self.getParameterValue(Clip.OVERLAY))

GEOS_EXCEPT = True
FEATURE_EXCEPT = True
vproviderA = vlayerA.dataProvider()
#allAttrsA = vproviderA.attributeIndexes()
#vproviderA.select( allAttrsA )
vproviderB = vlayerB.dataProvider()
#allAttrsB = vproviderB.attributeIndexes()
#vproviderB.select( allAttrsB )

# check for crs compatibility
crsA = vproviderA.crs()
crsB = vproviderB.crs()
if not crsA.isValid() or not crsB.isValid():
SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Intersection. Invalid CRS. Results might be unexpected")
else:
if crsA != crsB:
SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Intersection. Non-matching CRSs. Results might be unexpected")
fields = vproviderA.fields()
writer = QgsVectorFileWriter( output, systemEncoding,fields, vproviderA.geometryType(), vproviderA.crs() )

writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(layerA.pendingFields(),
layerA.dataProvider().geometryType(), layerA.dataProvider().crs())

inFeatA = QgsFeature()
inFeatB = QgsFeature()
outFeat = QgsFeature()
index = utils.createSpatialIndex(vlayerB)
#vproviderA.rewind()
nElement = 0
selectionA = QGisLayers.features(vlayerA)
nFeat = len(selectionA)

index = utils.createSpatialIndex(layerB)

selectionA = QGisLayers.features(layerA)

current = 0
total = 100.0 / float(len(selectionA))

for inFeatA in selectionA:
nElement += 1
progress.setPercentage(nElement/float(nFeat) * 100)
geom = QgsGeometry( inFeatA.geometry() )
int_geom = QgsGeometry( geom )
atMap = inFeatA.attributes()
intersects = index.intersects(geom.boundingBox())
found = False
geom = QgsGeometry(inFeatA.geometry())
attrs = inFeatA.attributes()
intersections = index.intersects(geom.boundingBox())
first = True
for id in intersects:
vlayerB.featureAtId(int(id), inFeatB)
tmpGeom = QgsGeometry(inFeatB.geometry())
if tmpGeom.intersects(geom):
found = True
if first:
outFeat.setGeometry(QgsGeometry(tmpGeom))
first = False
else:
found = False
if len(intersections) > 0:
for i in intersections:
request = QgsFeatureRequest().setFilterFid(i)
inFeatB = layerB.getFeatures(request).next()
tmpGeom = QgsGeometry(inFeatB.geometry())
if tmpGeom.intersects(geom):
found = True
if first:
outFeat.setGeometry(QgsGeometry(tmpGeom))
first = False
else:
try:
cur_geom = QgsGeometry(outFeat.geometry())
new_geom = QgsGeometry(cur_geom.combine(tmpGeom))
outFeat.setGeometry(QgsGeometry(new_geom))
except:
GEOS_EXCEPT = False
break
if found:
try:
cur_geom = QgsGeometry(outFeat.geometry())
new_geom = QgsGeometry(geom.intersection(cur_geom))
if new_geom.wkbType() == QGis.WKBNoGeometry :
int_com = QgsGeometry(geom.combine(cur_geom))
int_sym = QgsGeometry(geom.symDifference(cur_geom))
new_geom = QgsGeometry(int_com.difference(int_sym))
try:
cur_geom = QgsGeometry( outFeat.geometry() )
new_geom = QgsGeometry( cur_geom.combine( tmpGeom ) )
outFeat.setGeometry( QgsGeometry( new_geom ) )
outFeat.setGeometry(new_geom)
outFeat.setAttributes(attrs)
writer.addFeature(outFeat)
except:
GEOS_EXCEPT = False
break
if found:
try:
cur_geom = QgsGeometry( outFeat.geometry() )
new_geom = QgsGeometry( geom.intersection( cur_geom ) )
if new_geom.wkbType() == 7:
int_com = QgsGeometry( geom.combine( cur_geom ) )
int_sym = QgsGeometry( geom.symDifference( cur_geom ) )
new_geom = QgsGeometry( int_com.difference( int_sym ) )
try:
outFeat.setGeometry( new_geom )
outFeat.setAttributes( atMap )
writer.addFeature( outFeat )
FEAT_EXCEPT = False
continue
except:
FEATURE_EXCEPT = False
GEOS_EXCEPT = False
continue
except:
GEOS_EXCEPT = False
continue

current += 1
progress.setPercentage(int(current * total))

del writer

if not GEOS_EXCEPT:
SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Geometry exception while computing clip")
if not FEATURE_EXCEPT:
SextanteLog.addToLog(SextanteLog.LOG_WARNING, "Feature exception while computing clip")

def defineCharacteristics(self):
self.name = "Clip"
self.group = "Vector overlay tools"
self.addParameter(ParameterVector(Clip.INPUT, "Input layer", ParameterVector.VECTOR_TYPE_ANY))
self.addParameter(ParameterVector(Clip.INPUT2, "Clip layer", ParameterVector.VECTOR_TYPE_ANY))
self.addOutput(OutputVector(Clip.OUTPUT, "Clipped"))

0 comments on commit aa4f70c

Please sign in to comment.