Accessing GEOS functions in Python

Question

The issue: I am using Windows 7, 64-bit operating system with Python 3.2. I must continue to use Python 3.2 because of a dependency with another API. I need some functions in GEOS that are not currently supported by GDAL/OGR 1.8 and there is no port of Shapely for Python 3.2 on 64-bit Windows.

Spatialite has a beta version for 64-bits Windows but I can't get it to work and nobody in the main Stackoverflow list has a solution ( https://stackoverflow.com/questions/8590599/sqlite-load-extension-fail-for-spatialite-in-python ).

I don't want to have to go down the route of using PostGIS to get access occassinoally to one or two functions GDAL/OGR doesn't have....

The Question: SO: Can anybody explain how to either get Spatialite working or better still, how I can use CTypes to access the functionality in the GEOS dll (I know how to connect to it but I don't understand how to pass data, say wkt, in and get my result back out - https://stackoverflow.com/questions/8240674/using-geos-via-pythons-ctypes ).

I am VERY sorry for cross-posting but I really do need a solution. VERY many thanks.

Answer 1

I used the 2to3 tool (see geographika's answer) to convert Shapely-1.2.13 and then manually had to change 3 lines of code to handle Python 3.x's explicit handling of string and byte encoding. I have put this as a separate answer just for clarity and in case anybody else needs to do the same here are the lines I manually changed:

1. \geos.py (line 97)

   old: v = func().encode('ascii').split('-')[2]
   error: TypeError: Type str doesn't support the buffer API
   new: v = func().decode("utf-8").split('-')[2]

2. \geometry\base.py (line 60)

   old: EMPTY = wkb.deserialize('010700000000000000'.decode('hex'))
   error: AttributeError: 'str' object has no attribute 'decode'
   new: EMPTY = wkb.deserialize(bytes.fromhex('010700000000000000'))
   

3. \geometry\proxy.py (line 7)

   same as above

I then ran the same test as the developer and got the same result. With the above knowledge I could have done the conversion in less than 5 minutes... I just had to read up on how Python 3.x handles strings and bytes differently :).

Granted, this is not an exhaustive test but... so far so good and I hope this helps other Shapely people!

EDIT Running the doctests gives me a few more errors but they all seem to be of the same type.. so my fix will take a bit longer and I'll not document everything here but we're on the right track now!

• MappaGnosis: I have managed to eliminate all but four errors from the doctests (out of 13 or 16 if you count the three I caught prior to running doctest above). These four I think can be ignored. Three of them appear to due to differences in decimal precision between Python 3 and 2 (e.g. highly edited example: Expected (0.10000000000000001, 0.20000000000000001); Got (0.1, 0.2)). The last remaining "error" is from raising a topology exception. I get the correct result but a lot of extra verbage. Perhaps Python 3 is more verbose - I don't know but as long as the exception is correctly raised, I'm happy.
• geographika: Good stuff. Thanks for posting your approach back here.
• MappaGnosis: I have been using my ported version of Shapely for a few days now and the only issue I have encountered is the LineMerge process crashing when I have a special case of bad geometry (viz. a line of two identical points) - so, the port seems to be good... I just need to fix my geometry now :)

Answer 2

In regards to Shapely and Python 3.2:

We're not yet testing or deploying Shapely on Python 3 + Windows. There's probably some non-conforming code in the shapely.geometry.geo module that prevents its import.

June, 2011 - http://old.nabble.com/Shapely-package-installation-td31879542.html

So you could try getting the Shapely source from GitHub and using the automated 2to3 tool to translate the code. Then rebuild for Windows (all the required files are in GitHub).

To see how to use ctypes and GEOS have a look at the Shapely source:

https://github.com/sgillies/shapely/blob/master/shapely/geos.py

However it may be time better spent trying to fix some Python inconsistencies in the Shapely source than trying to create another wrapper.

• GeospatialPython.com: Are you using a large portion of the API or a smaller subset?
• MappaGnosis: Thanks for the comments. I'm only really wanting one function at the moment but it is quite crucial. I need the LineMerge tool for multi-linestrings. I can't just append geometry manually because the order of vertices MUST correspond to the line without jumping back and forth.
• MappaGnosis: @geographika: I looked at the Shapely code because I only want one function from GEOS. I understand how to connect to GEOS but I don't understand how to pass data in and out. I'm sorry for being a bit slow but I don't see the code in that module where this is happening. I'm very new to CTypes so a worked example of the LineMerge function would help me understand.

Answer 3

To follow up on MappaGnosis's answer from December 2011, here's how I did it with Shapely 1.2.17:

Download the source, extract, and cd into the source directory.

Use 2to3 to update all the Python code:

$ find ./ -type f -name '*.py' -exec 2to3 -w -n {} \;

Replace all the 'rb' with 'r' and 'wb' with 'w' in setup.py. This shouldn't have any consequences on a Unix-style platform, but I'm not sure about Windows.

Fix line 97 of shapely/geos.py:

- v = func()
+ v = func().decode("utf-8")

Fix line 87 of shapely/geometry/base.py:

- EMPTY = deserialize_wkb('010700000000000000'.decode('hex'))
+ EMPTY = deserialize_wkb(bytes.fromhex('010700000000000000'))

Install:

$ python setup.py install

I haven't done any testing of this hacked version to see how it works, but it does install, which is a start!

• MappaGnosis: Did you run the doctests? Have a look at my follow-up comment and you'll see that in the old version you needed to do a little more work to have Shapely fully functioning. I'm curious as to whether this is still necessary.
• cswingle: @MappaGnosis - Err, good question. I just ran the doctests and the results were littered with ERRORs and FAILs. Many places where bytes of integer address expected instead of str instance and some issues related to failed comparisons between strings and byte ('POINT' vs. b'POINT'). The final total was 15 failures and 8 errors out of 86 tests.
• MappaGnosis: It sounds a similar to what I got. When I fixed the errors/failures, Shapely has run perfectly and without problems ever since. So, it is very doable and worth the tiny amount of pain.