crossbind
GitHub
← All libraries

GEOS

Geometry predicates and spatial operations.

Install

Runs in browsers; Node.js and edge runtimes use the same package. Built for wasm32, single-threaded and multi-threaded.

shell
npm install @crossbind/port-geos-wasm@beta
crossbind.config.js
import geosWasm from '@crossbind/port-geos-wasm/crossbind.config.js';
 
export default {
dependencies: [geosWasm],
paths: { config: import.meta.url },
};
main.js
import { initNative, GEOSversion } from '@crossbind/port-geos/geos_c.h';
 
await initNative();

crossbind itself comes with your bundler plugin; Browser has Vite, Webpack, Rspack and Rollup.

The geometry engine behind PostGIS, in your browser

GEOS is the C++ library PostGIS, QGIS, GDAL and Shapely call for overlays, validity checks and simplification. These apps run it, compiled by crossbind, on shapes you type: fifteen operations with the DE-9IM relation of two shapes, a diagnosis and two repairs for broken polygons, and a map simplified without opening gaps between regions. GEOS is LGPL-2.1; every app links its source. The first run downloads 1.7 MB of WebAssembly once; every app on this page shares it, and nothing is uploaded.

APP 01

Fifteen geometry operations on shapes you can edit

Pick two shapes or type your own as WKT, then run an overlay, a buffer, a hull, a triangulation or a repair. The result is drawn over the inputs with its measurements, next to the DE-9IM matrix of how the two shapes relate. This is the engine PostGIS, QGIS, GDAL and Shapely call, built from the same C++.

Runs GEOS, which is LGPL-2.1: source · license · build recipe
Run an operation to draw its result over the two shapes.
SHOW THE CODE
src/native/geometry_lab.h
// src/native/geometry_lab.h (excerpt)
GEOSGeometry* apply(const std::string& operation, const GEOSGeometry* a,
const GEOSGeometry* b, double parameter) {
GEOSContextHandle_t context = geos.context;
if (operation == "intersection") return GEOSIntersection_r(context, a, b);
if (operation == "union") return GEOSUnion_r(context, a, b);
if (operation == "buffer") return GEOSBuffer_r(context, a, parameter, 8);
if (operation == "makeValid") return GEOSMakeValid_r(context, a);
if (operation == "maximumInscribedCircle") return GEOSMaximumInscribedCircle_r(context, a, parameter);
// ...ten more, then the ones that work on A and B together
const auto both = together(a, b);
if (operation == "voronoi") return GEOSVoronoiDiagram_r(context, both.get(), nullptr, 0, 0);
if (operation == "delaunay") return GEOSDelaunayTriangulation_r(context, both.get(), 0, 0);
throw std::invalid_argument("unknown operation " + operation);
}
main.js
const m = await initNative();
const lab = await new m.GeometryLab();
const park = 'POLYGON ((0 0, 60 0, 60 40, 0 40, 0 0), (35 15, 45 15, 45 25, 35 25, 35 15))';
const river = 'LINESTRING (-10 30, 20 30, 40 20, 70 20)';
 
const inside = JSON.parse(await lab.run('intersection', park, river, 0));
// inside.wkt: MULTILINESTRING ((0 30, 20 30, 35 22.5), (45 20, 60 20))
// inside.length: 51.77, the river's length in the park, pond excluded
const relation = JSON.parse(await lab.relate(park, river));
// relation.matrix: '1F20F1102', relation.crosses: true
APP 02

Find out why a polygon is invalid, then repair it two ways

Invalid polygons break overlays and spatial joins. GEOS names the problem and marks where it is, then repairs the shape with either of its two methods: one rebuilds the shape from all of its edges, the other keeps shells as shells and holes as holes. They agree on some shapes and not on others. jsts 2.12.1, the JavaScript port of the same Java library, has no GeometryFixer, the class behind the second method.

Runs GEOS, which is LGPL-2.1: source · license · build recipe
Pick a broken shape, or paste one, to see what is wrong and both repairs.
SHOW THE CODE
src/native/validity_doctor.h
// src/native/validity_doctor.h (excerpt)
const char valid = GEOSisValidDetail_r(geos.context, geometry.get(), flags, &reason, &location);
 
