Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Interpolating every X distance along multiline in shapely

Tags:

python

shapely

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:

enter image description here

like image 927
traggatmot Avatar asked Jan 20 '16 17:01

traggatmot


People also ask

What is Unary_union?

unary_union. Returns a geometry containing the union of all geometries in the GeoSeries .

What is LineString in Python?

Geometry types Point (*args) A geometry type that represents a single coordinate with x,y and possibly z values. LineString ([coordinates]) A geometry type composed of one or more line segments. LinearRing ([coordinates])

What is a shapely geometry?

Shapely is a Python package for set-theoretic analysis and manipulation of planar features using (via Python's ctypes module) functions from the well known and widely deployed GEOS library. GEOS, a port of the Java Topology Suite (JTS), is the geometry engine of the PostGIS spatial extension for the PostgreSQL RDBMS.


2 Answers

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)
like image 199
Mike T Avatar answered Sep 30 '22 15:09

Mike T


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)
like image 35
mgc Avatar answered Sep 30 '22 13:09

mgc