[GIS] Wrap GEOS functionality into python

geodjangogeospythonshapely

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! 😀

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!

Best Answer

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)]
Related Question