GEOSMakeValidParams* params = GEOSMakeValidParams_create_r(geos.context);
GEOSMakeValidParams_setMethod_r(geos.context, params,
method == "structure" ? GEOS_MAKE_VALID_STRUCTURE : GEOS_MAKE_VALID_LINEWORK);
GEOSMakeValidParams_setKeepCollapsed_r(geos.context, params, keepCollapsed ? 1 : 0);
GEOSGeometry* repaired = GEOSMakeValidWithParams_r(geos.context, geometry.get(), params);
GEOSMakeValidParams_destroy_r(geos.context, params);
main.js
const m = await initNative();
const doctor = await new m.ValidityDoctor();
const leaky = 'POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (5 5, 15 5, 15 15, 5 15, 5 5))';
 
const found = JSON.parse(await doctor.check(leaky, false));
// found.reason: 'Self-intersection', found.location: [5, 10]
const linework = JSON.parse(await doctor.repair(leaky, 'linework', false));
const structure = JSON.parse(await doctor.repair(leaky, 'structure', false));
// linework.area: 150 in 2 parts; structure.area: 75, the shell minus its hole
APP 03

Simplify a map without opening gaps between regions

Simplify each region of a map on its own and neighbours stop agreeing on their borders: slivers of gap and overlap open along every edge they share. GEOS can simplify the map as one coverage instead, so each border is simplified once and stays shared. jsts 2.12.1, the JavaScript port of the same Java library, has no coverage package.

The map is 36 regions in a 100 × 100 square, generated in the module. Every inner border is computed once and shared by its two regions, as in a real administrative map. Both methods get the same tolerance: GEOSTopologyPreserveSimplify for each region alone, GEOSCoverageSimplifyVW for the coverage.

Runs GEOS, which is LGPL-2.1: source · license · build recipe
Simplify the map to compare the two methods, border by border.
SHOW THE CODE
src/native/coverage_lab.h
// src/native/coverage_lab.h (excerpt)
// One region at a time: each shared border is simplified twice, once per neighbour.
for (int index = 0; index < geos.parts(map.get()); index += 1) {
each.push_back(geos.own(GEOSTopologyPreserveSimplify_r(
geos.context, GEOSGetGeometryN_r(geos.context, map.get(), index), tolerance)));
}
// As one coverage: each shared border is simplified once, for both neighbours.
const auto together = geos.own(GEOSCoverageSimplifyVW_r(geos.context, map.get(), tolerance, 0));
 
// The audit: what the regions no longer cover, and whether GEOS still sees a coverage.
const auto merged = geos.own(GEOSUnaryUnion_r(geos.context, regions));
const auto gaps = geos.own(GEOSDifference_r(geos.context, square.get(), merged.get()));
const int result = GEOSCoverageIsValid_r(geos.context, regions, 0, nullptr);
main.js
const m = await initNative();
const lab = await new m.CoverageLab(6); // 36 regions sharing wiggly borders
const result = JSON.parse(await lab.simplify(2));
 
// result.separately: 430 vertices, gaps 356.917, overlaps 186.433, 155 slivers
// result.coverage: 506 vertices, gaps 0, overlaps 0, coverageValid true

Examples

Each one is a small C++ header and the JavaScript that calls it, with what it printed when the site was built. The JavaScript only tab does the same with GEOS's own headers and no C++ file. The last is a whole C++ program, built for WASI.

Intersect, unite and subtract two polygons

Overlay is the most used part of GEOS: GEOSIntersection_r, GEOSUnion_r and GEOSDifference_r on shapes read from WKT, with GEOSArea_r to measure the result.

src/native/overlay.h
#pragma once
 
#include <geos_c.h>
 
#include <memory>
#include <stdexcept>
#include <string>
 
