Seeking ST_Subdivide implementation in Python?

Question

I'd like to utilize PostGIS' ST_Subdivide operation via Python. I would like to avoid interfacing directly with PostGIS, have a database up, etc. Instead, I would like to pass a geometry object and have the function evaluated in memory.

I can see that ST_Subdivide wraps over the function lwgeom_subdivide which itself wraps a recursive method defined here.

I am looking for either:

• An implementation of this algorithm in Python
• An isolated version of this method that has Python bindings... or some resources that might help me go about doing this
• A pattern for potentially accessing this function from Python (likely not a reasonable request)

Answer 1

You can use a dark art of Python called ctypes to interface with liblwgeom used by Spatialite and other non-PostgreSQL software. On my Windows computer, I have lwgeom.dll from OSGeo4W, but you can also look for liblwgeom.so or other names on other systems.

import ctypes

lib_path = 'liblwgeom.so'  # Linux, etc
lib_path = 'lwgeom.dll'  # Windows
lwgeom = ctypes.cdll.LoadLibrary(lib_path)

## See liblwgeom.h for public function definitions
# extern LWGEOM* lwgeom_from_wkt(const char *wkt, const char check);
lwgeom.lwgeom_from_wkt.argtypes = (ctypes.c_char_p, ctypes.c_byte)
lwgeom.lwgeom_from_wkt.restype = ctypes.c_void_p

# extern char *lwgeom_to_ewkt(const LWGEOM *lwgeom);
lwgeom.lwgeom_to_ewkt.argtypes = (ctypes.c_void_p,)
lwgeom.lwgeom_to_ewkt.restype = ctypes.c_char_p

# LWCOLLECTION *lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices);
lwgeom.lwgeom_subdivide.argtypes = (ctypes.c_void_p, ctypes.c_uint32)
lwgeom.lwgeom_subdivide.restype = ctypes.c_void_p

# Example
in_wkt = '''\
POLYGON((132 10,119 23,85 35,68 29,66 28,49 42,32 56,22 64,32 
110,40 119,36 150,
57 158,75 171,92 182,114 184,132 186,146 178,176 184,179 162,184 141,190 122,
190 100,185 79,186 56,186 52,178 34,168 18,147 13,132 10))'''
in_geom = lwgeom.lwgeom_from_wkt(in_wkt, 0)
print(lwgeom.lwgeom_to_ewkt(in_geom))  # check

# Try to process
out_geom = lwgeom.lwgeom_subdivide(in_geom, 8)
out_wkt = lwgeom.lwgeom_to_ewkt(out_geom)
print(out_wkt)

Should print a result:

GEOMETRYCOLLECTION(POLYGON((34.4285714286 54,64 54,64 29.6470588235,49 42,34.4285714286 54)),POLYGON((93.5 32,106 32,106 27.5882352941,93.5 32)),POLYGON((64 29.6470588235,64 32,76.5 32,68 29,66 28,64 29.6470588235)),POLYGON((64 32,64 54,85 54,85 35,76.5 32,64 32)),POLYGON((85 35,85 54,106 54,106 32,93.5 32,85 35)),POLYGON((22 64,29.3913043478 98,106 98,106 54,34.4285714286 54,32 56,22 64)),POLYGON((106 27.5882352941,106 32,127 32,127 15,119 23,106 27.5882352941)),POLYGON((127 15,127 32,148 32,148 13.2380952381,147 13,132 10,127 15)),POLYGON((106 32,106 54,148 54,148 32,106 32)),POLYGON((148 13.2380952381,148 54,186 54,186 52,178 34,168 18,148 13.2380952381)),POLYGON((106 54,106 98,189.523809524 98,185 79,186 56,186 54,106 54)),POLYGON((29.3913043478 98,32 110,40 119,37.0322580645 142,106 142,106 98,29.3913043478 98)),POLYGON((36 150,57 158,64 163.055555556,64 142,37.0322580645 142,36 150)),POLYGON((64 142,64 163.055555556,75 171,92 182,106 183.272727273,106 142,64 142)),POLYGON((106 98,106 142,148 142,148 98,106 98)),POLYGON((148 98,148 120,190 120,190 100,189.523809524 98,148 98)),POLYGON((148 120,148 142,183.761904762 142,184 141,190 122,190 120,148 120)),POLYGON((106 142,106 164,148 164,148 142,106 142)),POLYGON((106 164,106 183.272727273,114 184,127 185.444444444,127 164,106 164)),POLYGON((127 164,127 185.444444444,132 186,146 178,148 178.4,148 164,127 164)),POLYGON((148 142,148 178.4,176 184,179 162,183.761904762 142,148 142)))

---

Errors like:

The GEOS version this postgis binary was compiled against (34) doesn't support 'GEOSClipByRect' function (3.5.0+ required)

show up if you compiled liblwgeom with an older GEOS, which can only be fixed by recompiling liblwgeom with a new GEOS. As of writing, this example does not work with OSGeo4W...