Add Buffer operation

CHANGELOG.md

@@ -2,6 +2,13 @@
 
 ## Unreleased
 
+- Add `Buffer` function that computes the buffer of a geometry at a given
+  distance. Options are available for controlling quad segments, end cap style
+  (round, flat, square), join style (round, mitre, bevel), single-sided mode,
+  and simplify factor. The implementation is based on a port of JTS. This means
+  that `Buffer` is now available natively in Go without need for the GEOS
+  dependency.
+
 - Change `GeometryCollection.Dimension()` to return -1 for empty geometry
   collections (previously returned 0). This change is consistent with GEOS.
 

geom/alg_buffer.go

@@ -0,0 +1,126 @@
+package geom
+
+import (
+	"github.com/peterstace/simplefeatures/internal/jtsport/jts"
+)
+
+// Buffer returns a geometry that contains all points within the given distance
+// of the input geometry.
+//
+// In GIS, the positive (or negative) buffer of a geometry is defined as
+// the Minkowski sum (or difference) of the geometry with a circle of radius
+// equal to the absolute value of the buffer distance.
+//
+// The buffer operation always returns a polygonal result.
+// The negative or zero-distance buffer of lines and points is always an empty
+// Polygon.
+//
+// Since true buffer curves may contain circular arcs, computed buffer polygons
+// are only approximations to the true geometry. The user can control the
+// accuracy of the approximation by specifying the number of segments used to
+// approximate quarter circles (via [BufferQuadSegments]).
+//
+// An error may be returned in pathological cases of numerical degeneracy.
+func Buffer(g Geometry, distance float64, opts ...BufferOption) (Geometry, error) {
+	params := jts.OperationBuffer_NewBufferParameters()
+	for _, opt := range opts {
+		opt(params)
+	}
+
+	var result Geometry
+	err := catch(func() error {
+		wkbReader := jts.Io_NewWKBReader()
+		jtsG, err := wkbReader.ReadBytes(g.AsBinary())
+		if err != nil {
+			return wrap(err, "converting geometry to JTS")
+		}
+		jtsResult := jts.OperationBuffer_BufferOp_BufferOpWithParams(jtsG, distance, params)
+		wkbWriter := jts.Io_NewWKBWriter()
+		result, err = UnmarshalWKB(wkbWriter.Write(jtsResult), NoValidate{})
+		return wrap(err, "converting JTS buffer result to simplefeatures")
+	})
+	return result, err
+}
+
+// BufferOption allows the behaviour of the [Buffer] operation to be modified.
+type BufferOption func(*jts.OperationBuffer_BufferParameters)
+
+// BufferQuadSegments sets the number of segments used to approximate a quarter
+// circle. It defaults to 8.
+func BufferQuadSegments(quadSegments int) BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetQuadrantSegments(quadSegments)
+	}
+}
+
+// BufferEndCapRound sets the end cap style to 'round'. It is 'round' by
+// default.
+func BufferEndCapRound() BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetEndCapStyle(jts.OperationBuffer_BufferParameters_CAP_ROUND)
+	}
+}
+
+// BufferEndCapFlat sets the end cap style to 'flat'. It is 'round' by default.
+func BufferEndCapFlat() BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetEndCapStyle(jts.OperationBuffer_BufferParameters_CAP_FLAT)
+	}
+}
+
+// BufferEndCapSquare sets the end cap style to 'square'. It is 'round' by
+// default.
+func BufferEndCapSquare() BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetEndCapStyle(jts.OperationBuffer_BufferParameters_CAP_SQUARE)
+	}
+}
+
+// BufferJoinStyleRound sets the join style to 'round'. It is 'round' by
+// default.
+func BufferJoinStyleRound() BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetJoinStyle(jts.OperationBuffer_BufferParameters_JOIN_ROUND)
+	}
+}
+
+// BufferJoinStyleMitre sets the join style to 'mitre'. It is 'round' by
+// default. The mitreLimit controls how far a mitre join can extend from the
+// join point. Corners with a ratio which exceed the limit will be beveled.
+func BufferJoinStyleMitre(mitreLimit float64) BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetJoinStyle(jts.OperationBuffer_BufferParameters_JOIN_MITRE)
+		p.SetMitreLimit(mitreLimit)
+	}
+}
+
+// BufferJoinStyleBevel sets the join style to 'bevel'. It is 'round' by
+// default.
+func BufferJoinStyleBevel() BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetJoinStyle(jts.OperationBuffer_BufferParameters_JOIN_BEVEL)
+	}
+}
+
+// BufferSingleSided sets whether the computed buffer should be single-sided.
+// A single-sided buffer is constructed on only one side of each input line.
+// The side is determined by the sign of the buffer distance: positive
+// indicates the left-hand side, negative indicates the right-hand side.
+// The end cap style is ignored for single-sided buffers and forced to flat.
+func BufferSingleSided() BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetSingleSided(true)
+	}
+}
+
+// BufferSimplifyFactor sets the factor used to determine the simplify distance
+// tolerance for input simplification. The factor is multiplied by the buffer
+// distance to get the simplification tolerance. Simplifying can increase the
+// performance of computing buffers. Values between 0.01 and 0.1 produce
+// relatively good accuracy for the generated buffer. Larger values sacrifice
+// accuracy in return for performance. The default is 0.01.
+func BufferSimplifyFactor(factor float64) BufferOption {
+	return func(p *jts.OperationBuffer_BufferParameters) {
+		p.SetSimplifyFactor(factor)
+	}
+}

