4

I'm trying to simply get the Difference between two sets of lines, but the algorithm returns a copy of either one, depending of the combination of INPUT/DIFFERENCE layer I try. Never the difference.

Both of my shapefiles are in the same projection, and have been cleaned (v.clean).

I tried to dissolve them, but no luck.

Any ideas?

(Using QGIS 2.18.1 over Win7)


Here is a picture of the situation. I have 2 sets of roads (red and green). The green network overlaps everything from the red one, but has some other roads in it. I'm trying to get a shape of these extra segments (hence, the Difference). But the result gives me either an exact copy of the red or green network. enter image description here

  • Difference is meant to return the segments of a line which is not intersecting/overlapping part of another line (see e.g. postgis.net/docs/manual-2.0/ST_Difference.html). Is that what you're after or another overlay function? Could with have a look at a screenshot of your lines? – weiji14 Dec 22 '16 at 20:33
  • This is indeed what I'm trying to do (sorry for the explanation). I'm adding a picture of the situation right away. – Horizen Dec 22 '16 at 20:48
  • 2
    Can you zoom in really close and check if they are exactly overlapping? Sometimes there might be some small decimal point coordinate differences and that might affect the overlay. – weiji14 Dec 22 '16 at 22:07
  • You are right; when I zoom VERY close, they don't overlap. Is there a way to ask the tool to use a certain amount of tolerance for this? in command line I guess? – Horizen Dec 22 '16 at 22:21
  • The built-in difference tool doesn't have a tolerance option unfortunately. Coincidentally however, I am working on some python code now that is checking for minor changes between new/old road lines. If you're comfortable with python qgis, I could share some of the code :) It's not quite for the faint hearted though. – weiji14 Dec 22 '16 at 22:27
1

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.

0

You could try the Spatial Query plugin to do a select by overlap.

Select the green streets that overlap the red, then switch and save the selected features.

This assumes that the non-overlapping segments are unique features and not extensions of the overlapped features.

  • I'll have to look into it, but I'm afraid they are extensions. – Horizen Dec 22 '16 at 22:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.