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:
unary_union. Returns a geometry containing the union of all geometries in the GeoSeries .
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])
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.
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)
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)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With