// Overlay of two geometries written as WKT. Every result is normalised, so the same shape always
// prints the same WKT whatever order its vertices came in.
class Overlay {
public:
static std::string version() { return GEOSversion(); }
 
static std::string intersection(const std::string& a, const std::string& b) { return apply(GEOSIntersection_r, a, b); }
static std::string unite(const std::string& a, const std::string& b) { return apply(GEOSUnion_r, a, b); }
static std::string difference(const std::string& a, const std::string& b) { return apply(GEOSDifference_r, a, b); }
 
static double area(const std::string& wkt) {
Session geos;
double value = 0;
if (!GEOSArea_r(geos.context, geos.read(wkt).get(), &value)) geos.fail();
return value;
}
 
private:
using Operation = GEOSGeometry* (*)(GEOSContextHandle_t, const GEOSGeometry*, const GEOSGeometry*);
 
static std::string apply(Operation operation, const std::string& a, const std::string& b) {
Session geos;
const auto left = geos.read(a);
const auto right = geos.read(b);
return geos.result(operation(geos.context, left.get(), right.get()));
}
 
// GEOS's reentrant C API: each call gets its own context, which also holds the last error message.
struct Session {
struct Destroy {
GEOSContextHandle_t context;
void operator()(GEOSGeometry* geometry) const { GEOSGeom_destroy_r(context, geometry); }
};
using Geometry = std::unique_ptr<GEOSGeometry, Destroy>;
 
GEOSContextHandle_t context = GEOS_init_r();
std::string error;
 
Session() { GEOSContext_setErrorMessageHandler_r(context, remember, &error); }
~Session() { GEOS_finish_r(context); }
Session(const Session&) = delete;
Session& operator=(const Session&) = delete;
 
Geometry own(GEOSGeometry* geometry) {
if (!geometry) fail();
return Geometry(geometry, Destroy{context});
}
 
Geometry read(const std::string& wkt) {
GEOSWKTReader* reader = GEOSWKTReader_create_r(context);
GEOSGeometry* geometry = GEOSWKTReader_read_r(context, reader, wkt.c_str());
GEOSWKTReader_destroy_r(context, reader);
return own(geometry);
}
 
std::string write(const GEOSGeometry* geometry) {
GEOSWKTWriter* writer = GEOSWKTWriter_create_r(context);
char* text = GEOSWKTWriter_write_r(context, writer, geometry);
GEOSWKTWriter_destroy_r(context, writer);
if (!text) fail();
const std::string wkt = text;
GEOSFree_r(context, text);
return wkt;
}
 
// Takes ownership of a result and writes it normalised.
std::string result(GEOSGeometry* geometry) {
const Geometry owned = own(geometry);
if (GEOSNormalize_r(context, owned.get()) != 0) fail();
return write(owned.get());
}
 
[[noreturn]] void fail() const { throw std::runtime_error(error.empty() ? "GEOS operation failed" : error); }
 
static void remember(const char* message, void* error) { *static_cast<std::string*>(error) = message; }
};
};
main.js
import { initNative, Overlay } from './native/overlay.h';
 
await initNative();
const parcel = 'POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))';
const floodZone = 'POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))';
const flooded = await Overlay.intersection(parcel, floodZone);
console.log((await Overlay.version()).split('-CAPI-')[0]);
console.log(flooded, await Overlay.area(flooded));
console.log(await Overlay.area(await Overlay.unite(parcel, floodZone)), await Overlay.area(await Overlay.difference(parcel, floodZone)));
PRINTSfirst run downloads 1.7 MB
3.15.0
POLYGON ((5 5, 5 10, 10 10, 10 5, 5 5)) 25
175 75

Test points against a prepared polygon

Point in polygon is the most common spatial question. GEOSPrepare_r indexes the polygon once; GEOSPreparedContainsXY_r and GEOSPreparedIntersectsXY_r then answer per point and differ only on the boundary, and GEOSPreparedRelate_r gives the full DE-9IM matrix.

src/native/zone.h
#pragma once
 
#include <geos_c.h>
 
#include <memory>
#include <stdexcept>
#include <string>
 