geom/alg_buffer_test.go

@@ -0,0 +1,167 @@
+package geom_test
+
+import (
+	"strconv"
+	"testing"
+
+	"github.com/peterstace/simplefeatures/geom"
+	"github.com/peterstace/simplefeatures/internal/test"
+)
+
+func TestBuffer(t *testing.T) {
+	for i, tc := range []struct {
+		name  string
+		input string
+		dist  float64
+		opts  []geom.BufferOption
+		want  string
+	}{
+		{
+			name:  "point with defaults",
+			input: "POINT(0 0)",
+			dist:  1,
+			opts:  nil,
+			want:  "POLYGON((1 0,0.9807852804032304 -0.19509032201612825,0.9238795325112867 -0.3826834323650898,0.8314696123025453 -0.5555702330196022,0.7071067811865476 -0.7071067811865475,0.5555702330196023 -0.8314696123025452,0.38268343236508984 -0.9238795325112867,0.1950903220161283 -0.9807852804032304,0 -1,-0.1950903220161282 -0.9807852804032304,-0.3826834323650897 -0.9238795325112867,-0.555570233019602 -0.8314696123025453,-0.7071067811865475 -0.7071067811865476,-0.8314696123025453 -0.5555702330196021,-0.9238795325112867 -0.38268343236508984,-0.9807852804032303 -0.19509032201612858,-1 0,-0.9807852804032304 0.19509032201612836,-0.9238795325112868 0.38268343236508967,-0.8314696123025453 0.555570233019602,-0.7071067811865477 0.7071067811865475,-0.5555702330196022 0.8314696123025453,-0.38268343236509034 0.9238795325112865,-0.19509032201612866 0.9807852804032303,0 1,0.19509032201612828 0.9807852804032304,0.38268343236509 0.9238795325112866,0.5555702330196018 0.8314696123025455,0.7071067811865475 0.7071067811865477,0.8314696123025453 0.5555702330196022,0.9238795325112865 0.3826834323650904,0.9807852804032303 0.19509032201612872,1 0))",
+		},
+		{
+			name:  "distance of 2",
+			input: "POINT(0 0)",
+			dist:  2,
+			opts:  nil,
+			want:  "POLYGON((2 0,1.9615705608064609 -0.3901806440322565,1.8477590650225735 -0.7653668647301796,1.6629392246050907 -1.1111404660392044,1.4142135623730951 -1.414213562373095,1.1111404660392046 -1.6629392246050905,0.7653668647301797 -1.8477590650225735,0.3901806440322566 -1.9615705608064609,0 -2,-0.3901806440322564 -1.9615705608064609,-0.7653668647301795 -1.8477590650225735,-1.111140466039204 -1.6629392246050907,-1.414213562373095 -1.4142135623730951,-1.6629392246050907 -1.1111404660392041,-1.8477590650225735 -0.7653668647301797,-1.9615705608064606 -0.39018064403225716,-2 0,-1.9615705608064609 0.3901806440322567,-1.8477590650225737 0.7653668647301793,-1.6629392246050907 1.111140466039204,-1.4142135623730954 1.414213562373095,-1.1111404660392044 1.6629392246050907,-0.7653668647301807 1.847759065022573,-0.39018064403225733 1.9615705608064606,0 2,0.39018064403225655 1.9615705608064609,0.76536686473018 1.8477590650225733,1.1111404660392037 1.662939224605091,1.414213562373095 1.4142135623730954,1.6629392246050907 1.1111404660392044,1.847759065022573 0.7653668647301808,1.9615705608064606 0.39018064403225744,2 0))",
+		},
+		{
+			name:  "reduced quads",
+			input: "POINT(0 0)",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferQuadSegments(2)},
+			want:  "POLYGON((1 0,0.7071067811865476 -0.7071067811865475,0 -1,-0.7071067811865475 -0.7071067811865476,-1 0,-0.7071067811865477 0.7071067811865475,0 1,0.7071067811865475 0.7071067811865477,1 0))",
+		},
+		{
+			name:  "multipoint with distance apart greater than distance",
+			input: "MULTIPOINT(0 0,3 0)",
+			dist:  1,
+			opts:  nil,
+			want:  "MULTIPOLYGON(((4 0,3.9807852804032304 -0.19509032201612825,3.923879532511287 -0.3826834323650898,3.8314696123025453 -0.5555702330196022,3.7071067811865475 -0.7071067811865475,3.555570233019602 -0.8314696123025452,3.3826834323650896 -0.9238795325112867,3.1950903220161284 -0.9807852804032304,3 -1,2.804909677983872 -0.9807852804032304,2.6173165676349104 -0.9238795325112867,2.444429766980398 -0.8314696123025453,2.2928932188134525 -0.7071067811865476,2.1685303876974547 -0.5555702330196021,2.076120467488713 -0.38268343236508984,2.0192147195967696 -0.19509032201612858,2 0,2.0192147195967696 0.19509032201612836,2.076120467488713 0.38268343236508967,2.1685303876974547 0.555570233019602,2.292893218813452 0.7071067811865475,2.444429766980398 0.8314696123025453,2.6173165676349095 0.9238795325112865,2.804909677983871 0.9807852804032303,3 1,3.1950903220161284 0.9807852804032304,3.38268343236509 0.9238795325112866,3.5555702330196017 0.8314696123025455,3.7071067811865475 0.7071067811865477,3.8314696123025453 0.5555702330196022,3.9238795325112865 0.3826834323650904,3.9807852804032304 0.19509032201612872,4 0)),((1 0,0.9807852804032304 -0.19509032201612825,0.9238795325112867 -0.3826834323650898,0.8314696123025453 -0.5555702330196022,0.7071067811865476 -0.7071067811865475,0.5555702330196023 -0.8314696123025452,0.38268343236508984 -0.9238795325112867,0.1950903220161283 -0.9807852804032304,0 -1,-0.1950903220161282 -0.9807852804032304,-0.3826834323650897 -0.9238795325112867,-0.555570233019602 -0.8314696123025453,-0.7071067811865475 -0.7071067811865476,-0.8314696123025453 -0.5555702330196021,-0.9238795325112867 -0.38268343236508984,-0.9807852804032303 -0.19509032201612858,-1 0,-0.9807852804032304 0.19509032201612836,-0.9238795325112868 0.38268343236508967,-0.8314696123025453 0.555570233019602,-0.7071067811865477 0.7071067811865475,-0.5555702330196022 0.8314696123025453,-0.38268343236509034 0.9238795325112865,-0.19509032201612866 0.9807852804032303,0 1,0.19509032201612828 0.9807852804032304,0.38268343236509 0.9238795325112866,0.5555702330196018 0.8314696123025455,0.7071067811865475 0.7071067811865477,0.8314696123025453 0.5555702330196022,0.9238795325112865 0.3826834323650904,0.9807852804032303 0.19509032201612872,1 0)))",
+		},
+		{
+			name:  "multipoint with distance apart less than distance",
+			input: "MULTIPOINT(0 0,1 0)",
+			dist:  1,
+			opts:  nil,
+			want:  "POLYGON((0.49999999999999994 0.861172520678903,0.6173165676349097 0.9238795325112865,0.8049096779838714 0.9807852804032303,1 1,1.1950903220161282 0.9807852804032304,1.38268343236509 0.9238795325112866,1.5555702330196017 0.8314696123025455,1.7071067811865475 0.7071067811865477,1.8314696123025453 0.5555702330196022,1.9238795325112865 0.3826834323650904,1.9807852804032304 0.19509032201612872,2 0,1.9807852804032304 -0.19509032201612825,1.9238795325112867 -0.3826834323650898,1.8314696123025453 -0.5555702330196022,1.7071067811865475 -0.7071067811865475,1.5555702330196022 -0.8314696123025452,1.3826834323650898 -0.9238795325112867,1.1950903220161284 -0.9807852804032304,1 -1,0.8049096779838718 -0.9807852804032304,0.6173165676349103 -0.9238795325112867,0.5000000000000001 -0.861172520678903,0.38268343236508984 -0.9238795325112867,0.1950903220161283 -0.9807852804032304,0 -1,-0.1950903220161282 -0.9807852804032304,-0.3826834323650897 -0.9238795325112867,-0.555570233019602 -0.8314696123025453,-0.7071067811865475 -0.7071067811865476,-0.8314696123025453 -0.5555702330196021,-0.9238795325112867 -0.38268343236508984,-0.9807852804032303 -0.19509032201612858,-1 0,-0.9807852804032304 0.19509032201612836,-0.9238795325112868 0.38268343236508967,-0.8314696123025453 0.555570233019602,-0.7071067811865477 0.7071067811865475,-0.5555702330196022 0.8314696123025453,-0.38268343236509034 0.9238795325112865,-0.19509032201612866 0.9807852804032303,0 1,0.19509032201612828 0.9807852804032304,0.38268343236509 0.9238795325112866,0.49999999999999994 0.861172520678903))",
+		},
+		{
+			name:  "linestring with defaults",
+			input: "LINESTRING(0 0,1 0)",
+			dist:  1,
+			opts:  nil,
+			want:  "POLYGON((1 1,1.1950903220161284 0.9807852804032304,1.3826834323650898 0.9238795325112867,1.5555702330196022 0.8314696123025452,1.7071067811865475 0.7071067811865475,1.8314696123025453 0.5555702330196022,1.9238795325112867 0.3826834323650898,1.9807852804032304 0.19509032201612825,2 0,1.9807852804032304 -0.19509032201612825,1.9238795325112867 -0.3826834323650898,1.8314696123025453 -0.5555702330196022,1.7071067811865475 -0.7071067811865475,1.5555702330196022 -0.8314696123025452,1.3826834323650898 -0.9238795325112867,1.1950903220161284 -0.9807852804032304,1 -1,0 -1,-0.19509032201612866 -0.9807852804032303,-0.38268343236509034 -0.9238795325112865,-0.5555702330196022 -0.8314696123025453,-0.7071067811865477 -0.7071067811865475,-0.8314696123025453 -0.555570233019602,-0.9238795325112868 -0.38268343236508967,-0.9807852804032304 -0.19509032201612836,-1 0,-0.9807852804032303 0.19509032201612858,-0.9238795325112867 0.38268343236508984,-0.8314696123025453 0.5555702330196021,-0.7071067811865475 0.7071067811865476,-0.555570233019602 0.8314696123025453,-0.3826834323650897 0.9238795325112867,-0.1950903220161282 0.9807852804032304,0 1,1 1))",
+		},
+		{
+			name:  "endcap flat",
+			input: "LINESTRING(0 0,1 0)",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferEndCapFlat()},
+			want:  "POLYGON((1 1,1 -1,0 -1,0 1,1 1))",
+		},
+		{
+			name:  "endcap square",
+			input: "LINESTRING(0 0,1 0)",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferEndCapSquare()},
+			want:  "POLYGON((1 1,2 1,2 -1,0 -1,-1 -1,-1 1,1 1))",
+		},
+		{
+			name:  "single sided positive",
+			input: "LINESTRING(0 0,1 0)",
+			dist:  +1,
+			opts:  []geom.BufferOption{geom.BufferSingleSided()},
+			want:  "POLYGON((1 0,0 0,0 1,1 1,1 0))",
+		},
+		{
+			name:  "single sided negative",
+			input: "LINESTRING(0 0,1 0)",
+			dist:  -1,
+			opts:  []geom.BufferOption{geom.BufferSingleSided()},
+			want:  "POLYGON((0 0,1 0,1 -1,0 -1,0 0))",
+		},
+		{
+			name:  "join style miter below limit",
+			input: "LINESTRING(0 0,3 0,3 3)",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferEndCapFlat(), geom.BufferJoinStyleMitre(2)},
+			want:  "POLYGON((2 1,2 3,4 3,4 -1,0 -1,0 1,2 1))",
+		},
+		{
+			name:  "join style miter above limit",
+			input: "LINESTRING(0 0,3 0,3 3)",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferEndCapFlat(), geom.BufferJoinStyleMitre(1)},
+			want:  "POLYGON((2 1,2 3,4 3,4 -0.4142135623730951,3.414213562373095 -1,0 -1,0 1,2 1))",
+		},
+		{
+			name:  "join style bevel",
+			input: "LINESTRING(0 0,3 0,3 3)",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferEndCapFlat(), geom.BufferJoinStyleBevel()},
+			want:  "POLYGON((2 1,2 3,4 3,4 0,3 -1,0 -1,0 1,2 1))",
+		},
+		{
+			name:  "simplify factor high",
+			input: "LINESTRING(0 0,0.25 0.01,0.5 0,0.75 0.01,1 0)",
+			dist:  0.1,
+			opts:  []geom.BufferOption{geom.BufferSimplifyFactor(1.0), geom.BufferEndCapFlat()},
+			want:  "POLYGON((0.25 0.11,0.75 0.11,1.0039968038348872 0.09992009587217894,0.9960031961651128 -0.09992009587217894,0.5 -0.1,0.25 -0.09007996802557443,0.003996803834887157 -0.09992009587217894,-0.003996803834887157 0.09992009587217894,0.25 0.11))",
+		},
+		{
+			name:  "simplify factor low",
+			input: "LINESTRING(0 0,0.25 0.01,0.5 0,0.75 0.01,1 0)",
+			dist:  0.1,
+			opts:  []geom.BufferOption{geom.BufferSimplifyFactor(0.01), geom.BufferEndCapFlat()},
+			want:  "POLYGON((0.24600319616511285 0.10992009587217894,0.25399680383488715 0.10992009587217894,0.5 0.10007996802557442,0.7460031961651128 0.10992009587217894,0.7539968038348872 0.10992009587217894,1.0039968038348872 0.09992009587217894,0.9960031961651128 -0.09992009587217894,0.75 -0.09007996802557443,0.5039968038348872 -0.09992009587217894,0.49600319616511285 -0.09992009587217894,0.25 -0.09007996802557443,0.003996803834887157 -0.09992009587217894,-0.003996803834887157 0.09992009587217894,0.24600319616511285 0.10992009587217894))",
+		},
+		{
+			name:  "polygon with hole closed by buffer",
+			input: "POLYGON((0 0,10 0,10 10,0 10,0 0),(4 4,6 4,6 6,4 6,4 4))",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferQuadSegments(2)},
+			want:  "POLYGON((-1 0,-1 10,-0.7071067811865475 10.707106781186548,0 11,10 11,10.707106781186548 10.707106781186548,11 10,11 0,10.707106781186548 -0.7071067811865475,10 -1,0 -1,-0.7071067811865475 -0.7071067811865476,-1 0))",
+		},
+		{
+			name:  "polygon with hole not closed by buffer",
+			input: "POLYGON((0 0,10 0,10 10,0 10,0 0),(3 3,7 3,7 7,3 7,3 3))",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferQuadSegments(2)},
+			want:  "POLYGON((-1 0,-1 10,-0.7071067811865475 10.707106781186548,0 11,10 11,10.707106781186548 10.707106781186548,11 10,11 0,10.707106781186548 -0.7071067811865475,10 -1,0 -1,-0.7071067811865475 -0.7071067811865476,-1 0),(4 4,6 4,6 6,4 6,4 4))",
+		},
+		{
+			name:  "negative buffer",
+			input: "POLYGON((0 0,10 0,10 10,0 10,0 0))",
+			dist:  -1,
+			opts:  []geom.BufferOption{geom.BufferQuadSegments(2)},
+			want:  "POLYGON((1 1,1 9,9 9,9 1,1 1))",
+		},
+		{
+			name:  "negative buffer causing disappearance",
+			input: "POLYGON((0 0,2 0,2 2,0 2,0 0))",
+			dist:  -2,
+			opts:  []geom.BufferOption{geom.BufferQuadSegments(2)},
+			want:  "POLYGON EMPTY",
+		},
+		{
+			name:  "geometry collection",
+			input: "GEOMETRYCOLLECTION(POLYGON((0 0,4 0,4 4,0 4,0 0)),LINESTRING(6 2,10 2),POINT(12 2))",
+			dist:  1,
+			opts:  []geom.BufferOption{geom.BufferQuadSegments(2)},
+			want:  "MULTIPOLYGON(((-1 0,-1 4,-0.7071067811865475 4.707106781186548,0 5,4 5,4.707106781186548 4.707106781186548,5 4,5 2,5 0,4.707106781186548 -0.7071067811865475,4 -1,0 -1,-0.7071067811865475 -0.7071067811865476,-1 0)),((5 2,5.292893218813452 2.7071067811865475,6 3,10 3,10.707106781186548 2.7071067811865475,11 2,10.707106781186548 1.2928932188134525,10 1,6 1,5.292893218813452 1.2928932188134525,5 2)),((11 2,11.292893218813452 2.7071067811865475,12 3,12.707106781186548 2.707106781186548,13 2,12.707106781186548 1.2928932188134525,12 1,11.292893218813452 1.2928932188134525,11 2)))",
+		},
+	} {
+		t.Run(strconv.Itoa(i), func(t *testing.T) {
+			in := test.FromWKT(t, tc.input)
+			got, err := geom.Buffer(in, tc.dist, tc.opts...)
+			test.NoErr(t, err)
+			test.ExactEqualsWKT(t, got, tc.want, geom.ToleranceXY(1e-10))
+		})
+	}
+}