PostGIS geometry to GEOS geometry without lwgeom dependency?

Question

I've been looking for ways to convert PostGIS to GEOS geometries without linking to or dependency on lwgeom/liblwgeom.

The closest thing I can find is a mention of the pointcloud approach here:

For an alternate example of an extension with a lighter dependency on PostGIS, check out pgpointcloud, which has it's own structure for spatial data (a point patch) and exchanges data with PostGIS via well-known-binary. This removes the liblwgeom dependency, which means it's possible to compile and use pgpointcloud without PostGIS installed, which is not entirely uncommon.

But I do not want to reinvent geometry types, and instead want to use what GEOS already have (, which I understand, corresponds to be subset of PostGIS/lwgeom types). So I am thinking about the path

PostGIS geometry  ==> (step1. via ST_AsEWKB or its like) EWKB 
==> (step2. via GEOSGeomFromWKB_buf) GEOS geometry ==> call geos functions...

But I am not familiar with the code structures of the packages.

What I know so far are:

1. ST_AsEWKB for step1, but I guess I'll need a C/C++ equivalent.

2. Some examples reading WKB in step2.

   /* WKB input */
   g2 = GEOSGeomFromWKB_buf(uptr, size); free(uptr);
   

Can anyone with the expertise give some pointers to completing the loop? (e.g. about specific pointcloud functions, better work flows, etc.)

(My target environment is a PostgreSQL extension (PostgreSQL 11) with GEOS 3.7 and PostGIS 2.5 as of Jan 2019.)

Other references:

How to use GEOS/C++ to efficiently find all point pairs closer than a threshold?

Answer 1

It is not clear to me whether you are looking specifically for information on loading point clouds, or geometries in general.

GEOS may not be the ideal tool for manipulating point clouds, but you can certainly transport a subset as point geometries. My example below uses a polygon, but the general strategy is geometry-agnostic.

(I hadn't tried this before, but it seemed like a good opportunity to learn it -- it's possible that I've missed some shortcuts or clean-up steps. Also, I used the versions of GEOS (3.5) and PostGIS (2.5) on my system, but it should be similar -- the GEOS API changes I've noticed are in C++).

The first gotcha is that you want to load the EWKB string from PostGIS and encode it as HEX so it can be consumed by the GEOS reader. You have to specifically encode the binary output using the PostgreSQL encode function (see below). GEOS does have a binary reader, but I didn't get it to work.

Here's a brief example in C that loads a square polygon and prints out the points. A working example is on github.

// Connect to PG and execute the query.
PGconn* conn = PQconnectdb("host=localhost dbname=test user=rob");
PGresult* res = PQexec(conn, "SELECT encode(ST_AsEWKB(geom), 'hex') FROM test LIMIT 1");

// Grab the geometry and its length.
const unsigned char* geom = PQgetvalue(res, 0, 0);
int len = strlen(geom);

// End the PG connection.
PQfinish(conn);

// Initialize GEOS. Don't forget this. You'll be scratching your head.
initGEOS(&notice, &notice);

// Create a reader and use it to read the geometry.
GEOSWKBReader* reader = GEOSWKBReader_create();
GEOSGeometry* g = GEOSWKBReader_readHEX(reader, geom, len);
GEOSWKBReader_destroy(reader);

// Retrieve the exterior ring and its coordinates.
const GEOSGeometry* r = GEOSGetExteriorRing(g);
const GEOSCoordSequence* s = GEOSGeom_getCoordSeq(r);

// Iterate over the points in the ring.
int n = GEOSGeomGetNumPoints(r);
double x, y;

for(int i = 0; i < n; ++i) {
    GEOSCoordSeq_getX(s, i, &x);
    GEOSCoordSeq_getY(s, i, &y);
    printf("Coord %d: %f %f\n", i, x, y);
}

// Cleanup.
GEOSGeom_destroy(g);
finishGEOS();

The output looks like this:

Number of points: 5
Coord 0: 1000.000000 1000.000000
Coord 1: 2000.000000 1000.000000
Coord 2: 2000.000000 2000.000000
Coord 3: 1000.000000 2000.000000
Coord 4: 1000.000000 1000.000000

The original geometry was created thusly:

insert into test (geom) values (st_geomfromtext('polygon((1000 1000, 2000 1000, 2000 2000, 1000 2000, 1000 1000))', 26910));