// A polygon prepared once and then asked about many points: GEOSPrepare_r indexes its edges on the
// first question, so the later ones do not walk every edge again.
class Zone {
public:
explicit Zone(const std::string& wkt) : shape(geos.read(wkt)), prepared(GEOSPrepare_r(geos.context, shape.get())) {
if (!prepared) geos.fail();
}
~Zone() { GEOSPreparedGeom_destroy_r(geos.context, prepared); }
Zone(const Zone&) = delete;
Zone& operator=(const Zone&) = delete;
 
// Strictly inside: a point on the boundary is not contained.
bool contains(double x, double y) { return answer(GEOSPreparedContainsXY_r(geos.context, prepared, x, y)); }
 
// Inside or on the boundary.
bool intersects(double x, double y) { return answer(GEOSPreparedIntersectsXY_r(geos.context, prepared, x, y)); }
 
// The DE-9IM matrix: the dimension where the interior, boundary and exterior of this shape meet
// those of the other, row by row.
std::string relate(const std::string& wkt) {
char* matrix = GEOSPreparedRelate_r(geos.context, prepared, geos.read(wkt).get());
if (!matrix) geos.fail();
const std::string text = matrix;
GEOSFree_r(geos.context, matrix);
return text;
}
 
private:
bool answer(char result) {
if (result == 2) geos.fail();
return result == 1;
}
 
// GEOS's reentrant C API: each call gets its own context, which also holds the last error message.
struct Session {
struct Destroy {
GEOSContextHandle_t context;
void operator()(GEOSGeometry* geometry) const { GEOSGeom_destroy_r(context, geometry); }
};
using Geometry = std::unique_ptr<GEOSGeometry, Destroy>;
 
GEOSContextHandle_t context = GEOS_init_r();
std::string error;
 
Session() { GEOSContext_setErrorMessageHandler_r(context, remember, &error); }
~Session() { GEOS_finish_r(context); }
Session(const Session&) = delete;
Session& operator=(const Session&) = delete;
 
Geometry own(GEOSGeometry* geometry) {
if (!geometry) fail();
return Geometry(geometry, Destroy{context});
}
 
Geometry read(const std::string& wkt) {
GEOSWKTReader* reader = GEOSWKTReader_create_r(context);
GEOSGeometry* geometry = GEOSWKTReader_read_r(context, reader, wkt.c_str());
GEOSWKTReader_destroy_r(context, reader);
return own(geometry);
}
 
std::string write(const GEOSGeometry* geometry) {
GEOSWKTWriter* writer = GEOSWKTWriter_create_r(context);
char* text = GEOSWKTWriter_write_r(context, writer, geometry);
GEOSWKTWriter_destroy_r(context, writer);
if (!text) fail();
const std::string wkt = text;
GEOSFree_r(context, text);
return wkt;
}
 
// Takes ownership of a result and writes it normalised.
std::string result(GEOSGeometry* geometry) {
const Geometry owned = own(geometry);
if (GEOSNormalize_r(context, owned.get()) != 0) fail();
return write(owned.get());
}
 
[[noreturn]] void fail() const { throw std::runtime_error(error.empty() ? "GEOS operation failed" : error); }
 
static void remember(const char* message, void* error) { *static_cast<std::string*>(error) = message; }
};
 
Session geos;
Session::Geometry shape;
const GEOSPreparedGeometry* prepared;
};
main.js
import { initNative, Zone } from './native/zone.h';
 
await initNative();
const zone = await new Zone('POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (4 4, 6 4, 6 6, 4 6, 4 4))');
const points = [[2, 2], [5, 5], [10, 5], [12, 5]]; // inside, in the hole, on the edge, outside
for (const [x, y] of points) {
console.log(`${x} ${y}: contains ${await zone.contains(x, y)}, intersects ${await zone.intersects(x, y)}`);
}
console.log(await zone.relate('POLYGON ((5 5, 15 5, 15 15, 5 15, 5 5))'));
PRINTSfirst run downloads 1.7 MB
2 2: contains true, intersects true
5 5: contains false, intersects false
10 5: contains false, intersects true
12 5: contains false, intersects false
212101212

Measure area, length and distance

GEOSArea_r, GEOSLength_r and GEOSDistance_r in the units of the coordinates, GEOSNearestPoints_r for where two shapes come closest, and GEOSGetCentroid_r. GEOS works in the plane, so project longitude and latitude before measuring.

src/native/measure.h
#pragma once
 
#include <geos_c.h>
 
#include <memory>
#include <stdexcept>
#include <string>
 
