0

I am trying to speed the spatial distance calculation between polygons and points/polygons by calling the GEOS library directly. However I couldn't find any help how to call this function correctly. Can anyone please point me to the location where I can find the reference for this function or point out where I have done incorrectly?

working example:

from shapely.geos import lgeos

points_geom = np.array([x._geom for x in points])
polygons_geom = np.array([x._geom for x in polygons])  
lgeos._lgeos.GEOSContains_r(lgeos.geos_handle,polygons_geom[0],points_geom[0])

Not working:

lgeos._lgeos.GEOSDistance_r(lgeos.geos_handle,polygons_geom[0],points_geom[0])


TypeError                                 Traceback (most recent call last)
<ipython-input-138-392cb700cfbc> in <module>()
----> 1 lgeos._lgeos.GEOSDistance_r(lgeos.geos_handle,polygons_geom[0],points_geom[0])

TypeError: this function takes at least 4 arguments (3 given)
ngoldbaum
  • 5,430
  • 3
  • 28
  • 35
Karen Chen
  • 65
  • 1
  • 7

1 Answers1

1

GEOSDistance_r takes 4 arguments, you're only passing three:

extern int GEOS_DLL GEOSDistance_r(GEOSContextHandle_t handle,
                                   const GEOSGeometry* g1,
                                   const GEOSGeometry* g2, double *dist);

(from https://github.com/OSGeo/geos/blob/5a730fc50dab2610a9e6c037b521accc66b7777b/capi/geos_c.h.in#L1100)

You're using shapely's private interface to GEOS, which it looks like uses ctypes, so you'll need to use the ctypes invocation to pass a double by reference:

import ctypes
dist = ctypes.c_double()
lgeos._lgeos.GEOSContains_r(
    lgeos.geos_handle, polygons_geom[0], points_geom[0], ctypes.byref(dist))
dist = dist.value
ngoldbaum
  • 5,430
  • 3
  • 28
  • 35