Skip to main content
Turned question in answer into statement
Source Link
PolyGeo
  • 65.4k
  • 29
  • 114
  • 345

It'sIt is not clear to me: are 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.

It's not clear to me: are you 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.

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.

Note about lib versions.
Source Link
Rob Skelly
  • 3.1k
  • 16
  • 14

(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++).

(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.)

(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++).

Source Link
Rob Skelly
  • 3.1k
  • 16
  • 14

It's not clear to me: are you 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.)

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));