// Measurements in the plane, in the units of the coordinates: metres for most projected data.
// GEOS does not measure on the ellipsoid, so longitude and latitude need projecting first.
class Measure {
public:
static double area(const std::string& wkt) { return measure(GEOSArea_r, wkt); }
 
// A line's length or a polygon's perimeter. Not named length: every JavaScript function, a
// bound class included, already has a length property, and a static method cannot replace it.
static double lengthOf(const std::string& wkt) { return measure(GEOSLength_r, wkt); }
 
static double distance(const std::string& a, const std::string& b) {
Session geos;
const auto left = geos.read(a);
const auto right = geos.read(b);
double value = 0;
if (!GEOSDistance_r(geos.context, left.get(), right.get(), &value)) geos.fail();
return value;
}
 
// The shortest line between the two shapes, from its end on a to its end on b.
static std::string nearestPoints(const std::string& a, const std::string& b) {
Session geos;
const auto left = geos.read(a);
const auto right = geos.read(b);
GEOSCoordSequence* ends = GEOSNearestPoints_r(geos.context, left.get(), right.get());
if (!ends) geos.fail();
return geos.write(geos.own(GEOSGeom_createLineString_r(geos.context, ends)).get());
}
 
static std::string centroid(const std::string& wkt) {
Session geos;
return geos.result(GEOSGetCentroid_r(geos.context, geos.read(wkt).get()));
}
 
private:
using Metric = int (*)(GEOSContextHandle_t, const GEOSGeometry*, double*);
 
static double measure(Metric metric, const std::string& wkt) {
Session geos;
double value = 0;
if (!metric(geos.context, geos.read(wkt).get(), &value)) geos.fail();
return value;
}
 
// GEOS's reentrant C API: each call gets its own context, which also holds the last error message.
struct Session {
struct Destroy {
GEOSContextHandle_t context;
void operator()(GEOSGeometry* geometry) const { GEOSGeom_destroy_r(context, geometry); }
};
using Geometry = std::unique_ptr<GEOSGeometry, Destroy>;
 
GEOSContextHandle_t context = GEOS_init_r();
std::string error;
 
Session() { GEOSContext_setErrorMessageHandler_r(context, remember, &error); }
~Session() { GEOS_finish_r(context); }
Session(const Session&) = delete;
Session& operator=(const Session&) = delete;
 
Geometry own(GEOSGeometry* geometry) {
if (!geometry) fail();
return Geometry(geometry, Destroy{context});
}
 
Geometry read(const std::string& wkt) {
GEOSWKTReader* reader = GEOSWKTReader_create_r(context);
GEOSGeometry* geometry = GEOSWKTReader_read_r(context, reader, wkt.c_str());
GEOSWKTReader_destroy_r(context, reader);
return own(geometry);
}
 
std::string write(const GEOSGeometry* geometry) {
GEOSWKTWriter* writer = GEOSWKTWriter_create_r(context);
char* text = GEOSWKTWriter_write_r(context, writer, geometry);
GEOSWKTWriter_destroy_r(context, writer);
if (!text) fail();
const std::string wkt = text;
GEOSFree_r(context, text);
return wkt;
}
 
// Takes ownership of a result and writes it normalised.
std::string result(GEOSGeometry* geometry) {
const Geometry owned = own(geometry);
if (GEOSNormalize_r(context, owned.get()) != 0) fail();
return write(owned.get());
}
 
[[noreturn]] void fail() const { throw std::runtime_error(error.empty() ? "GEOS operation failed" : error); }
 
static void remember(const char* message, void* error) { *static_cast<std::string*>(error) = message; }
};
};
main.js
import { initNative, Measure } from './native/measure.h';
 
await initNative();
const field = 'POLYGON ((0 0, 40 0, 40 30, 0 30, 0 0))';
const well = 'POINT (52 39)';
const track = 'LINESTRING (0 0, 30 40, 30 50)';
console.log(await Measure.area(field), await Measure.lengthOf(field), await Measure.lengthOf(track));
console.log(await Measure.distance(field, well), await Measure.nearestPoints(field, well));
console.log(await Measure.centroid(field));
PRINTSfirst run downloads 1.7 MB
1200 140 60
15 LINESTRING (40 30, 52 39)
POINT (20 15)

Buffer points, lines and polygons

GEOSBuffer_r turns a point into a circle of straight segments, GEOSBufferWithStyle_r sets the caps and joins of a line's corridor, a negative distance shrinks a polygon, and GEOSOffsetCurve_r draws a parallel line.

src/native/buffer.h
#pragma once
 
#include <geos_c.h>
 
#include <memory>
#include <stdexcept>
#include <string>
 
