Interpolating every X distance along multiline in shapely

Question

If I have a shapely multiline object that contains many lines whose total length each is 50km (when traced from the origin), and I want to interpolate along the multiline every X meters (let's say 100m), returning shapely point objects every 100m, how can I achieve this?

Here is what I have so far, but it only returns one distinct point (when I know it should return several thousand, as tested in ArcMap):

points = []

for x in range(100,50000,100):
    x,y = multiline.interpolate(x).xy
    xy = (x[0],y[0])
    points.append(xy)

trim = list(set(points))

And here is what trim contains:

[(-90.5864707030599, 38.4688716729703)]

If I break the multiline into individual lines, I can get more points (162), but still not the 1000s I should generate:

lines_shape = []

# breaking out individual lines into shapely line strings
for line in multiline:
    lines_shapes.append(gm.asLineString(line))

points = []

# interpolating points along each line
for line in lines_shapes:
    for x in range(100,50000,100):
        points.append(line.interpolate(x))

point_list = []
for point in points:
    x,y = point.xy
    xy = (x[0],y[0])
    point_list.append(xy)

trim = list(set(point_list))

and trim contains:

[(-90.4766827033434, 38.5972414054466), (-90.5461231698478, 38.5000688058116), (-90.4968889030998, 38.5464768729898), (-90.6189469029766, 38.4718782063951), (-90.513144169767, 38.5129150061034), (-90.5810937033111, 38.5323612724354), (-90.5396967035394, 38.5864442726652), (-90.6868857696619, 38.447798672452), (-90.6213193028429, 38.5146844725911), (-90.509758103553, 38.5536796727591), (-90.476078103521, 38.546296205486), (-90.6427893693333, 38.5555682724382), (-90.6217561026614, 38.5147302723649), (-90.630436369842, 38.5000560723107), (-90.4760629697296, 38.5914768059899), (-90.6849423697939, 38.4566312056467), (-90.6538383698443, 38.4599348058277), (-90.4699029698455, 38.5837950057152), (-90.6694679700701, 38.4474470726055), (-90.6829395031549, 38.4198408724397), (-90.6732927030951, 38.4360894058827), (-90.5866617028747, 38.5332390062557), (-90.6961129028875, 38.4531938061303), (-90.5947731029956, 38.4539120730648), (-90.5142317702794, 38.5249724721824), (-90.5903123694197, 38.5034458725098), (-90.5755567699284, 38.5000690063604), (-90.56752770364, 38.5589170059954), (-90.6271141698824, 38.5430544723665), (-90.6201437702182, 38.4550620064875), (-90.6180047030539, 38.5000570723568), (-90.6589975699812, 38.500065872223), (-90.685946302779, 38.4505508056563), (-90.6840849030951, 38.4547842059075), (-90.6854941029699, 38.4134104059469), (-90.5918703027716, 38.5000662058715), (-90.5471115031927, 38.4808140060424), (-90.463451370184, 38.5525954725834), (-90.6780159694608, 38.4460308724081), (-90.5896267029094, 38.5346980060887), (-90.5807949701111, 38.5325790054973), (-90.7004947695174, 38.4727768728354), (-90.5955769035506, 38.5373834725497), (-90.6569789697031, 38.4542776061077), (-90.6253669031594, 38.521868606032), (-90.6565429028321, 38.4538942062337), (-90.4636589696851, 38.552869806277), (-90.634338969767, 38.546211672811), (-90.4360581033981, 38.5162080059754), (-90.478937169714, 38.604583072231), (-90.6238845029668, 38.5403794056656), (-90.5115069701699, 38.5170084061937), (-90.6607629031848, 38.4788458064697), (-90.6182993029698, 38.4996076730365), (-90.6795531698423, 38.4464060057138), (-90.5715115033526, 38.5319224725247), (-90.6601255698386, 38.4660196727507), (-90.5944027028219, 38.4399328061353), (-90.6658113032484, 38.5630904727274), (-90.657297702926, 38.4486316730944), (-90.6076871697265, 38.5411054058703), (-90.5189131030566, 38.4958198724699), (-90.4633435036991, 38.5571938060298), (-90.4550633700778, 38.5794588057595), (-90.6551231026484, 38.5026692062996), (-90.6304847695558, 38.4814012057819), (-90.6293281694588, 38.5510622058365), (-90.5337045036593, 38.5480590727484), (-90.5636725699325, 38.5031146728847), (-90.5893179035978, 38.5349642063136), (-90.5502847700255, 38.5531414725813), (-90.6577657029242, 38.4673132728702), (-90.6842995029191, 38.45499127301), (-90.596995102941, 38.536759005804), (-90.4590855033927, 38.579707805551), (-90.618116169525, 38.5138608060121), (-90.662823570151, 38.4436198061903), (-90.5324573031614, 38.5290608063917), (-90.4761187699646, 38.5460198060497), (-90.4365065035716, 38.5163978060946), (-90.5993421032549, 38.4517628723399), (-90.6580597695421, 38.4672556730918), (-90.4697125698785, 38.5713498724789), (-90.6246789694582, 38.5225716726252), (-90.5623807028065, 38.539724406037), (-90.5255411029389, 38.5592538722488), (-90.6779729027267, 38.4309624725214), (-90.5757351702413, 38.5712982055281), (-90.4349003036055, 38.5470764060318), (-90.4881969699136, 38.5968882057062), (-90.5654833701608, 38.5641550721791), (-90.6243467697869, 38.5407896062358), (-90.6833993696845, 38.4554192064123), (-90.5125843031225, 38.5353802722622), (-90.4343933036075, 38.5320694055591), (-90.5949243033133, 38.4437996058484), (-90.6531879693487, 38.5465062727268), (-90.4678253704438, 38.5769334060321), (-90.6560765692776, 38.4512860059256), (-90.6705015033411, 38.4387672065226), (-90.4510205703424, 38.5293378056758), (-90.5965189694742, 38.453478006286), (-90.5158851037139, 38.5293438059525), (-90.5153219032826, 38.5000608728918), (-90.571685370183, 38.5322450728324), (-90.5907731028958, 38.5035282063425), (-90.562387970228, 38.5396736725824), (-90.5632127699527, 38.5000698058576), (-90.4351277034801, 38.5467304062636), (-90.6846297033974, 38.4247822729648), (-90.4300920979354, 38.5121736658553), (-90.6077071031996, 38.5414430059705), (-90.700472969951, 38.4725166728872), (-90.6083531698623, 38.5440140058222), (-90.6611483697009, 38.478763472637), (-90.5148001696915, 38.5340326723505), (-90.6261891029441, 38.4987850056043), (-90.6077131034763, 38.5415200724739), (-90.5714903701838, 38.5588274721905), (-90.487086969987, 38.5953794059149), (-90.6309363029674, 38.4813284731115), (-90.4375485030603, 38.5798804061361), (-90.4348689037762, 38.5318386727971), (-90.4636395029601, 38.5528440056267), (-90.5537693695601, 38.5918534726398), (-90.686189369743, 38.45065887269), (-90.5137219032426, 38.5459956061927), (-90.4708333032137, 38.5588582056221), (-90.562663570167, 38.5400910056759), (-90.5898827695735, 38.4573650056716), (-90.6541333699584, 38.4598728056666), (-90.5201981695093, 38.5275174726339), (-90.5279559031418, 38.5197458058049), (-90.4969485029704, 38.5468012728409), (-90.6963449693442, 38.4523242057786), (-90.5869187695849, 38.5329812056987), (-90.7090465028932, 38.4391836726696), (-90.4722723695735, 38.5676876055788), (-90.6076811029, 38.541104472374), (-90.6281757035453, 38.4421816724275), (-90.5341799698291, 38.5482110060135), (-90.5948001033415, 38.494925272764), (-90.5864707030599, 38.4688716729703), (-90.5419337032772, 38.5685326058754), (-90.7019565698392, 38.4481112057488), (-90.6253121695204, 38.521868606032), (-90.5138205031131, 38.5456764727715), (-90.5529913031055, 38.455285806376), (-90.5293321032933, 38.5986526054135), (-90.6309013697018, 38.5463032723585), (-90.4975365696563, 38.6037882720911), (-90.6184813032689, 38.4720362730366), (-90.6781275699308, 38.4462692727899), (-90.5651579702635, 38.4796824727488), (-90.6532305027848, 38.5024180058671), (-90.6545907030978, 38.5581974062657), (-90.5950401027179, 38.4950944055632), (-90.6468105700511, 38.5573522728695), (-90.6698367029009, 38.4476676727061), (-90.4501431701705, 38.5848148728883), (-90.577996903134, 38.503304272455), (-90.6713375697724, 38.545723072942)]

Maybe the problem is in the argument I am passing to interpolate? How do I know what units it expects/is assuming the distance is in?

Edit:

Definition of a multiline - this code opens up a geojson feature collection object in a json file, and then pulls out the specific geojson multilinestring and converts it to a shapely multilinestring object.

from shapely import geometry as gm
import json

file = 'Geometries_GEOJSON.json'

with open(file) as input_file:
    data = json.load(input_file)

geojson_multilines = data['features'][1]['geometry']

multiline = gm.shape(geojson_multilines)

Here is what it looks like:

Before conversion, as the geojson object:

{u'type': u'MultiLineString', u'coordinates': [[[-90.4360581033981, 38.5162080059754], [-90.4356193034874, 38.5151156723311], [-90.4349665036006, 38.5145500059587], [-90.4312643032044, 38.5129080723304], [-90.4300920979354, 38.5121736658553]], [[-90.4348689037762, 38.5318386727971], [-90.4353155700569, 38.5236814729023], [-90.4360581033981, 38.5162080059754]].......

After conversion, as the geojson object:

MULTILINESTRING ((-90.4360581033981 38.5162080059754, -90.4356193034874 38.5151156723311, -90.43496650360061 38.5145500059587, -90.4312643032044 38.5129080723304, -90.4300920979354 38.5121736658553), (-90.43486890377621 38.5318386727971, -90.4353155700569 38.5236814729023, -90.4360581033981 38.5162080059754),........

An example of one such object mapped in ArcMap:

A quick point to add to the answers below - it is good to bear in mind that Shapely intentionally has no notion of units or different coordinate systems, not mentioning geographic coordinates. — Jan Šimbera, Aug 12 '21 at 10:14

score 13 · Answer 1 · answered Jan 26 '16 at 22:29

As @mgc has pointed out, you must project the data to a Cartesian coordinate reference system that has distance units in metres.

I use a function to redistribute the vertices, which uses the interpolate linear referencing method. The distance argument is only approximate, since the total distance of the linestring may not be a multiple of the preferred distance. This function works on only [Multi]LineString geometry types.

from shapely import wkt
from shapely.geometry import LineString

def redistribute_vertices(geom, distance):
    if geom.geom_type == 'LineString':
        num_vert = int(round(geom.length / distance))
        if num_vert == 0:
            num_vert = 1
        return LineString(
            [geom.interpolate(float(n) / num_vert, normalized=True)
             for n in range(num_vert + 1)])
    elif geom.geom_type == 'MultiLineString':
        parts = [redistribute_vertices(part, distance)
                 for part in geom]
        return type(geom)([p for p in parts if not p.is_empty])
    else:
        raise ValueError('unhandled geometry %s', (geom.geom_type,))

# E.g., every 100 m on a projected MultiLineString
multiline_r = redistribute_vertices(multiline, 100)

mgc · Answer 2 · 2016-01-20T21:54:34.630

1

In your first snippet of code you are iterating 150m by 150m but your dataset seems to be in longitude-latitude. You should reproject it in meters (in the appropriate CRS for your zone) as shapely take the distance parameter in the unit of the objet. Exemple taken from the shapely documentation :

from functools import partial
import pyproj

project = partial(
    pyproj.transform,
    pyproj.Proj(init=’espg:4326’),
    pyproj.Proj(init=’epsg:26913’)) # Replace by the appropriate CRS
g2 = transform(project, g1) # g1 is a shapely geometry

You can also avoid the transformation by using normalized value in the interpolate function :

# Progress of 10% of the length:
line.interpolate(0.10, normalized=True)

edited Jan 20 '16 at 21:54

answered Jan 20 '16 at 21:48

mgc

5,223
1
24
37

I projected the multiline before applying the rest of the code, and the results are still the same. Trying normalization now. – traggatmot Jan 20 '16 at 22:48
In addition I saw that in your second snippet of code (when doing `for line in lines_shapes: for x in range(100,50000,100): points.append(line.interpolate(100))` your forget to pass the x value to interpolate, so you arent progressing on the line. – mgc Jan 20 '16 at 23:08
Changed that, but it doesn't change the results. So, now, using a projected multiline as a starting point and correctly using X in the for loop, but no difference in the results. – traggatmot Jan 21 '16 at 04:37

Interpolating every X distance along multiline in shapely

2 Answers2

Linked