Like I said, not for the faint hearted, and probably far from pythonic, but maybe you'll get some ideas. It was designed to handle fairly large datasets (20000-ish road segments) in a few minutes. Tried to remove most of the bells and whistles, so hopefully it's generic enough to apply to most situations. Run it in QGIS Python Console (find it under the Plugins dropdown).
#Sample pyqgis code to select/highlight features in the 'new' layer that do not have a corresponding match to the the 'old' layer
import itertools
tolrDist=2 #Distance between old and new features that is tolerable
lddLrs = qgis.utils.iface.legendInterface().layers() #get all loaded layers
for lyr in lddLrs:
if (lyr.type()==QgsMapLayer.VectorLayer and lyr.name()[-24:]=='newLyrName'): #change newLyrName to the newer/upper layer in the legend
newLyr = lyr
if (lyr.type()==QgsMapLayer.VectorLayer and lyr.name()[-24:]=='oldLyrName'): #change oldLyrName to the older/lower layer in the legend
oldAdLyr = lyr
##Get all the features to start
try:
featDict
except NameError:
featDict={}
featDict['old'] = {feature.id(): feature for (feature) in oldLyr.getFeatures()} #old Features
featDict['new'] = {feature.id(): feature for (feature) in newLyr.getFeatures()} #new Features
##Load features into QgsSpatialIndex (speeds up processing time)
try:
fidxDict
except NameError:
fidxDict={}
fidxDict['old'] = QgsSpatialIndex(oldLyr.getFeatures(QgsFeatureRequest().setFilterFids(featDict['old'].keys()))); #old Index
fidxDict['new'] = QgsSpatialIndex(newLyr.getFeatures(QgsFeatureRequest().setFilterFids(featDict['new'].keys()))); #new Index
print r"{0} features in older layer".format(len(featDict['old'].keys())
print r"{0} features in newer layer".format(len(featDict['new'].keys())
##***Finally do the intersection check
#Absolute check (detects changes to floating point precision)
actualNewFeatures=[] #place to store ids of all the 'new' features that do not overlap older ones. Note that this doesn't yet account for minor coordinate changes.
idxIntersector = lambda fidx, ftDict: map(fidx.intersects, [ft.geometry().boundingBox() for ft in ftDict.itervalues()]) #function to check whether feature intersects with any features in QgsSpatialIndex
newLyr.removeSelection() #clears selection from newLyr
if newLyr.geometryType() == QGis.Line:
for i, ids in itertools.izip(featDict['new'].keys(), idxIntersector(fidxDict['old'], featDict['new'])): #checks if new features' bounding box intersects that of the old features
if sum([featDict['old'][id].geometry().equals(featDict['new'][i].geometry()) for id in ids])==0:
actualNewFeatures.append(i) #if there are no intersections, that means it's a new feature!
newLyr.selectByIds(actualNewFeatures) #selects the list of all new features
#Fuzzy-matching check
fuzzyNewFeatures=[] #place to store ids of all the fuzzy 'new' features that do not overlap older ones. Account for minor coordinate changes as set by tolrDist.
newOldPair={} #Dictionary to store all new feature's (key) closest old pair (value)
for f in newLyr.selectedFeatures():
nids = fidxDict[sName]['old'].intersects(f.geometry().boundingBox())
try:
distanceToClosestOldFeat = sorted(enumerate(map(f.geometry().distance, [featDict[sName]['old'][nid].geometry() for nid in nids])))[0][1]
if distanceToClosestOldFeat > tolrDist:
fuzzyNewFeatures.append(f.id())
except IndexError: #error is raised when we can't find a close enough pair within the bounding box, probably means it's far enough
fuzzyNewFeatures.append(f.id())
newLyr.removeSelection() #clears selection from newLyr
newLyr.selectByIds(fuzzyNewFeatures) #selects the list of fuzzy new features
After the script finishes, you should be able to see the highlighted new features.