// Buffers grow a shape by a distance, or shrink it with a negative one; an offset curve runs parallel
// to a line. Curves become straight segments: quadrantSegments of them for each quarter circle.
class BufferOp {
public:
static std::string around(const std::string& wkt, double distance, int quadrantSegments) {
Session geos;
return geos.result(GEOSBuffer_r(geos.context, geos.read(wkt).get(), distance, quadrantSegments));
}
 
// cap: "round", "flat" or "square"; join: "round", "mitre" or "bevel".
static std::string withStyle(const std::string& wkt, double distance, const std::string& cap, const std::string& join) {
Session geos;
const auto geometry = geos.read(wkt);
return geos.result(GEOSBufferWithStyle_r(geos.context, geometry.get(), distance, 8, capStyle(cap), joinStyle(join), 5.0));
}
 
// A positive distance offsets to the left of the line's direction, a negative one to the right.
static std::string offsetCurve(const std::string& wkt, double distance) {
Session geos;
return geos.result(GEOSOffsetCurve_r(geos.context, geos.read(wkt).get(), distance, 8, GEOSBUF_JOIN_ROUND, 5.0));
}
 
private:
static int capStyle(const std::string& name) {
if (name == "round") return GEOSBUF_CAP_ROUND;
if (name == "flat") return GEOSBUF_CAP_FLAT;
if (name == "square") return GEOSBUF_CAP_SQUARE;
throw std::invalid_argument("cap must be round, flat or square");
}
 
static int joinStyle(const std::string& name) {
if (name == "round") return GEOSBUF_JOIN_ROUND;
if (name == "mitre") return GEOSBUF_JOIN_MITRE;
if (name == "bevel") return GEOSBUF_JOIN_BEVEL;
throw std::invalid_argument("join must be round, mitre or bevel");
}
 
// GEOS's reentrant C API: each call gets its own context, which also holds the last error message.
struct Session {
struct Destroy {
GEOSContextHandle_t context;
void operator()(GEOSGeometry* geometry) const { GEOSGeom_destroy_r(context, geometry); }
};
using Geometry = std::unique_ptr<GEOSGeometry, Destroy>;
 
GEOSContextHandle_t context = GEOS_init_r();
std::string error;
 
Session() { GEOSContext_setErrorMessageHandler_r(context, remember, &error); }
~Session() { GEOS_finish_r(context); }
Session(const Session&) = delete;
Session& operator=(const Session&) = delete;
 
Geometry own(GEOSGeometry* geometry) {
if (!geometry) fail();
return Geometry(geometry, Destroy{context});
}
 
Geometry read(const std::string& wkt) {
GEOSWKTReader* reader = GEOSWKTReader_create_r(context);
GEOSGeometry* geometry = GEOSWKTReader_read_r(context, reader, wkt.c_str());
GEOSWKTReader_destroy_r(context, reader);
return own(geometry);
}
 
std::string write(const GEOSGeometry* geometry) {
GEOSWKTWriter* writer = GEOSWKTWriter_create_r(context);
char* text = GEOSWKTWriter_write_r(context, writer, geometry);
GEOSWKTWriter_destroy_r(context, writer);
if (!text) fail();
const std::string wkt = text;
GEOSFree_r(context, text);
return wkt;
}
 
// Takes ownership of a result and writes it normalised.
std::string result(GEOSGeometry* geometry) {
const Geometry owned = own(geometry);
if (GEOSNormalize_r(context, owned.get()) != 0) fail();
return write(owned.get());
}
 
[[noreturn]] void fail() const { throw std::runtime_error(error.empty() ? "GEOS operation failed" : error); }
 
static void remember(const char* message, void* error) { *static_cast<std::string*>(error) = message; }
};
};
main.js
import { initNative, BufferOp } from './native/buffer.h';
import { Measure } from './native/measure.h';
 
await initNative();
const circle = await BufferOp.around('POINT (0 0)', 10, 8); // 8 segments per quarter circle
console.log((await Measure.area(circle)).toFixed(4));
const road = 'LINESTRING (0 0, 100 0)';
const flat = await BufferOp.withStyle(road, 5, 'flat', 'round');
const round = await BufferOp.withStyle(road, 5, 'round', 'round');
console.log(await Measure.area(flat), (await Measure.area(round)).toFixed(4));
const setback = await BufferOp.around('POLYGON ((0 0, 20 0, 20 20, 0 20, 0 0))', -2, 8);
console.log(setback, await Measure.area(setback));
console.log(await BufferOp.offsetCurve(road, 5));
PRINTSfirst run downloads 1.7 MB
312.1445
1000 1078.0361
POLYGON ((2 2, 2 18, 18 18, 18 2, 2 2)) 256
LINESTRING (0 5, 100 5)

Find why a polygon is invalid and repair it

GEOSisValid_r and GEOSisValidReason_r say whether a shape breaks the OGC rules, what is wrong and where. GEOSMakeValidWithParams_r repairs it two ways, which agree on a bowtie and differ on a hole that leaks out of its shell.

src/native/validity.h
#pragma once
 
#include <geos_c.h>
 
#include <memory>
#include <stdexcept>
#include <string>
 
