3

I would like to perform the following operation using geopandas.

enter image description here

The segments are delimited by red points and the blue items are attribute information. My inputs are the first and second line segments and my output is the third line segment.

Initially I thought this would be an intersection operation, but I soon learned that geopandas can only intersect polygons, therefore something like:

intersection = geopandas.overlay(split_lines, original_lines, how='intersection')

returns the following error:

raise TypeError("overlay only takes GeoDataFrames with (multi)polygon "
TypeError: overlay only takes GeoDataFrames with (multi)polygon

This to me looks like a standard geoprocessing operation and I really hope I won't have to code this up from scratch. Are there any simplified ways to come up with the following result without having to code a custom function?

EDIT Unfortunately I cannot get this to work for more complicated geometries such as

ORIGINAL_LINES

                                            geometry property
0  (LINESTRING (0 0, 6.656423206909781 4.43757029...        a
1  (LINESTRING (6.656423206909781 4.4375702913320...        b
2  (LINESTRING (8.070636769282876 5.8517838537051...        c
3  (LINESTRING (6.656423206909781 4.4375702913320...        d
4  (LINESTRING (10.98655022583197 1.9375702913320...        e
5  (LINESTRING (13.68293236472948 0.6224568509648...        a
6  (LINESTRING (17.54663566988575 -0.412819329445...        a

SPLIT_LINES

                                            geometry  susc
0  LINESTRING (0 0, 4.160264504318614 2.773481432...     1
1  LINESTRING (4.160264504318614 2.77348143208253...     2
2  LINESTRING (6.656423206909781 4.43757029133205...     3
3  LINESTRING (9.950815132334437 8.18950268263064...     4
4  LINESTRING (13.08444573742037 12.0857007308397...     5
5  LINESTRING (6.656423206909781 4.43757029133205...     4
6  LINESTRING (10.98655022583197 1.93757029133205...     3
7  LINESTRING (15.61478401730761 0.10481876075978...     2

The output seem to be a 1D line... which is incorrect for this application.

  • 1
    Have you tried sjoin() with op='intersection'? I can't remember if it works for multi-linestrings but I remember it would often work for me when overlay failed. – Jon Sep 27 '18 at 13:55
  • 1
    Are you working with actual 1-D lines like your illustration? – Jon Sep 27 '18 at 14:06
  • 1
    I don't know if there's an easy gpd solution, but you could just make a Nx2 numpy array; first column is the x-coordinate, second column is the segment ID. Then just sort the array rows according to the first column. You'd need to do a little interfacing with shapely to get coordinates and re-build the final line, but that's pretty easy. – Jon Sep 27 '18 at 14:10
  • 1
    Good point. I would make three columns, then. First one is x-coordinate, second one is line 1 segment id, third one is line 2 segment id. Then sort based on first row. Then loop through first column, and for each x-coordinate, check if the segment id has changed for both attributes. If not, store the last-used attribute, and if so, store the new attribute. The second two columns could have nan's wherever they don't change. Sorry if that's confusing. Maybe someone can still help with an easier solution. – Jon Sep 27 '18 at 14:16
  • 1
    stackoverflow.com/questions/2828059/… has some solutions; Steve's is the most compact but difficult to understand. – Jon Sep 27 '18 at 15:37
4

EDIT: The following answer only works for the 1-D case. To extend it to 2-D, you will need to parameterize your links by the along-road distance, and replace the x-coordinate with the parameterized length. However, I'm fairly confident this is doable much simpler with Geopandas.

It would be too hard to give hints in the comments, so here's a script that should give you what you want. It's not written for efficiency--probably you could get geopandas to do what you want with some finegaling, but here ya go. It's also not written very generally, but that could be done if you have more than one attribute.

import geopandas as gpd
from shapely.geometry import LineString
import numpy as np

slgeom = [[(0,0),(7,0)],[(7,0),(13,0)],[(13,0),(15,0)],[(15,0),(19,0)]]
geoms = []
for s in slgeom:
    geoms.append(LineString(s))
properti = ['a','b','c','d']
split_lines = gpd.GeoDataFrame(geometry=geoms)
split_lines['property'] = properti

olgeom = [[(0,0),(5,0)],[(5,0),(7,0), (10,0)],[(10,0),(13,0),(15,0)],[(15,0),(19,0)]]
geoms = []
for o in olgeom:
    geoms.append(LineString(o))
susc = [1,2,3,4]
original_lines = gpd.GeoDataFrame(geometry=geoms)
original_lines['susc'] = susc


# Do split lines
xs1 = []
attrib1 = []
for g, a in zip(split_lines.geometry.values, split_lines.property.values):
    x = g.coords.xy[0].tolist()
    xs1.extend(x)
    try:
        attrib1[-1] = a
    except:
        pass
    attrib1.extend([a for l in range(len(x))])

# Do originals
xs2 = []
attrib2 = []
for g, a in zip(original_lines.geometry.values, original_lines.susc.values):
    x = g.coords.xy[0].tolist()
    xs2.extend(x)
    try:
        attrib2[-1] = a
    except:
        pass
    attrib2.extend([a for l in range(len(x))])

# Create numpy array for sorting
allxs = list(set(xs1 + xs2))
x_forsort = []
a1_forsort = []
a2_forsort = []
for x in allxs:
    try:
        idx = xs1.index(x)
        a1_forsort.append(attrib1[idx])
    except:
        a1_forsort.append(None)
    try:
        idx = xs2.index(x)
        a2_forsort.append(attrib2[idx])
    except:
        a2_forsort.append(None)
forsort = np.transpose(np.array([allxs, a1_forsort, a2_forsort]))

# Now sort based on x value (1st column)
sorteds = forsort[forsort[:,0].argsort()]

# Work through the sorted lists to create segments with the appropriate attributes
# Store results in a dictionary
output = dict()
output['geometry'] = [] 
output['attrib_1'] = []
output['attrib_2'] = []
for i in range(len(sorteds)-1):

    # Store shapely linestring
    output['geometry'].append(LineString([(sorteds[i,0],0),(sorteds[i+1,0],0)]))

    # Store attributes
    if i == 0:
        output['attrib_1'].append(sorteds[i,1])
        output['attrib_2'].append(sorteds[i,2])
    else:
        if sorteds[i,1] is None:
            output['attrib_1'].append(output['attrib_1'][-1])
        else:
            output['attrib_1'].append(sorteds[i,1])

        if sorteds[i,2] is None:
            output['attrib_2'].append(output['attrib_2'][-1])
        else:
            output['attrib_2'].append(sorteds[i,2])

# Convert back to geopandas dataframe
out_gdf = gpd.GeoDataFrame(output)
out_gdf.crs = original_lines.crs

Result:

out_gdf
Out[185]: 
  attrib_1  attrib_2                 geometry
0        a         1    LINESTRING (0 0, 5 0)
1        a         2    LINESTRING (5 0, 7 0)
2        b         2   LINESTRING (7 0, 10 0)
3        b         3  LINESTRING (10 0, 13 0)
4        c         3  LINESTRING (13 0, 15 0)
5        d         4  LINESTRING (15 0, 19 0)
  • Thanks a lot for your effort, but using the example above I am getting an incorrect result. (Please see edit to original question) – user32882 Sep 27 '18 at 16:47
  • @user32882 Try edited code. I get the expected output now. I had to add the "try" statements in the attrib loops. – Jon Sep 27 '18 at 17:09
  • Really nice. It does work I have to admit. Now I have to look through it to make sure I understand how it does the trick. Do you think this would also work for larger GeoDataFrames with LineStrings that are not necessarily contiguous? where intersections happen and what not? – user32882 Sep 27 '18 at 17:22
  • Glad it works! Haha, I was afraid you'd want something more general. With a little modification it could work for gapped lines. I think there's also a much simpler way (might not be as efficient but would be easier to read) to do the whole thing, but I'm too busy to re-code it :P – Jon Sep 27 '18 at 17:35
  • Shame, this won't work for 2D geometries.... it seems only adapted for 1D lines and actually only generates 1D lines for the result. Could this possibly be extended to 2D lines with multiple segments? Please take a look at the edit to original question above – user32882 Sep 28 '18 at 11:27

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.