6

I've been looking on how to get GEOSExtractLine available in python.

I looked into Shapely source to if it was available for it, but it's not. I'm trying to wrap this funcionality into Django, much like other linear referencing methods (https://code.djangoproject.com/ticket/11948).

I changed the code properly but I'm getting a "GEOSExtractLine" not found error when trying to load the function. Here's the traceback:

  File "C:\Python27\lib\site-packages\django\contrib\gis\geos\__init__.py", line 6, in <module>
    from django.contrib.gis.geos.geometry import GEOSGeometry, wkt_regex, hex_regex   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\geometry.py", line 15, in <module>
    from django.contrib.gis.geos.coordseq import GEOSCoordSeq   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\coordseq.py", line 10, in <module>
    from django.contrib.gis.geos import prototypes as capi   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\prototypes\__init__.py", line 13, in <module>
    from django.contrib.gis.geos.prototypes.geom import from_hex, from_wkb, from_wkt, \   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\prototypes\geom.py", line 134, in <module>
    geos_extract_line = geom_output(GEOSFunc('GEOSExtractLine'),[c_double,c_double])   File "C:\Python27\lib\site-packages\django\contrib\gis\geos\prototypes\threadsafe.py", line 39, in __init__
    self.cfunc = getattr(lgeos, func_name)   File "C:\Python27\lib\ctypes\__init__.py", line 366, in __getattr__
    func = self.__getitem__(name)   File "C:\Python27\lib\ctypes\__init__.py", line 371, in __getitem__
    func = self._FuncPtr((name_or_ordinal, self)) AttributeError: function 'GEOSExtractLine' not found

How do I find out the correct name that GEOS publishes for this method? I tried to follow the convention used by Shapely and Django's GEOS wrapper, found here http://geos.osgeo.org/doxygen/classgeos_1_1linearref_1_1LengthIndexedLine.html#a20e3c3e6faa7183e23e3f9a312893aca

What I'm truly after is this:

def test_extract_line(self):
    line = LineString((0,0),(10,0))
    extract = line.extract_line(.5,1)

    from_point = Point(5,0)
    to_point = Point(10,0)

    self.assertEqual(extract[0][0],from_point[0])
    self.assertEqual(extract[0][1],from_point[1])
    self.assertEqual(extract[1][0],to_point[0])
    self.assertEqual(extract[1][1],to_point[1])

Any tips?

Sean Gillies come to rescue? Help! :D

EDIT: based on this ticket https://code.djangoproject.com/ticket/11948 and the patches it describes I've added the following code to the respective files:

geometry.py (the actual wrapper):

def extract_line(self,from_distance,to_distance):
    if not isinstance(self,(LineString,MultiLineString)):
        raise TypeError('extract line only works on LineString and MultiLineString geometries')
    return capi.geos_extract_line(from_distance,to_distance)

geom.py (where I point to geos ctype library)

geos_extract_line = geom_output(GEOSFunc('GEOSExtractLine'),[c_double,c_double])

__init__.py

I've added the import instruction for geos_extract_line.

When I do this, I get the traceback written above "function GEOSEXtractLine not found".

Digging a bit further on GEOS docs, it looks like the correct function I'm after it's called Extract Line By Location (http://geos.osgeo.org/doxygen/classgeos_1_1linearref_1_1ExtractLineByLocation.html)

Lastest EDIT:

My guess is that this function is not published by GEOS capi, because:

>>> from django.contrib.gis.geos import libgeos
>>> libgeos.lgeos
<CDLL 'C:\OSGeo4W\bin\geos_c.dll', handle 4200000 at 2ad9870>
>>> dir(libgeos.lgeos)
['GEOSArea_r', 'GEOSBoundary_r', 'GEOSBuffer_r', 'GEOSContains_r', 'GEOSConvexHull_r', 'GEOSCoordSeq_clone_r', 'GEOSCoordSeq_create_r', 'GEOSCoordSeq_getDimensions_r', 'GEOSCoordSeq_getOrdinate_r', 'GEOSCoordSeq_getSize_r', 'GEOSCoordSeq_getX_r', 'GEOSCoordSeq_getY_r', 'GEOSCoordSeq_getZ_r', 'GEOSCoordSeq_setOrdinate_r', 'GEOSCoordSeq_setX_r', 'GEOSCoordSeq_setY_r', 'GEOSCoordSeq_setZ_r', 'GEOSCrosses_r', 'GEOSDifference_r', 'GEOSDisjoint_r', 'GEOSDistance_r', 'GEOSEnvelope_r', 'GEOSEqualsExact_r', 'GEOSEquals_r', 'GEOSFree_r', 'GEOSGeomFromHEX_buf_r', 'GEOSGeomFromWKB_buf_r', 'GEOSGeomFromWKT_r', 'GEOSGeomToHEX_buf_r', 'GEOSGeomToWKB_buf_r', 'GEOSGeomToWKT_r', 'GEOSGeomTypeId_r', 'GEOSGeomType_r', 'GEOSGeom_clone_r', 'GEOSGeom_createCollection_r', 'GEOSGeom_createLineString_r', 'GEOSGeom_createLinearRing_r', 'GEOSGeom_createPoint_r', 'GEOSGeom_createPolygon_r', 'GEOSGeom_destroy_r', 'GEOSGeom_getCoordSeq_r', 'GEOSGeom_getDimensions_r', 'GEOSGetCentroid_r', 'GEOSGetExteriorRing_r', 'GEOSGetGeometryN_r', 'GEOSGetInteriorRingN_r', 'GEOSGetNumCoordinates_r', 'GEOSGetNumGeometries_r', 'GEOSGetNumInteriorRings_r', 'GEOSGetSRID_r', 'GEOSHasZ_r', 'GEOSInterpolateNormalized_r', 'GEOSInterpolate_r', 'GEOSIntersection_r', 'GEOSIntersects_r', 'GEOSLength_r', 'GEOSLineMerge_r', 'GEOSNormalize_r', 'GEOSOverlaps_r', 'GEOSPointOnSurface_r', 'GEOSPrepare_r', 'GEOSPreparedContainsProperly_r', 'GEOSPreparedContains_r', 'GEOSPreparedCovers_r', 'GEOSPreparedGeom_destroy_r', 'GEOSPreparedIntersects_r', 'GEOSProjectNormalized_r', 'GEOSProject_r', 'GEOSRelatePattern_r', 'GEOSRelate_r', 'GEOSSetSRID_r', 'GEOSSimplify_r', 'GEOSSymDifference_r', 'GEOSTopologyPreserveSimplify_r', 'GEOSTouches_r', 'GEOSUnionCascaded_r', 'GEOSUnion_r', 'GEOSWKBReader_create_r', 'GEOSWKBReader_destroy_r', 'GEOSWKBReader_readHEX_r', 'GEOSWKBReader_read_r', 'GEOSWKBWriter_create_r', 'GEOSWKBWriter_destroy_r', 'GEOSWKBWriter_getByteOrder_r', 'GEOSWKBWriter_getIncludeSRID_r', 'GEOSWKBWriter_getOutputDimension_r', 'GEOSWKBWriter_setByteOrder_r', 'GEOSWKBWriter_setIncludeSRID_r', 'GEOSWKBWriter_setOutputDimension_r', 'GEOSWKBWriter_writeHEX_r', 'GEOSWKBWriter_write_r', 'GEOSWKTReader_create_r', 'GEOSWKTReader_destroy_r', 'GEOSWKTReader_read_r', 'GEOSWKTWriter_create_r', 'GEOSWKTWriter_destroy_r', 'GEOSWKTWriter_write_r', 'GEOSWithin_r', 'GEOSisEmpty_r', 'GEOSisRing_r', 'GEOSisSimple_r', 'GEOSisValidReason_r', 'GEOSisValid_r', 'GEOSversion', '_FuncPtr', '__class__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattr__', '__getattribute__', '__getitem__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_func_flags_', '_func_restype_', '_handle', '_name', 'finishGEOS_r', 'initGEOS_r']

With this I've got a glimpse of all the functions published by GEOS API and there is no Extract, ExtractLine or extractByIndex, or anything.

In fact, if you look at CAPI (the published API) there is truly no reference of these functions:

http://trac.osgeo.org/geos/browser/trunk/capi/geos_c.h.in

Still, if anyone can help, I appreciate.

Another edit: PostGIS has the funcionality. How did they do it?

http://postgis.net/docs/ST_Line_Substring.html

Thanks!

8
  • 1
    I have not attempted this (that's why I don't put this as an answer), but after reading a comment in the GEOS source code svn.osgeo.org/geos/trunk/swig/python/README.txt I assume that you can just add a few ctype declarations docs.python.org/library/ctypes.html in the shapely code to add the methods that are missing. Try it here first: github.com/sgillies/shapely/blob/master/shapely/… Mar 26, 2012 at 3:24
  • I actually tried to do that! But not with the shapely code, with GeoDjango code. I did not posted the full code here because it spans a few classes, but I will. Mar 26, 2012 at 3:26
  • ExtractLine is not in geos_c.h, which means you can't access these methods from Shapely, or other C API clients (PostGIS, etc.)
    – Mike T
    Mar 26, 2012 at 8:09
  • 1
    @MikeToews You are right. George: if you really want that particular function (or any other ones that are not in the official C API), you'll have to expose them. basically you will need to edit svn.osgeo.org/geos/trunk/capi/geos_c.h.in and svn.osgeo.org/geos/trunk/capi/geos_c.cpp as well as the file that contains the implementation and recompile. That means expose first to the C-API in GEOS, then expose to Shapely, then use it. Not sure if this is worth the trouble, but may be handy for other similar problems in the future. Mar 26, 2012 at 16:42
  • 1
    @George: PostGIS's ST_Line_Substring does not use GEOS; it uses their own light-weight geometry library, liblwgeom. It is implemented in lwgeom_functions_analytic.c.
    – Mike T
    Mar 26, 2012 at 22:05

1 Answer 1

7

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)]
6
  • Hi Mike! That'a great news. In Django the coordinates are slicable too. What I'm truly after, unless I understood correctly, is that extract line will let me give a certain from-to distance and it will return me the sliced line, not entirely based on the vertexes. Mar 25, 2012 at 21:35
  • Ah, so like extract after the first 10 metres, to the point before 10 metres at the end, like line.interpolate(10) and line.interpolate(-10). It is possible to make a python function to extract the correct subset of coordinates, then modify the start and end coordinates.
    – Mike T
    Mar 25, 2012 at 21:42
  • Hi Mike. Exactly. I just think this would be too much work. I'll take a stab at it :D AFAIK the extract line functionality is present in GEOS. Mar 26, 2012 at 3:09
  • @George, I've updated my answer with a pure Python function, which should work with your example.
    – Mike T
    Mar 26, 2012 at 4:26
  • great. Nice function :D Mar 26, 2012 at 11:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.