// OGC validity: whether a shape is valid, what is wrong and where if it is not, and a repair.
class Validity {
public:
static bool isValid(const std::string& wkt) {
Session geos;
const char valid = GEOSisValid_r(geos.context, geos.read(wkt).get());
if (valid == 2) geos.fail();
return valid == 1;
}
 
// "Valid Geometry", or the problem and its location, such as "Self-intersection[5 5]".
static std::string reason(const std::string& wkt) {
Session geos;
char* text = GEOSisValidReason_r(geos.context, geos.read(wkt).get());
if (!text) geos.fail();
const std::string reason = text;
GEOSFree_r(geos.context, text);
return reason;
}
 
// "linework" rebuilds the shape from all of its noded edges; "structure" repairs each ring,
// then keeps shells as shells and subtracts the holes from them.
static std::string makeValid(const std::string& wkt, const std::string& method) {
if (method != "linework" && method != "structure") throw std::invalid_argument("method must be linework or structure");
Session geos;
const auto geometry = geos.read(wkt);
GEOSMakeValidParams* params = GEOSMakeValidParams_create_r(geos.context);
GEOSMakeValidParams_setMethod_r(geos.context, params, method == "structure" ? GEOS_MAKE_VALID_STRUCTURE : GEOS_MAKE_VALID_LINEWORK);
GEOSGeometry* repaired = GEOSMakeValidWithParams_r(geos.context, geometry.get(), params);
GEOSMakeValidParams_destroy_r(geos.context, params);
return geos.result(repaired);
}
 
private:
// GEOS's reentrant C API: each call gets its own context, which also holds the last error message.
struct Session {
struct Destroy {
GEOSContextHandle_t context;
void operator()(GEOSGeometry* geometry) const { GEOSGeom_destroy_r(context, geometry); }
};
using Geometry = std::unique_ptr<GEOSGeometry, Destroy>;
 
GEOSContextHandle_t context = GEOS_init_r();
std::string error;
 
Session() { GEOSContext_setErrorMessageHandler_r(context, remember, &error); }
~Session() { GEOS_finish_r(context); }
Session(const Session&) = delete;
Session& operator=(const Session&) = delete;
 
Geometry own(GEOSGeometry* geometry) {
if (!geometry) fail();
return Geometry(geometry, Destroy{context});
}
 
Geometry read(const std::string& wkt) {
GEOSWKTReader* reader = GEOSWKTReader_create_r(context);
GEOSGeometry* geometry = GEOSWKTReader_read_r(context, reader, wkt.c_str());
GEOSWKTReader_destroy_r(context, reader);
return own(geometry);
}
 
std::string write(const GEOSGeometry* geometry) {
GEOSWKTWriter* writer = GEOSWKTWriter_create_r(context);
char* text = GEOSWKTWriter_write_r(context, writer, geometry);
GEOSWKTWriter_destroy_r(context, writer);
if (!text) fail();
const std::string wkt = text;
GEOSFree_r(context, text);
return wkt;
}
 
// Takes ownership of a result and writes it normalised.
std::string result(GEOSGeometry* geometry) {
const Geometry owned = own(geometry);
if (GEOSNormalize_r(context, owned.get()) != 0) fail();
return write(owned.get());
}
 
[[noreturn]] void fail() const { throw std::runtime_error(error.empty() ? "GEOS operation failed" : error); }
 
static void remember(const char* message, void* error) { *static_cast<std::string*>(error) = message; }
};
};
main.js
import { initNative, Validity } from './native/validity.h';
import { Measure } from './native/measure.h';
 
await initNative();
const bowtie = 'POLYGON ((0 0, 10 10, 10 0, 0 10, 0 0))';
console.log(await Validity.isValid(bowtie), await Validity.reason(bowtie));
const repaired = await Validity.makeValid(bowtie, 'linework');
console.log(repaired, await Validity.isValid(repaired), await Measure.area(repaired));
const leaky = 'POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (5 5, 15 5, 15 15, 5 15, 5 5))';
console.log(await Validity.reason(leaky));
for (const method of ['linework', 'structure']) {
console.log(method, await Measure.area(await Validity.makeValid(leaky, method)));
}
PRINTSfirst run downloads 1.7 MB
false Self-intersection[5 5]
MULTIPOLYGON (((5 5, 10 10, 10 0, 5 5)), ((0 0, 0 10, 5 5, 0 0))) true 50
Self-intersection[5 10]
linework 150
structure 75

A GEOS command-line tool

WASI has no JavaScript bindings: the program is a main() that links the library, built into one .wasm and run with wasmtime. This one reads shapes from WKT files and prints what GEOS makes of them.

crossbind.config.js
import geosWasi from '@crossbind/port-geos-wasi/crossbind.config.js';
 
export default {
general: { name: 'geos-tool' },
dependencies: [geosWasi],
paths: { config: import.meta.url },
};
src/native/main.cpp
// A command-line GEOS for WASI: reads shapes from WKT files and prints results as normalised WKT.
// geos-tool intersection|union|difference <a.wkt> <b.wkt>
// geos-tool validate <a.wkt>
#include <geos_c.h>
 
