As of Shapely version 1.2.14, coordinates are slicable. This looks very similar to GEOSExtractLine, where a subset of the LineString can be extracted.
Here are some examples how you can slice coordinates to extract a new line object:
from shapely.geometry import LineString, Point
# Original LineString used for examples
line = LineString([(30, 50), (60, 100), (100, 70), (150, 120), (190, 80)])
line.length # 235.58873956203155
# Extract a LineString from first two coordinates
LineString(line.coords[:2]).length # 58.309518948453004
# .. and from the last two coordinates
LineString(line.coords[-2:]).length # 56.568542494923804
# Everything but the first and last coordinate
LineString(line.coords[1:-1]).length # 120.71067811865476
# Every second coordinate
LineString(line.coords[::2]).length # 163.35495027417937
# Reversed/flipped
LineString(line.coords[::-1]).length # 235.58873956203155
etc.
As for the general topic of getting GEOS functionality to Shapely, that is indeed a @sgillies answer. Generally, the functionality needs to also be supported in the GEOS C API.
To extract a new LineString using linear referencing methods, here is an add-on function in pure Python. If normalized=True, then uses normalized start/end linear referencing (i.e., within [0.0, 1.0]); otherwise the default is absolute linear referencing (i.e., within [0.0, total length]). Negative end is supported as the reverse distance from the end of the linestring.
def extract_line(line, start, end, normalized=False):
"""Extract LineString from `start` to `end`"""
assert start >= 0, '`start` must be greater than or equal to zero'
if end < 0:
if normalized:
end += 1.0
else:
end += line.length
assert end > start, '`end` must be greater than `start`'
start_c = line.interpolate(start, normalized).coords[0]
end_c = line.interpolate(end, normalized).coords[0]
for start_i, c in enumerate(line.coords):
if c == start_c:
start_c = None
break
if line.project(Point(c), normalized) >= start:
break
for end_i, c in enumerate(line.coords):
if c == end_c:
end_c = None
end_i += 1
break
if line.project(Point(c), normalized) >= end:
break
extracted = line.coords[start_i:end_i]
if start_c:
extracted.insert(0, start_c)
if end_c:
extracted.append(end_c)
return LineString(extracted)
Examples:
# Your example in the question:
extract_line(LineString([(0,0), (10,0)]), 0.5, 1, True).coords[:]
# [(5.0, 0.0), (10.0, 0.0)]
# My example, length of 10.0 removed from each end:
extract_line(line, 10, -10).coords[:]
# [(35.144957554275265, 58.57492925712544), (60.0, 100.0),
# (100.0, 70.0), (150.0, 120.0), (182.92893218813452, 87.071067811865476)]
It depends on your needs. If you need to display a relatively small number of relatively simple features (such as all locations for shop) then you can eliminate Geoserver and simply use GeoDjango as an ersatz web feature server by serializing data into kml or geojson and serving it to OpenLayers for display, you can also do plenty of CRUD operations easily in GeoDjango.
If you have relatively larger and complex data (such as land parcels or building outlines for a city for example), then you will want the 'proper' WFS features that you get with Geoserver. You can still use GeoDjango to create/query/update/ data in the db, but GeoServer would operate entirely separately from GeoDjango, using data from the same db to fulfill requests coming in from OpenLayers.
For me GeoDjango's biggest strength is being able to use spatial queries with basically the same great ORM syntax as Django, this makes it incredibly easy to write and test complex spatial operations that perform well. The ability to chain queries together is superb. Although you can use it as a simple feature service, it would be difficult to make it perform well under heavy use compared to a purpose-built WFS like Geoserver.
Best Answer
Here's a potential answer (at least, the one I've chosen to use).
My problem comes from the fact that current geodjango API does not use linearref features from GEOS. So, I've patched my django code (using linked patch in my question) to get the 'project()' and 'interpolate()' methods.
Then, I simply loop to interpolate the needed points. The key here is to use the linearref feature of GEOS:
http://geos.osgeo.org/doxygen/classgeos_1_1linearref_1_1LengthIndexedLine.html#7f1fd2d3467a3bae6e28e2da893506f9