Commit 325b6154 authored by abhishek20123g's avatar abhishek20123g

geo/geomfn: Implements ST_Segmentize for geometry

Fixes https://github.com/cockroachdb/cockroach/issues/49029

This PR implements ST_Segmentize({geometry, float8}) builtin
function, which allows modify given geometry such that no
segment longer than the given max_segment_length.

Also this PR refactors and add some extra test cases for
ST_Segmentize for geography.

Release note (sql change): This PR implements ST_Segmentize({geometry,
float8}) builtin function.
parent 1be1ec06
......@@ -1114,6 +1114,8 @@ given Geometry.</p>
<p>The calculations are done on a sphere.</p>
<p>This function utilizes the S2 library for spherical calculations.</p>
</span></td></tr>
<tr><td><a name="st_segmentize"></a><code>st_segmentize(geometry: geometry, max_segment_length: <a href="float.html">float</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns a modified Geometry having no segment longer than the given max_segment_length. Length units are in units of spatial reference.</p>
</span></td></tr>
<tr><td><a name="st_setsrid"></a><code>st_setsrid(geography: geography, srid: <a href="int.html">int</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>Sets a Geography to a new SRID without transforming the coordinates.</p>
</span></td></tr>
<tr><td><a name="st_setsrid"></a><code>st_setsrid(geometry: geometry, srid: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Sets a Geometry to a new SRID without transforming the coordinates.</p>
......
......@@ -15,6 +15,7 @@ import (
"github.com/cockroachdb/cockroach/pkg/geo"
"github.com/cockroachdb/cockroach/pkg/geo/geographiclib"
"github.com/cockroachdb/cockroach/pkg/geo/geosegmentize"
"github.com/cockroachdb/errors"
"github.com/golang/geo/s2"
"github.com/twpayne/go-geom"
......@@ -40,7 +41,7 @@ func Segmentize(geography *geo.Geography, segmentMaxLength float64) (*geo.Geogra
// Convert segmentMaxLength to Angle with respect to earth sphere as
// further calculation is done considering segmentMaxLength as Angle.
segmentMaxAngle := segmentMaxLength / spheroid.SphereRadius
segGeometry, err := segmentizeGeom(geometry, segmentMaxAngle)
segGeometry, err := geosegmentize.SegmentizeGeom(geometry, segmentMaxAngle, segmentizeCoords)
if err != nil {
return nil, err
}
......@@ -48,123 +49,35 @@ func Segmentize(geography *geo.Geography, segmentMaxLength float64) (*geo.Geogra
}
}
// segmentizeGeom returns a modified geom.T having no segment longer than
// the given maximum segment length.
func segmentizeGeom(geometry geom.T, segmentMaxAngle float64) (geom.T, error) {
if geometry.Empty() {
return geometry, nil
}
switch geometry := geometry.(type) {
case *geom.Point, *geom.MultiPoint:
return geometry, nil
case *geom.LineString:
var allFlatCoordinates []float64
for pointIdx := 1; pointIdx < geometry.NumCoords(); pointIdx++ {
allFlatCoordinates = append(
allFlatCoordinates,
segmentizeCoords(geometry.Coord(pointIdx-1), geometry.Coord(pointIdx), segmentMaxAngle)...,
)
}
// Appending end point as it wasn't included in the iteration of coordinates.
allFlatCoordinates = append(allFlatCoordinates, geometry.Coord(geometry.NumCoords()-1)...)
return geom.NewLineStringFlat(geom.XY, allFlatCoordinates).SetSRID(geometry.SRID()), nil
case *geom.MultiLineString:
segMultiLine := geom.NewMultiLineString(geom.XY).SetSRID(geometry.SRID())
for lineIdx := 0; lineIdx < geometry.NumLineStrings(); lineIdx++ {
l, err := segmentizeGeom(geometry.LineString(lineIdx), segmentMaxAngle)
if err != nil {
return nil, err
}
err = segMultiLine.Push(l.(*geom.LineString))
if err != nil {
return nil, err
}
}
return segMultiLine, nil
case *geom.LinearRing:
var allFlatCoordinates []float64
for pointIdx := 1; pointIdx < geometry.NumCoords(); pointIdx++ {
allFlatCoordinates = append(
allFlatCoordinates,
segmentizeCoords(geometry.Coord(pointIdx-1), geometry.Coord(pointIdx), segmentMaxAngle)...,
)
}
// Appending end point as it wasn't included in the iteration of coordinates.
allFlatCoordinates = append(allFlatCoordinates, geometry.Coord(geometry.NumCoords()-1)...)
return geom.NewLinearRingFlat(geom.XY, allFlatCoordinates).SetSRID(geometry.SRID()), nil
case *geom.Polygon:
segPolygon := geom.NewPolygon(geom.XY).SetSRID(geometry.SRID())
for loopIdx := 0; loopIdx < geometry.NumLinearRings(); loopIdx++ {
l, err := segmentizeGeom(geometry.LinearRing(loopIdx), segmentMaxAngle)
if err != nil {
return nil, err
}
err = segPolygon.Push(l.(*geom.LinearRing))
if err != nil {
return nil, err
}
}
return segPolygon, nil
case *geom.MultiPolygon:
segMultiPolygon := geom.NewMultiPolygon(geom.XY).SetSRID(geometry.SRID())
for polygonIdx := 0; polygonIdx < geometry.NumPolygons(); polygonIdx++ {
p, err := segmentizeGeom(geometry.Polygon(polygonIdx), segmentMaxAngle)
if err != nil {
return nil, err
}
err = segMultiPolygon.Push(p.(*geom.Polygon))
if err != nil {
return nil, err
}
}
return segMultiPolygon, nil
case *geom.GeometryCollection:
segGeomCollection := geom.NewGeometryCollection().SetSRID(geometry.SRID())
for geoIdx := 0; geoIdx < geometry.NumGeoms(); geoIdx++ {
g, err := segmentizeGeom(geometry.Geom(geoIdx), segmentMaxAngle)
if err != nil {
return nil, err
}
err = segGeomCollection.Push(g)
if err != nil {
return nil, err
}
}
return segGeomCollection, nil
}
return nil, errors.Newf("unknown type: %T", geometry)
}
// segmentizeCoords inserts multiple points between given two-coordinate and
// return resultant points as flat []float64. Such that distance between any two
// segmentizeCoords inserts multiple points between given two coordinates and
// return resultant point as flat []float64. Such that distance between any two
// points is less than given maximum segment's length, the total number of
// segments is the power of 2, and all the segments are of the same length.
// NOTE: List of points does not consist of end point.
// Note: List of points does not consist of end point.
func segmentizeCoords(a geom.Coord, b geom.Coord, segmentMaxAngle float64) []float64 {
// Converted geom.Coord into s2.Point so we can segmentize the coordinates.
pointA := s2.PointFromLatLng(s2.LatLngFromDegrees(a.Y(), a.X()))
pointB := s2.PointFromLatLng(s2.LatLngFromDegrees(b.Y(), b.X()))
allSegmentizedCoordinates := a.Clone()
chordAngleBetweenPoints := s2.ChordAngleBetweenPoints(pointA, pointB).Angle().Radians()
if segmentMaxAngle <= chordAngleBetweenPoints {
// This calculation is to determine the total number of segment between given
// 2 coordinates, ensuring that the segments are divided into parts divisible by
// a power of 2.
//
// For that fraction by segment must be less than or equal to
// the fraction of max segment length to distance between point, since the
// total number of segment must be power of 2 therefore we can write as
// 1 / (2^n)[numberOfSegmentToCreate] <= segmentMaxLength / distanceBetweenPoints < 1 / (2^(n-1))
// (2^n)[numberOfSegmentToCreate] >= distanceBetweenPoints / segmentMaxLength > 2^(n-1)
// therefore n = ceil(log2(segmentMaxLength/distanceBetweenPoints)). Hence
// numberOfSegmentToCreate = 2^(ceil(log2(segmentMaxLength/distanceBetweenPoints))).
numberOfSegmentToCreate := int(math.Pow(2, math.Ceil(math.Log2(chordAngleBetweenPoints/segmentMaxAngle))))
for pointInserted := 1; pointInserted < numberOfSegmentToCreate; pointInserted++ {
newPoint := s2.Interpolate(float64(pointInserted)/float64(numberOfSegmentToCreate), pointA, pointB)
latLng := s2.LatLngFromPoint(newPoint)
allSegmentizedCoordinates = append(allSegmentizedCoordinates, latLng.Lng.Degrees(), latLng.Lat.Degrees())
}
// This calculation is to determine the total number of segment between given
// 2 coordinates, ensuring that the segments are divided into parts divisible by
// a power of 2.
//
// For that fraction by segment must be less than or equal to
// the fraction of max segment length to distance between point, since the
// total number of segment must be power of 2 therefore we can write as
// 1 / (2^n)[numberOfSegmentToCreate] <= segmentMaxLength / distanceBetweenPoints < 1 / (2^(n-1))
// (2^n)[numberOfSegmentToCreate] >= distanceBetweenPoints / segmentMaxLength > 2^(n-1)
// therefore n = ceil(log2(segmentMaxLength/distanceBetweenPoints)). Hence
// numberOfSegmentToCreate = 2^(ceil(log2(segmentMaxLength/distanceBetweenPoints))).
numberOfSegmentToCreate := int(math.Pow(2, math.Ceil(math.Log2(chordAngleBetweenPoints/segmentMaxAngle))))
allSegmentizedCoordinates := make([]float64, 0, 2*(1+numberOfSegmentToCreate))
allSegmentizedCoordinates = append(allSegmentizedCoordinates, a.Clone()...)
for pointInserted := 1; pointInserted < numberOfSegmentToCreate; pointInserted++ {
newPoint := s2.Interpolate(float64(pointInserted)/float64(numberOfSegmentToCreate), pointA, pointB)
latLng := s2.LatLngFromPoint(newPoint)
allSegmentizedCoordinates = append(allSegmentizedCoordinates, latLng.Lng.Degrees(), latLng.Lat.Degrees())
}
return allSegmentizedCoordinates
}
......@@ -97,7 +97,7 @@ func TestSegmentize(t *testing.T) {
},
}
for _, test := range segmentizeTestCases {
t.Run(fmt.Sprintf("%s, maximum segment length: %v", test.wkt, test.maxSegmentLength), func(t *testing.T) {
t.Run(fmt.Sprintf("%s, maximum segment length: %f", test.wkt, test.maxSegmentLength), func(t *testing.T) {
geog, err := geo.ParseGeography(test.wkt)
require.NoError(t, err)
modifiedGeog, err := Segmentize(geog, test.maxSegmentLength)
......@@ -107,9 +107,9 @@ func TestSegmentize(t *testing.T) {
require.Equal(t, expectedGeog, modifiedGeog)
})
}
// Test for segment maximum length as negative for geometry collection.
t.Run(fmt.Sprintf("%s, maximum segment length: %v", segmentizeTestCases[9].wkt, 0.0), func(t *testing.T) {
geog, err := geo.ParseGeography(segmentizeTestCases[9].wkt)
// Test for segment maximum length as negative.
t.Run("Error when maximum segment length is less than 0", func(t *testing.T) {
geog, err := geo.ParseGeography("MULTILINESTRING ((0 0, 1 1, 5 5), (5 5, 0 0))")
require.NoError(t, err)
_, err = Segmentize(geog, 0)
require.EqualError(t, err, "maximum segment length must be positive")
......@@ -152,6 +152,20 @@ func TestSegmentizeCoords(t *testing.T) {
segmentMaxLength: -1,
resultantCoordinates: []float64{85, 85},
},
{
desc: `Coordinate(0, 0) to Coordinate(0, 0), 0.29`,
a: geom.Coord{0, 0},
b: geom.Coord{0, 0},
segmentMaxLength: 0.29,
resultantCoordinates: []float64{0, 0},
},
{
desc: `Coordinate(85, 85) to Coordinate(0, 0), 1.563200444168918`,
a: geom.Coord{85, 85},
b: geom.Coord{0, 0},
segmentMaxLength: 1.563200444168918,
resultantCoordinates: []float64{85, 85},
},
}
for _, test := range testCases {
t.Run(test.desc, func(t *testing.T) {
......
// Copyright 2020 The Cockroach Authors.
//
// Use of this software is governed by the Business Source License
// included in the file licenses/BSL.txt.
//
// As of the Change Date specified in that file, in accordance with
// the Business Source License, use of this software will be governed
// by the Apache License, Version 2.0, included in the file
// licenses/APL.txt.
package geomfn
import (
"math"
"github.com/cockroachdb/cockroach/pkg/geo"
"github.com/cockroachdb/cockroach/pkg/geo/geosegmentize"
"github.com/cockroachdb/errors"
"github.com/twpayne/go-geom"
)
// Segmentize return modified Geometry having no segment longer
// that given maximum segment length.
// This works by inserting the extra points in such a manner that
// minimum number of new segments with equal length is created,
// between given two-points such that each segment has length less
// than or equal to given maximum segment length.
func Segmentize(g *geo.Geometry, segmentMaxLength float64) (*geo.Geometry, error) {
geometry, err := g.AsGeomT()
if err != nil {
return nil, err
}
switch geometry := geometry.(type) {
case *geom.Point, *geom.MultiPoint:
return g, nil
default:
if segmentMaxLength <= 0 {
return nil, errors.Newf("maximum segment length must be positive")
}
segGeometry, err := geosegmentize.SegmentizeGeom(geometry, segmentMaxLength, segmentizeCoords)
if err != nil {
return nil, err
}
return geo.NewGeometryFromGeom(segGeometry)
}
}
// segmentizeCoords inserts multiple points between given two coordinates and
// return resultant point as flat []float64. Points are inserted in such a
// way that they create minimum number segments of equal length such that each
// segment has a length less than or equal to given maximum segment length.
// Note: List of points does not consist of end point.
func segmentizeCoords(a geom.Coord, b geom.Coord, maxSegmentLength float64) []float64 {
distanceBetweenPoints := math.Sqrt(math.Pow(a.X()-b.X(), 2) + math.Pow(b.Y()-a.Y(), 2))
// numberOfSegmentToCreate represent the total number of segments
// in which given two coordinates will be divided.
numberOfSegmentToCreate := int(math.Ceil(distanceBetweenPoints / maxSegmentLength))
// segmentFraction represent the fraction of length each segment
// has with respect to total length between two coordinates.
allSegmentizedCoordinates := make([]float64, 0, 2*(1+numberOfSegmentToCreate))
allSegmentizedCoordinates = append(allSegmentizedCoordinates, a.Clone()...)
segmentFraction := 1.0 / float64(numberOfSegmentToCreate)
for pointInserted := 1; pointInserted < numberOfSegmentToCreate; pointInserted++ {
allSegmentizedCoordinates = append(
allSegmentizedCoordinates,
b.X()*float64(pointInserted)*segmentFraction+a.X()*(1-float64(pointInserted)*segmentFraction),
b.Y()*float64(pointInserted)*segmentFraction+a.Y()*(1-float64(pointInserted)*segmentFraction),
)
}
return allSegmentizedCoordinates
}
// Copyright 2020 The Cockroach Authors.
//
// Use of this software is governed by the Business Source License
// included in the file licenses/BSL.txt.
//
// As of the Change Date specified in that file, in accordance with
// the Business Source License, use of this software will be governed
// by the Apache License, Version 2.0, included in the file
// licenses/APL.txt.
package geomfn
import (
"fmt"
"testing"
"github.com/cockroachdb/cockroach/pkg/geo"
"github.com/stretchr/testify/require"
"github.com/twpayne/go-geom"
)
func TestSegmentize(t *testing.T) {
segmentizeTestCases := []struct {
wkt string
maxSegmentLength float64
expectedWKT string
}{
{
wkt: "POINT (1.0 1.0)",
maxSegmentLength: 1,
expectedWKT: "POINT (1.0 1.0)",
},
{
wkt: "LINESTRING (1.0 1.0, 2.0 2.0, 3.0 3.0)",
maxSegmentLength: 1,
expectedWKT: "LINESTRING (1.0 1.0, 1.5 1.5, 2.0 2.0, 2.5 2.5, 3.0 3.0)",
},
{
wkt: "LINESTRING (1.0 1.0, 2.0 2.0, 3.0 3.0)",
maxSegmentLength: 0.33333,
expectedWKT: "LINESTRING (1.0 1.0, 1.2000000000000002 1.2000000000000002, 1.4 1.4, 1.6 1.6, 1.8 1.8, 2.0 2.0, 2.2 2.2, 2.4000000000000004 2.4000000000000004, 2.5999999999999996 2.5999999999999996, 2.8000000000000003 2.8000000000000003, 3 3)",
},
{
wkt: "LINESTRING EMPTY",
maxSegmentLength: 1,
expectedWKT: "LINESTRING EMPTY",
},
{
wkt: "LINESTRING (1.0 1.0, 2.0 2.0, 3.0 3.0)",
maxSegmentLength: 2,
expectedWKT: "LINESTRING (1.0 1.0, 2.0 2.0, 3.0 3.0)",
},
{
wkt: "LINESTRING (0.0 0.0, 0.0 10.0, 0.0 16.0)",
maxSegmentLength: 3,
expectedWKT: "LINESTRING (0.0 0.0,0.0 2.5,0.0 5.0,0.0 7.5,0.0 10.0,0.0 13.0,0.0 16.0)",
},
{
wkt: "POLYGON ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0))",
maxSegmentLength: 0.8,
expectedWKT: "POLYGON ((0.0 0.0, 0.5 0.0, 1.0 0.0, 1.0 0.5, 1.0 1.0, 0.5 0.5, 0.0 0.0))",
},
{
wkt: "POLYGON ((0.0 0.0, 1.0 0.0, 3.0 3.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))",
maxSegmentLength: 1,
expectedWKT: "POLYGON ((0.0 0.0, 1.0 0.0, 1.5 0.75, 2.0 1.5, 2.5 2.25, 3.0 3.0, 2.4000000000000004 2.4000000000000004, 1.7999999999999998 1.7999999999999998, 1.1999999999999997 1.1999999999999997, 0.5999999999999999 0.5999999999999999, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))",
},
{
wkt: "POLYGON ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))",
maxSegmentLength: 5,
expectedWKT: "POLYGON ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1))",
},
{
wkt: "POLYGON EMPTY",
maxSegmentLength: 1,
expectedWKT: "POLYGON EMPTY",
},
{
wkt: "MULTIPOINT ((1.0 1.0), (2.0 2.0))",
maxSegmentLength: 1,
expectedWKT: "MULTIPOINT ((1.0 1.0), (2.0 2.0))",
},
{
wkt: "MULTILINESTRING ((1.0 1.0, 2.0 2.0, 3.0 3.0), (6.0 6.0, 7.0 6.0))",
maxSegmentLength: 1,
expectedWKT: "MULTILINESTRING ((1.0 1.0, 1.5 1.5, 2.0 2.0, 2.5 2.5, 3.0 3.0), (6.0 6.0, 7.0 6.0))",
},
{
wkt: "MULTILINESTRING (EMPTY, (1.0 1.0, 2.0 2.0, 3.0 3.0), (6.0 6.0, 7.0 6.0))",
maxSegmentLength: 1,
expectedWKT: "MULTILINESTRING (EMPTY, (1.0 1.0, 1.5 1.5, 2.0 2.0, 2.5 2.5, 3.0 3.0), (6.0 6.0, 7.0 6.0))",
},
{
wkt: "MULTIPOLYGON (((3.0 3.0, 4.0 3.0, 4.0 4.0, 3.0 3.0)), ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1)))",
maxSegmentLength: 1,
expectedWKT: "MULTIPOLYGON (((3.0 3.0, 4.0 3.0, 4.0 4.0, 3.5 3.5, 3.0 3.0)), ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.5 0.5, 0.0 0.0), (0.1 0.1, 0.2 0.1, 0.2 0.2, 0.1 0.1)))",
},
{
wkt: "GEOMETRYCOLLECTION (POINT (40.0 10.0), LINESTRING (10.0 10.0, 20.0 20.0, 10.0 40.0), POLYGON ((40.0 40.0, 20.0 45.0, 45.0 30.0, 40.0 40.0)))",
maxSegmentLength: 10,
expectedWKT: "GEOMETRYCOLLECTION (POINT (40.0 10.0), LINESTRING (10.0 10.0, 15.0 15.0, 20.0 20.0, 16.666666666666668 26.666666666666668, 13.333333333333334 33.33333333333333, 10.0 40.0), POLYGON ((40.0 40.0, 33.333333333333336 41.66666666666667, 26.666666666666668 43.333333333333336, 20.0 45.0, 28.333333333333336 40.0, 36.66666666666667 35.0, 45.0 30.0, 42.5 35.0, 40.0 40.0)))",
},
{
wkt: "MULTIPOINT ((0.0 0.0), (1.0 1.0))",
maxSegmentLength: -1,
expectedWKT: "MULTIPOINT ((0.0 0.0), (1.0 1.0))",
},
}
for _, test := range segmentizeTestCases {
t.Run(fmt.Sprintf("%s, maximum segment length: %f", test.wkt, test.maxSegmentLength), func(t *testing.T) {
geom, err := geo.ParseGeometry(test.wkt)
require.NoError(t, err)
modifiedGeom, err := Segmentize(geom, test.maxSegmentLength)
require.NoError(t, err)
expectedGeom, err := geo.ParseGeometry(test.expectedWKT)
require.NoError(t, err)
require.Equal(t, expectedGeom, modifiedGeom)
})
}
// Test for segment maximum length as negative.
t.Run("Error when maximum segment length is less than 0", func(t *testing.T) {
geom, err := geo.ParseGeometry("MULTILINESTRING ((0 0, 1 1, 5 5), (5 5, 0 0))")
require.NoError(t, err)
_, err = Segmentize(geom, 0)
require.EqualError(t, err, "maximum segment length must be positive")
})
}
func TestSegmentizeCoords(t *testing.T) {
testCases := []struct {
desc string
a geom.Coord
b geom.Coord
segmentMaxLength float64
resultantCoordinates []float64
}{
{
desc: `Coordinate(0, 0) to Coordinate(1, 1), 1`,
a: geom.Coord{0, 0},
b: geom.Coord{1, 1},
segmentMaxLength: 1,
resultantCoordinates: []float64{0, 0, 0.5, 0.5},
},
{
desc: `Coordinate(0, 0) to Coordinate(1, 1), 0.3`,
a: geom.Coord{0, 0},
b: geom.Coord{1, 1},
segmentMaxLength: 0.3,
resultantCoordinates: []float64{0, 0, 0.2, 0.2, 0.4, 0.4, 0.6000000000000001, 0.6000000000000001, 0.8, 0.8},
},
{
desc: `Coordinate(0, 0) to Coordinate(1, 0), 0.49999999999999`,
a: geom.Coord{0, 0},
b: geom.Coord{1, 0},
segmentMaxLength: 0.49999999999999,
resultantCoordinates: []float64{0, 0, 0.3333333333333333, 0, 0.6666666666666666, 0},
},
{
desc: `Coordinate(1, 1) to Coordinate(0, 0), -1`,
a: geom.Coord{1, 1},
b: geom.Coord{0, 0},
segmentMaxLength: -1,
resultantCoordinates: []float64{1, 1},
},
{
desc: `Coordinate(1, 1) to Coordinate(0, 0), 2`,
a: geom.Coord{1, 1},
b: geom.Coord{0, 0},
segmentMaxLength: 2,
resultantCoordinates: []float64{1, 1},
},
{
desc: `Coordinate(0, 0) to Coordinate(0, 0), 1`,
a: geom.Coord{0, 0},
b: geom.Coord{0, 0},
segmentMaxLength: 1,
resultantCoordinates: []float64{0, 0},
},
}
for _, test := range testCases {
t.Run(test.desc, func(t *testing.T) {
convertedPoints := segmentizeCoords(test.a, test.b, test.segmentMaxLength)
require.Equal(t, test.resultantCoordinates, convertedPoints)
})
}
}
// Copyright 2020 The Cockroach Authors.
//
// Use of this software is governed by the Business Source License
// included in the file licenses/BSL.txt.
//
// As of the Change Date specified in that file, in accordance with
// the Business Source License, use of this software will be governed
// by the Apache License, Version 2.0, included in the file
// licenses/APL.txt.
package geosegmentize
import (
"github.com/cockroachdb/errors"
"github.com/twpayne/go-geom"
)
// SegmentizeGeom returns a modified geom.T having no segment longer
// than the given maximum segment length.
// segmentMaxAngleOrLength represents two different things depending
// on the object, which is about to segmentize as in case of geography
// it represents maximum segment angle whereas, in case of geometry it
// represents maximum segment distance.
// segmentizeCoords represents the function's definition which allows
// us to segmentize given two-points. We have to specify segmentizeCoords
// explicitly, as the algorithm for segmentization is significantly
// different for geometry and geography.
func SegmentizeGeom(
geometry geom.T,
segmentMaxAngleOrLength float64,
segmentizeCoords func(geom.Coord, geom.Coord, float64) []float64,
) (geom.T, error) {
if geometry.Empty() {
return geometry, nil
}
switch geometry := geometry.(type) {
case *geom.Point, *geom.MultiPoint:
return geometry, nil
case *geom.LineString:
var allFlatCoordinates []float64
for pointIdx := 1; pointIdx < geometry.NumCoords(); pointIdx++ {
allFlatCoordinates = append(
allFlatCoordinates,
segmentizeCoords(geometry.Coord(pointIdx-1), geometry.Coord(pointIdx), segmentMaxAngleOrLength)...,
)
}
// Appending end point as it wasn't included in the iteration of coordinates.
allFlatCoordinates = append(allFlatCoordinates, geometry.Coord(geometry.NumCoords()-1)...)
return geom.NewLineStringFlat(geom.XY, allFlatCoordinates).SetSRID(geometry.SRID()), nil
case *geom.MultiLineString:
segMultiLine := geom.NewMultiLineString(geom.XY).SetSRID(geometry.SRID())
for lineIdx := 0; lineIdx < geometry.NumLineStrings(); lineIdx++ {
l, err := SegmentizeGeom(geometry.LineString(lineIdx), segmentMaxAngleOrLength, segmentizeCoords)
if err != nil {
return nil, err
}
err = segMultiLine.Push(l.(*geom.LineString))
if err != nil {
return nil, err
}
}
return segMultiLine, nil
case *geom.LinearRing:
var allFlatCoordinates []float64
for pointIdx := 1; pointIdx < geometry.NumCoords(); pointIdx++ {
allFlatCoordinates = append(
allFlatCoordinates,
segmentizeCoords(geometry.Coord(pointIdx-1), geometry.Coord(pointIdx), segmentMaxAngleOrLength)...,
)
}
// Appending end point as it wasn't included in the iteration of coordinates.
allFlatCoordinates = append(allFlatCoordinates, geometry.Coord(geometry.NumCoords()-1)...)
return geom.NewLinearRingFlat(geom.XY, allFlatCoordinates).SetSRID(geometry.SRID()), nil
case *geom.Polygon:
segPolygon := geom.NewPolygon(geom.XY).SetSRID(geometry.SRID())
for loopIdx := 0; loopIdx < geometry.NumLinearRings(); loopIdx++ {
l, err := SegmentizeGeom(geometry.LinearRing(loopIdx), segmentMaxAngleOrLength, segmentizeCoords)
if err != nil {
return nil, err
}
err = segPolygon.Push(l.(*geom.LinearRing))
if err != nil {
return nil, err
}
}
return segPolygon, nil
case *geom.MultiPolygon:
segMultiPolygon := geom.NewMultiPolygon(geom.XY).SetSRID(geometry.SRID())
for polygonIdx := 0; polygonIdx < geometry.NumPolygons(); polygonIdx++ {
p, err := SegmentizeGeom(geometry.Polygon(polygonIdx), segmentMaxAngleOrLength, segmentizeCoords)
if err != nil {
return nil, err
}
err = segMultiPolygon.Push(p.(*geom.Polygon))
if err != nil {
return nil, err
}
}
return segMultiPolygon, nil
case *geom.GeometryCollection:
segGeomCollection := geom.NewGeometryCollection().SetSRID(geometry.SRID())
for geoIdx := 0; geoIdx < geometry.NumGeoms(); geoIdx++ {
g, err := SegmentizeGeom(geometry.Geom(geoIdx), segmentMaxAngleOrLength, segmentizeCoords)
if err != nil {
return nil, err
}
err = segGeomCollection.Push(g)
if err != nil {
return nil, err
}
}
return segGeomCollection, nil
}
return nil, errors.Newf("unknown type: %T", geometry)
}
......@@ -1434,6 +1434,33 @@ MULTIPOINT (0 0, 1 1)
statement error st_segmentize\(\): maximum segment length must be positive
SELECT ST_Segmentize('POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 1.0, 0.0 0.0))'::geography, 0)
query TTT
SELECT
dsc,
ST_AsText(ST_Segmentize(geom, 1)),
ST_AsText(ST_Segmentize(geom, 0.3))
FROM geom_operators_test
ORDER BY dsc
----
Empty GeometryCollection GEOMETRYCOLLECTION EMPTY GEOMETRYCOLLECTION EMPTY
Empty LineString LINESTRING EMPTY LINESTRING EMPTY
Faraway point POINT (5 5) POINT (5 5)
Line going through left and right square LINESTRING (-0.5 0.5, 0.5 0.5) LINESTRING (-0.5 0.5, -0.25 0.5, 0 0.5, 0.25 0.5, 0.5 0.5)
NULL NULL NULL
Point middle of Left Square POINT (-0.5 0.5) POINT (-0.5 0.5)
Point middle of Right Square POINT (0.5 0.5) POINT (0.5 0.5)
Square (left) POLYGON ((-1 0, 0 0, 0 1, -1 1, -1 0)) POLYGON ((-1 0, -0.75 0, -0.5 0, -0.25 0, 0 0, 0 0.25, 0 0.5, 0 0.75, 0 1, -0.25 1, -0.5 1, -0.75 1, -1 1, -1 0.75, -1 0.5, -1 0.25, -1 0))
Square (right) POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0)) POLYGON ((0 0, 0.25 0, 0.5 0, 0.75 0, 1 0, 1 0.25, 1 0.5, 1 0.75, 1 1, 0.75 1, 0.5 1, 0.25 1, 0 1, 0 0.75, 0 0.5, 0 0.25, 0 0))
Square overlapping left and right square POLYGON ((-0.1 0, 0.45 0, 1 0, 1 1, 0.45 1, -0.1 1, -0.1 0)) POLYGON ((-0.1 0, 0.175 0, 0.45 0, 0.725 0, 1 0, 1 0.25, 1 0.5, 1 0.75, 1 1, 0.725 1, 0.45 1, 0.175 1, -0.1 1, -0.1 0.75, -0.1 0.5, -0.1 0.25, -0.1 0))
query T
SELECT ST_AsText(ST_Segmentize('MULTIPOINT (0 0, 1 1)'::geometry, -1))
----
MULTIPOINT (0 0, 1 1)
statement error st_segmentize\(\): maximum segment length must be positive
SELECT ST_Segmentize('POLYGON ((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 1.0, 0.0 0.0))'::geometry, -1)
subtest pg_extension
statement ok
......
......@@ -1995,11 +1995,11 @@ Note If the result has zero or one points, it will be returned as a POINT. If it
Fn: func(_ *tree.EvalContext, args tree.Datums) (tree.Datum, error) {
g := args[0].(*tree.DGeography)
segmentMaxLength := float64(*args[1].(*tree.DFloat))
segGeometry, err := geogfn.Segmentize(g.Geography, segmentMaxLength)
segGeography, err := geogfn.Segmentize(g.Geography, segmentMaxLength)
if err != nil {
return nil, err
}
return tree.NewDGeography(segGeometry), nil
return tree.NewDGeography(segGeography), nil
},
Info: infoBuilder{
info: `Returns a modified Geography having no segment longer than the given max_segment_length meters.
......@@ -2009,6 +2009,27 @@ The calculations are done on a sphere.`,
}.String(),
Volatility: tree.VolatilityImmutable,
},
tree.Overload{
Types: tree.ArgTypes{
{"geometry", types.Geometry},
{"max_segment_length", types.Float},
},
ReturnType: tree.FixedReturnType(types.Geometry),
Fn: func(_ *tree.EvalContext, args tree.Datums) (tree.Datum, error) {
g := args[0].(*tree.DGeometry)
segmentMaxLength := float64(*args[1].(*tree.DFloat))
segGeometry, err := geomfn.Segmentize(g.Geometry, segmentMaxLength)