#include <cstdio>
#include <cstring>
#include <fstream>
#include <memory>
#include <sstream>
#include <string>
 
namespace {
 
GEOSContextHandle_t context;
 
struct Destroy {
void operator()(GEOSGeometry* geometry) const { GEOSGeom_destroy_r(context, geometry); }
};
using Geometry = std::unique_ptr<GEOSGeometry, Destroy>;
 
void report(const char* message, void*) { std::fprintf(stderr, "geos: %s\n", message); }
 
Geometry read(const char* path) {
std::ifstream file(path);
if (!file) {
std::fprintf(stderr, "cannot open %s\n", path);
return nullptr;
}
std::stringstream text;
text << file.rdbuf();
GEOSWKTReader* reader = GEOSWKTReader_create_r(context);
Geometry geometry(GEOSWKTReader_read_r(context, reader, text.str().c_str()));
GEOSWKTReader_destroy_r(context, reader);
return geometry;
}
 
// Prints the result normalised, with its area; false when GEOS failed (its message is on stderr).
bool print(const std::string& label, GEOSGeometry* result) {
if (!result) return false;
const Geometry owned(result);
GEOSNormalize_r(context, owned.get());
GEOSWKTWriter* writer = GEOSWKTWriter_create_r(context);
char* wkt = GEOSWKTWriter_write_r(context, writer, owned.get());
GEOSWKTWriter_destroy_r(context, writer);
double area = 0;
GEOSArea_r(context, owned.get(), &area);
std::printf("%s: %s, area %g\n", label.c_str(), wkt, area);
GEOSFree_r(context, wkt);
return true;
}
 
int validate(const char* path) {
const Geometry shape = read(path);
if (!shape) return 1;
char* reason = GEOSisValidReason_r(context, shape.get());
if (!reason) return 1;
std::printf("validate %s: %s\n", path, reason);
const bool valid = std::strcmp(reason, "Valid Geometry") == 0;
GEOSFree_r(context, reason);
return valid || print("repaired", GEOSMakeValid_r(context, shape.get())) ? 0 : 1;
}
 
int overlay(const std::string& operation, const char* first, const char* second) {
const Geometry a = read(first);
const Geometry b = read(second);
if (!a || !b) return 1;
GEOSGeometry* result = operation == "intersection" ? GEOSIntersection_r(context, a.get(), b.get())
: operation == "union" ? GEOSUnion_r(context, a.get(), b.get())
: GEOSDifference_r(context, a.get(), b.get());
return print(operation + " " + first + " " + second, result) ? 0 : 1;
}
 
int run(int argc, char** argv) {
const std::string operation = argc >= 2 ? argv[1] : "";
if ((operation == "intersection" || operation == "union" || operation == "difference") && argc >= 4) return overlay(operation, argv[2], argv[3]);
if (operation == "validate" && argc >= 3) return validate(argv[2]);
std::fprintf(stderr, "usage: geos-tool intersection|union|difference <a.wkt> <b.wkt>\n geos-tool validate <a.wkt>\n");
return 2;
}
 
} // namespace
 
int main(int argc, char** argv) {
context = GEOS_init_r();
GEOSContext_setErrorMessageHandler_r(context, report, nullptr);
const int code = run(argc, argv);
GEOS_finish_r(context);
return code;
}
shell
npx crossbind build -p wasi -e wasi -b release
wasmtime run --dir=. .crossbind/build/geos-tool-wasi-wasm32-st-release.wasm intersection parcel.wkt flood.wkt
wasmtime run --dir=. .crossbind/build/geos-tool-wasi-wasm32-st-release.wasm difference parcel.wkt flood.wkt
wasmtime run --dir=. .crossbind/build/geos-tool-wasi-wasm32-st-release.wasm validate bowtie.wkt
output
intersection parcel.wkt flood.wkt: POLYGON ((5 5, 5 10, 10 10, 10 5, 5 5)), area 25
difference parcel.wkt flood.wkt: POLYGON ((0 0, 0 10, 5 10, 5 5, 10 5, 10 0, 0 0)), area 75
validate bowtie.wkt: Self-intersection[5 5]
repaired: MULTIPOLYGON (((5 5, 10 10, 10 0, 5 5)), ((0 0, 0 10, 5 5, 0 0))), area 50
Type to search every guide page and section.
↑↓ navigate↵ openesc close