I'm trying to match small segments with a larger segment they are the most probably related to: relatively close, similar bearing, and facing each other.

Here is a typical example of the data I have:


Here I'd need to match segment 652 to 198969, while having 711 and 707 not matching anything.

I've looked for different methods, in particular the Hausdorff distance (based on the answers here). I computed it using PostGIS but I'm getting odd results: the shortest distance I get is between 707 and 198985, and 652 has a greater distance to 198969 than to 198985 for example (I can add the query and results if needed).

Is Hausdorff actually the correct method to solve this? Are there other approaches? I thought of simply creating a set of checks on the parameters I mentioned (distance, bearing, etc.) but I'm afraid of having to add a whole bunch of conditions to handle edge cases or things like thresholding on how much they're facing each other.

Update: I found a method which seems like an acceptable compromise:

  • I first find the 10 nearest black segments from the blue one I'm trying to match (using the PostGIS <-> operator) which are less than 10 meters away.
  • I then create a new segment by finding the nearest points to the ends of the blue segment on each of the black ones (using ST_ClosestPoint) and filter out the results whose length is less than 90% of the blue one (meaning the segments aren't facing, or that the bearing difference is more than ~20°)
  • Then I get the first result sorted by distance and Hausdorff distance, if any.

There might be some fine tuning to do but it seems to do an acceptable job for now. Still looking for any other methods or additional checks to run if I missed some edge cases.

  • 1
    Why not just use the endpoints of the (blue) segments to identify potential matches among the target (black) segments? The segment matches you are looking for occur when both endpoints are close to a common target, which is a simple query to execute. This method handles discrepancies that are due to errors in the locations of the segment endpoints, which otherwise might be challenging to cope with: notice that a very short (blue) segment has a much less precise bearing than a longer segment.
    – whuber
    Commented Jan 13, 2017 at 18:11
  • 1
    Yep I'm actually trying something along these lines, I'll update the question with the details.
    – Jukurrpa
    Commented Jan 13, 2017 at 18:22
  • 1
    Not sure if I understand you correctly, but have you tried to create centroids for blue lines and then check the distances from them to the nearest lines, leaving the shortest distance as a result? Commented Jan 8, 2019 at 19:38
  • 1
    Hi Cyril, I'm not working on this problem anymore, but the issue was also to match the blue segments based on their orientation, and how much they are "facing" the black segments. Which means in this case that 711 shouldn't match with anything even though he's pretty close to the black segments
    – Jukurrpa
    Commented Jan 8, 2019 at 22:01

1 Answer 1


Here are a couple functions I wrote that should let you do what you need to do. Check and see if the line is polyline or segment if polyline explode the line, then compare azimuth and inverse of first point and last point on lines, set you acceptable criteria for code to make decision for you. These are arcpy methods but could be modified.

return azimuth from line segment ESRI @shape

def returnAzimuth(shape):
    point1 = shape.firstPoint
    point2 = shape.lastPoint
    dX = point2.X-point1.X
    dY = point2.Y-point1.Y
    az = math.atan2(dX,dY)*180/math.pi
    if az<0:
        az = az+360
        return az
    return az

return inverse from ESRI points

def returnInverse(point1,point2):
    dX = point2.X-point1.X
    dY = point2.Y-point1.Y
    dis = sqrt(dX**2+dY**2)
    az = math.atan2(dX,dY)*180/math.pi
    if az<0:
        az = az+360
        return az,dis
    return az,dis

explode polyline into line segment

def explodePolyline(shape):
    sr = "insert your spatial reference"
    pnts = shape.getPart(0)
    for x in range(len(pnts)-1):
        array = arcpy.Array()
        point1 = arcpy.Point(pnts.getObject(x).X,pnts.getObject(x).Y,pnts.getObject(x).Z)
        point2 = arcpy.Point(pnts.getObject(x+1).X,pnts.getObject(x+1).Y,pnts.getObject(x+1).Z)
        line = arcpy.Polyline(array,sr,True,True)
    return lines

run through your table like this

for shape in Table:
    for shape2 in Table:
         check if shape is polyline:
            if yes, explode and compare returned line segments with shape2
                for seg in returnedlinesegments
                    if seg.firstpoint=shape2.lastpoint and shape.az=shape2.az
                        do stuff.....
            if no, compare shape and shape2
                if shape.firstpoint=shape2.lastpoint and shape.az=shape2.az
                   do stuff.....

these properties should be available in postgis - firstpoint, lastpoint, pointarray I assume esri properties above because it is what I know best but the above can easily be changed to work with postgis.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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