Create mapping between overlapping linestrings

Question

What I have is two shapefiles, describing roads in a city. The shapefiles are made of linestrings, each covering a section of road.

• The road segments in shapefile 1 are generally several times larger than the segments in shapefile 2.
• The road segments break at road junctions for both shapefiles.
• The shapefiles do not describe every road in the city.
• Both shapefiles describe the same city, although there is some difference in the roads they cover. Some roads decribed by segments in shapefile 1 are not present in shapefile 2, and vice versa.

What I want is to create a many to one map between the shapefiles. I would like a table that lists each segment in shapefile 2, and it's closest match in shapefile 1. By closest match I mean closest in similarity, not necessarily in distance. In situations where a segment in shapefile 2 overlaps two segments in shapefile 1, the segment it overlaps the most would be chosen. Where there is no matching segment in shapefile 1, the table entry would be null.

I'm not sure how to best approach this task. My initial plan is to create a new polygon layer using the buffer tool in QGIS, from the features in shapefile 1, and then somehow find which buffer the linestrings from shapefile 2 predominantly fall in. I don't know if there are any tools which would help with this second step, or tackles the whole problem. I know how to use QGIS, and Python.

Answer 1

So, I thought about this a bit. Your problem not a conventional one, and thus it requires a non-obvious solution. The problem is, it's impossible to tailor a solution specifically for it without a good knowledge of your data. So, instead of giving you a set solution, I'll try pointing you to a few, but you'll have to try them yourself and see what works.

The problem. Similarity between lines is not easily definable, but given the configuration of an urban network, I could think of three parameters that affect it:

• Azimuth: to avoid one line from being matched to streets perpendicular to it, it's important to filter by the azimuths of the lines. The azimuths can't be a perfect match, because of vectorization (or measuring) imperfections, so a given margin θ is necessary. This margin, however, must not exceed the smallest angle between two connected road segments in your network (if it's a perfect grid, that'd be 90 degrees).

• Distance: to avoid two lines with similar azimuths, but in opposite sides of the city, from matching, distance must be taken into account as well. Since traditional distance algorithms only take into account the two closest points between two features, to avoid two connecting streets from matching to each other, I suggest matching both initial and final vertices. This, however, should have a smaller weight than azimuth, because, given that the two shapefiles have lines of different lengths, it's possible that smaller distance be matched to adjascent streets.

  • Direction: if you have information on both lanes, this must be taken into account when matching distance and azimuths. If not, you can calculate matching azimuths and counter-azimuths for the same line.

Doing this in QGIS is possible, but quite complicated - it simply wasn't made for this sort of multi-step processing. Since you said you're familiar with python, you can try using GDAL (PyQGIS would also be a possibility, but I'm not as familiar with it).

GDAL is already installed with QGIS. All you'll need to do is download and configure the python bindings for it. To open a shapefile in GDAL:

from osgeo import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
data_source = driver.Open(<filepath>)
layer = data_source.GetLayer()

for i in range(layer.GetFeatureCount()):
    feat = layer.GetFeature(i)
    geom = feat.GetGeometryRef()

You can read more about it here. Now, with your data open, you can start doing the checks. Your idea of comparing buffer polygons is not a bad one, and is quite uncomplicated. It may generate incorrect results depending on your network shape, but it's worth giving it a try. You can start by making a search on the neighborhood of your linestring, to limit which features you'll be comparing to. Then checking for intersects:

di_intersect = dict()

buffer_geom1 = geom1.Buffer(radius)
buffer_geom2 = geom2.Buffer(radius)

if buffer_geom1.Intersects(buffer_geom2):
    intersect = buffer_geom1.Intersection(buffer_geom2)
    area_intersect = intersect.Area()
    feat_id = feat2.GetFieldAsString('FID')
    di_intersects[feat_id] = area_intersect

Then you compare which of the features had the largest intersect area, and put that feature's id in a list, or something. Some of these methods might seem confusing if you've never messed with GDAL before, but you can check the documentation here.

If, however, this does not solve the problem, you can try the more laborious way. I can't write the whole code out, but can give you a few snippets to guide your way:

Finding the azimuth (Source)

angle = math.atan2(point2.GetX() - point1.GetX(), point2.GetY() - point1.GetY())
azimuth = math.degrees(angle) if angle > 0 else math.degrees(angle) + 180

Finding closest distance

start1 = geom1.GetPoint(0)
start2 = geom2.GetPoint(0)
end1 = geom1.GetPoint(geom1.GetPointCount() - 1)
end2 = geom2.GetPoint(geom2.GetPointCount() - 1)

dist = start1.Distance(start2) + end1.Distance(end2)

OBS: In the above code, if the features are inverted (i.e. direction doesn't matter), you'll have to match start1 to end2, and start2 to end1.

With this, you can work it in a more complete code, and tweek the search radius to fit your model. The main idea is to iterate through each line in your shapefile 1, and test (either buffer intersection, or azimuth/distance/direction) it for each of the lines in shapefile 2 that are close to it. Then, check which is the best fit, and write that to a list. Final part is to export this list as a CSV, presumibly.

Hope this helps!