Thursday 29 December 2016

python - Intersecting MultiLinestring based geodataframes with geopandas?


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.




Answer



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)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...