st_interiorringn(geometry: geometry, n: int) → geometryReturns the n-th (1-indexed) interior ring of a Polygon as a LineString. Returns NULL if the shape is not a Polygon, or the ring does not exist.
st_intersection(geometry_a: geometry, geometry_b: geometry) → geometryReturns the point intersections of the given geometries.
+This function utilizes the GEOS module.
+st_intersects(geography_a: geography, geography_b: geography) → boolReturns true if geography_a shares any portion of space with geography_b.
The calculations performed are have a precision of 1cm.
This function utilizes the S2 library for spherical calculations.
@@ -1088,6 +1091,9 @@ given Geometry.st_pointn(geometry: geometry, n: int) → geometryReturns the n-th Point of a LineString (1-indexed). Returns NULL if out of bounds or not a LineString.
st_pointonsurface(geometry: geometry) → geometryReturns a point that intersects with the given Geometry.
+This function utilizes the GEOS module.
+st_polyfromtext(str: string, srid: int) → geometryReturns the Geometry from a WKT or EWKT representation with an SRID. If the shape underneath is not Polygon, NULL is returned. If the SRID is present in both the EWKT and the argument, the argument value is used.
st_polyfromtext(val: string) → geometryReturns the Geometry from a WKT or EWKT representation. If the shape underneath is not Polygon, NULL is returned.
@@ -1160,6 +1166,9 @@ given Geometry.st_transform(geometry: geometry, to_proj_text: string) → geometryTransforms a geometry into the coordinate reference system referenced by the projection text by projecting its coordinates.
This function utilizes the PROJ library for coordinate projections.
st_union(geometry_a: geometry, geometry_b: geometry) → geometryReturns the union of the given geometries as a single Geometry object.
+This function utilizes the GEOS module.
+st_within(geometry_a: geometry, geometry_b: geometry) → boolReturns true if geometry_a is completely inside geometry_b.
This function utilizes the GEOS module.
This function will automatically use any available index.
diff --git a/pkg/geo/geo.go b/pkg/geo/geo.go index d835c787eac7..9bd731aa9d17 100644 --- a/pkg/geo/geo.go +++ b/pkg/geo/geo.go @@ -241,6 +241,12 @@ func (g *Geometry) BoundingBoxIntersects(o *Geometry) bool { return g.SpatialObject.BoundingBox.Intersects(o.SpatialObject.BoundingBox) } +// Layout returns the geom layout of the given geometry. +func (g *Geometry) Layout() geom.Layout { + // We are currently always 2D. + return geom.XY +} + // // Geography // diff --git a/pkg/geo/geomfn/topology_operations.go b/pkg/geo/geomfn/topology_operations.go new file mode 100644 index 000000000000..4d80056345be --- /dev/null +++ b/pkg/geo/geomfn/topology_operations.go @@ -0,0 +1,65 @@ +// 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 ( + "github.com/cockroachdb/cockroach/pkg/geo" + "github.com/cockroachdb/cockroach/pkg/geo/geos" + "github.com/twpayne/go-geom" +) + +// Centroid returns the Centroid of a given Geometry. +func Centroid(g *geo.Geometry) (*geo.Geometry, error) { + if g.Empty() { + return geo.NewGeometryFromGeom(geom.NewPointEmpty(g.Layout())) + } + centroidEWKB, err := geos.Centroid(g.EWKB()) + if err != nil { + return nil, err + } + return geo.ParseGeometryFromEWKB(centroidEWKB) +} + +// PointOnSurface returns the PointOnSurface of a given Geometry. +func PointOnSurface(g *geo.Geometry) (*geo.Geometry, error) { + if g.Empty() { + return geo.NewGeometryFromGeom(geom.NewPointEmpty(g.Layout())) + } + pointOnSurfaceEWKB, err := geos.PointOnSurface(g.EWKB()) + if err != nil { + return nil, err + } + return geo.ParseGeometryFromEWKB(pointOnSurfaceEWKB) +} + +// Intersection returns the geometries of intersection between A and B. +func Intersection(a *geo.Geometry, b *geo.Geometry) (*geo.Geometry, error) { + if a.SRID() != b.SRID() { + return nil, geo.NewMismatchingSRIDsError(a, b) + } + retEWKB, err := geos.Intersection(a.EWKB(), b.EWKB()) + if err != nil { + return nil, err + } + return geo.ParseGeometryFromEWKB(retEWKB) +} + +// Union returns the geometries of intersection between A and B. +func Union(a *geo.Geometry, b *geo.Geometry) (*geo.Geometry, error) { + if a.SRID() != b.SRID() { + return nil, geo.NewMismatchingSRIDsError(a, b) + } + retEWKB, err := geos.Union(a.EWKB(), b.EWKB()) + if err != nil { + return nil, err + } + return geo.ParseGeometryFromEWKB(retEWKB) +} diff --git a/pkg/geo/geomfn/topology_operations_test.go b/pkg/geo/geomfn/topology_operations_test.go new file mode 100644 index 000000000000..832685a05f84 --- /dev/null +++ b/pkg/geo/geomfn/topology_operations_test.go @@ -0,0 +1,149 @@ +// 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 TestCentroid(t *testing.T) { + testCases := []struct { + wkt string + expected string + }{ + {"POINT(1.0 1.0)", "POINT (1.0 1.0)"}, + {"SRID=4326;POINT(1.0 1.0)", "SRID=4326;POINT (1.0 1.0)"}, + {"LINESTRING(1.0 1.0, 2.0 2.0, 3.0 3.0)", "POINT (2.0 2.0)"}, + {"POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0))", "POINT (0.666666666666667 0.333333333333333)"}, + {"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))", "POINT (0.671717171717172 0.335353535353535)"}, + {"MULTIPOINT((1.0 1.0), (2.0 2.0))", "POINT (1.5 1.5)"}, + {"MULTILINESTRING((1.0 1.0, 2.0 2.0, 3.0 3.0), (6.0 6.0, 7.0 6.0))", "POINT (3.17541743733684 3.04481549985497)"}, + {"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)))", "POINT (2.17671691792295 1.84187604690117)"}, + {"GEOMETRYCOLLECTION (POINT (40 10),LINESTRING (10 10, 20 20, 10 40),POLYGON ((40 40, 20 45, 45 30, 40 40)))", "POINT (35 38.3333333333333)"}, + } + + for _, tc := range testCases { + t.Run(tc.wkt, func(t *testing.T) { + g, err := geo.ParseGeometry(tc.wkt) + require.NoError(t, err) + ret, err := Centroid(g) + require.NoError(t, err) + + retAsGeomT, err := ret.AsGeomT() + require.NoError(t, err) + + expected, err := geo.ParseGeometry(tc.expected) + require.NoError(t, err) + expectedAsGeomT, err := expected.AsGeomT() + require.NoError(t, err) + + // Ensure points are close in terms of precision. + require.InEpsilon(t, expectedAsGeomT.(*geom.Point).X(), retAsGeomT.(*geom.Point).X(), 2e-10) + require.InEpsilon(t, expectedAsGeomT.(*geom.Point).Y(), retAsGeomT.(*geom.Point).Y(), 2e-10) + require.Equal(t, expected.SRID(), ret.SRID()) + }) + } +} + +func TestPointOnSurface(t *testing.T) { + testCases := []struct { + wkt string + expected string + }{ + {"POINT(1.0 1.0)", "POINT (1.0 1.0)"}, + {"SRID=4326;POINT(1.0 1.0)", "SRID=4326;POINT (1.0 1.0)"}, + {"LINESTRING(1.0 1.0, 2.0 2.0, 3.0 3.0)", "POINT (2.0 2.0)"}, + {"POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0))", "POINT(0.75 0.5)"}, + {"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))", "POINT(0.8 0.6)"}, + {"MULTIPOINT((1.0 1.0), (2.0 2.0))", "POINT (1 1)"}, + {"MULTILINESTRING((1.0 1.0, 2.0 2.0, 3.0 3.0), (6.0 6.0, 7.0 6.0))", "POINT (2 2)"}, + {"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)))", "POINT(3.75 3.5)"}, + {"GEOMETRYCOLLECTION (POINT (40 10),LINESTRING (10 10, 20 20, 10 40),POLYGON ((40 40, 20 45, 45 30, 40 40)))", "POINT(39.5833333333333 35)"}, + } + + for _, tc := range testCases { + t.Run(tc.wkt, func(t *testing.T) { + g, err := geo.ParseGeometry(tc.wkt) + require.NoError(t, err) + ret, err := PointOnSurface(g) + require.NoError(t, err) + + retAsGeomT, err := ret.AsGeomT() + require.NoError(t, err) + + expected, err := geo.ParseGeometry(tc.expected) + require.NoError(t, err) + expectedAsGeomT, err := expected.AsGeomT() + require.NoError(t, err) + + // Ensure points are close in terms of precision. + require.InEpsilon(t, expectedAsGeomT.(*geom.Point).X(), retAsGeomT.(*geom.Point).X(), 2e-10) + require.InEpsilon(t, expectedAsGeomT.(*geom.Point).Y(), retAsGeomT.(*geom.Point).Y(), 2e-10) + require.Equal(t, expected.SRID(), ret.SRID()) + }) + } +} + +func TestIntersection(t *testing.T) { + testCases := []struct { + a *geo.Geometry + b *geo.Geometry + expected *geo.Geometry + }{ + {rightRect, rightRect, geo.MustParseGeometry("POLYGON ((1 0, 0 0, 0 1, 1 1, 1 0))")}, + {rightRect, rightRectPoint, rightRectPoint}, + {rightRectPoint, rightRectPoint, rightRectPoint}, + } + + for i, tc := range testCases { + t.Run(fmt.Sprintf("tc:%d", i), func(t *testing.T) { + g, err := Intersection(tc.a, tc.b) + require.NoError(t, err) + require.Equal(t, tc.expected, g) + }) + } + + t.Run("errors if SRIDs mismatch", func(t *testing.T) { + _, err := Intersection(mismatchingSRIDGeometryA, mismatchingSRIDGeometryB) + requireMismatchingSRIDError(t, err) + }) +} + +func TestUnion(t *testing.T) { + testCases := []struct { + a *geo.Geometry + b *geo.Geometry + expected *geo.Geometry + }{ + {rightRect, rightRect, geo.MustParseGeometry("POLYGON ((1 0, 0 0, 0 1, 1 1, 1 0))")}, + {rightRect, rightRectPoint, geo.MustParseGeometry("POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))")}, + {rightRectPoint, rightRectPoint, rightRectPoint}, + {leftRect, rightRect, geo.MustParseGeometry("POLYGON ((0 0, -1 0, -1 1, 0 1, 1 1, 1 0, 0 0))")}, + } + + for i, tc := range testCases { + t.Run(fmt.Sprintf("tc:%d", i), func(t *testing.T) { + g, err := Union(tc.a, tc.b) + require.NoError(t, err) + require.Equal(t, tc.expected, g) + }) + } + + t.Run("errors if SRIDs mismatch", func(t *testing.T) { + _, err := Union(mismatchingSRIDGeometryA, mismatchingSRIDGeometryB) + requireMismatchingSRIDError(t, err) + }) +} diff --git a/pkg/geo/geomfn/unary_operators.go b/pkg/geo/geomfn/unary_operators.go index 94d23c720898..1219966eb4bf 100644 --- a/pkg/geo/geomfn/unary_operators.go +++ b/pkg/geo/geomfn/unary_operators.go @@ -19,25 +19,6 @@ import ( "github.com/twpayne/go-geom/encoding/ewkb" ) -// Centroid returns the Centroid of a given Geometry. -func Centroid(g *geo.Geometry) (*geo.Geometry, error) { - // Empty geometries do not react well in GEOS, so we have to - // convert and check beforehand. - // Remove after #49209 is resolved. - t, err := g.AsGeomT() - if err != nil { - return nil, err - } - if t.Empty() { - return geo.NewGeometryFromGeom(geom.NewPointEmpty(geom.XY)) - } - centroidEWKB, err := geos.Centroid(g.EWKB()) - if err != nil { - return nil, err - } - return geo.ParseGeometryFromEWKB(centroidEWKB) -} - // Length returns the length of a given Geometry. // Note only (MULTI)LINESTRING objects have a length. // (MULTI)POLYGON objects should use Perimeter. diff --git a/pkg/geo/geomfn/unary_operators_test.go b/pkg/geo/geomfn/unary_operators_test.go index 1ee9903e6857..9a837f0b53fa 100644 --- a/pkg/geo/geomfn/unary_operators_test.go +++ b/pkg/geo/geomfn/unary_operators_test.go @@ -15,48 +15,8 @@ import ( "github.com/cockroachdb/cockroach/pkg/geo" "github.com/stretchr/testify/require" - "github.com/twpayne/go-geom" ) -func TestCentroid(t *testing.T) { - testCases := []struct { - wkt string - expected string - }{ - {"POINT(1.0 1.0)", "POINT (1.0 1.0)"}, - {"SRID=4326;POINT(1.0 1.0)", "SRID=4326;POINT (1.0 1.0)"}, - {"LINESTRING(1.0 1.0, 2.0 2.0, 3.0 3.0)", "POINT (2.0 2.0)"}, - {"POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 0.0))", "POINT (0.666666666666667 0.333333333333333)"}, - {"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))", "POINT (0.671717171717172 0.335353535353535)"}, - {"MULTIPOINT((1.0 1.0), (2.0 2.0))", "POINT (1.5 1.5)"}, - {"MULTILINESTRING((1.0 1.0, 2.0 2.0, 3.0 3.0), (6.0 6.0, 7.0 6.0))", "POINT (3.17541743733684 3.04481549985497)"}, - {"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)))", "POINT (2.17671691792295 1.84187604690117)"}, - {"GEOMETRYCOLLECTION (POINT (40 10),LINESTRING (10 10, 20 20, 10 40),POLYGON ((40 40, 20 45, 45 30, 40 40)))", "POINT (35 38.3333333333333)"}, - } - - for _, tc := range testCases { - t.Run(tc.wkt, func(t *testing.T) { - g, err := geo.ParseGeometry(tc.wkt) - require.NoError(t, err) - ret, err := Centroid(g) - require.NoError(t, err) - - retAsGeomT, err := ret.AsGeomT() - require.NoError(t, err) - - expected, err := geo.ParseGeometry(tc.expected) - require.NoError(t, err) - expectedAsGeomT, err := expected.AsGeomT() - require.NoError(t, err) - - // Ensure points are close in terms of precision. - require.InEpsilon(t, expectedAsGeomT.(*geom.Point).X(), retAsGeomT.(*geom.Point).X(), 2e-10) - require.InEpsilon(t, expectedAsGeomT.(*geom.Point).Y(), retAsGeomT.(*geom.Point).Y(), 2e-10) - require.Equal(t, expected.SRID(), ret.SRID()) - }) - } -} - func TestLength(t *testing.T) { testCases := []struct { wkt string diff --git a/pkg/geo/geos/geos.cc b/pkg/geo/geos/geos.cc index 8c2180626596..6cfb48427115 100644 --- a/pkg/geo/geos/geos.cc +++ b/pkg/geo/geos/geos.cc @@ -81,7 +81,12 @@ typedef CR_GEOS_Geometry (*CR_GEOS_BufferWithParams_r)(CR_GEOS_Handle, CR_GEOS_G typedef int (*CR_GEOS_Area_r)(CR_GEOS_Handle, CR_GEOS_Geometry, double*); typedef int (*CR_GEOS_Length_r)(CR_GEOS_Handle, CR_GEOS_Geometry, double*); + typedef CR_GEOS_Geometry (*CR_GEOS_Centroid_r)(CR_GEOS_Handle, CR_GEOS_Geometry); +typedef CR_GEOS_Geometry (*CR_GEOS_Union_r)(CR_GEOS_Handle, CR_GEOS_Geometry, CR_GEOS_Geometry); +typedef CR_GEOS_Geometry (*CR_GEOS_Intersection_r)(CR_GEOS_Handle, CR_GEOS_Geometry, + CR_GEOS_Geometry); +typedef CR_GEOS_Geometry (*CR_GEOS_PointOnSurface_r)(CR_GEOS_Handle, CR_GEOS_Geometry); typedef CR_GEOS_Geometry (*CR_GEOS_Interpolate_r)(CR_GEOS_Handle, CR_GEOS_Geometry, double); @@ -147,7 +152,11 @@ struct CR_GEOS { CR_GEOS_Area_r GEOSArea_r; CR_GEOS_Length_r GEOSLength_r; + CR_GEOS_Centroid_r GEOSGetCentroid_r; + CR_GEOS_Union_r GEOSUnion_r; + CR_GEOS_PointOnSurface_r GEOSPointOnSurface_r; + CR_GEOS_Intersection_r GEOSIntersection_r; CR_GEOS_Interpolate_r GEOSInterpolate_r; @@ -209,6 +218,9 @@ struct CR_GEOS { INIT(GEOSArea_r); INIT(GEOSLength_r); INIT(GEOSGetCentroid_r); + INIT(GEOSUnion_r); + INIT(GEOSPointOnSurface_r); + INIT(GEOSIntersection_r); INIT(GEOSInterpolate_r); INIT(GEOSDistance_r); INIT(GEOSCovers_r); @@ -442,6 +454,10 @@ CR_GEOS_Status CR_GEOS_Length(CR_GEOS* lib, CR_GEOS_Slice a, double* ret) { return CR_GEOS_UnaryOperator(lib, lib->GEOSLength_r, a, ret); } +// +// Topology operators. +// + CR_GEOS_Status CR_GEOS_Centroid(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String* centroidEWKB) { std::string error; auto handle = initHandleWithErrorBuffer(lib, &error); @@ -460,34 +476,107 @@ CR_GEOS_Status CR_GEOS_Centroid(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String* c return toGEOSString(error.data(), error.length()); } +CR_GEOS_Status CR_GEOS_Union(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, + CR_GEOS_String* unionEWKB) { + std::string error; + auto handle = initHandleWithErrorBuffer(lib, &error); + *unionEWKB = {.data = NULL, .len = 0}; + + auto geomA = CR_GEOS_GeometryFromSlice(lib, handle, a); + auto geomB = CR_GEOS_GeometryFromSlice(lib, handle, b); + if (geomA != nullptr && geomB != nullptr) { + auto unionGeom = lib->GEOSUnion_r(handle, geomA, geomB); + if (unionGeom != nullptr) { + auto srid = lib->GEOSGetSRID_r(handle, geomA); + CR_GEOS_writeGeomToEWKB(lib, handle, unionGeom, unionEWKB, srid); + lib->GEOSGeom_destroy_r(handle, unionGeom); + } + } + if (geomA != nullptr) { + lib->GEOSGeom_destroy_r(handle, geomA); + } + if (geomB != nullptr) { + lib->GEOSGeom_destroy_r(handle, geomB); + } + + lib->GEOS_finish_r(handle); + return toGEOSString(error.data(), error.length()); +} + +CR_GEOS_Status CR_GEOS_PointOnSurface(CR_GEOS* lib, CR_GEOS_Slice a, + CR_GEOS_String* pointOnSurfaceEWKB) { + std::string error; + auto handle = initHandleWithErrorBuffer(lib, &error); + auto geom = CR_GEOS_GeometryFromSlice(lib, handle, a); + *pointOnSurfaceEWKB = {.data = NULL, .len = 0}; + if (geom != nullptr) { + auto pointOnSurfaceGeom = lib->GEOSPointOnSurface_r(handle, geom); + if (pointOnSurfaceGeom != nullptr) { + auto srid = lib->GEOSGetSRID_r(handle, geom); + CR_GEOS_writeGeomToEWKB(lib, handle, pointOnSurfaceGeom, pointOnSurfaceEWKB, srid); + lib->GEOSGeom_destroy_r(handle, pointOnSurfaceGeom); + } + lib->GEOSGeom_destroy_r(handle, geom); + } + lib->GEOS_finish_r(handle); + return toGEOSString(error.data(), error.length()); +} + +CR_GEOS_Status CR_GEOS_Intersection(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, + CR_GEOS_String* intersectionEWKB) { + std::string error; + auto handle = initHandleWithErrorBuffer(lib, &error); + *intersectionEWKB = {.data = NULL, .len = 0}; + + auto geomA = CR_GEOS_GeometryFromSlice(lib, handle, a); + auto geomB = CR_GEOS_GeometryFromSlice(lib, handle, b); + if (geomA != nullptr && geomB != nullptr) { + auto intersectionGeom = lib->GEOSIntersection_r(handle, geomA, geomB); + if (intersectionGeom != nullptr) { + auto srid = lib->GEOSGetSRID_r(handle, geomA); + CR_GEOS_writeGeomToEWKB(lib, handle, intersectionGeom, intersectionEWKB, srid); + lib->GEOSGeom_destroy_r(handle, intersectionGeom); + } + } + if (geomA != nullptr) { + lib->GEOSGeom_destroy_r(handle, geomA); + } + if (geomB != nullptr) { + lib->GEOSGeom_destroy_r(handle, geomB); + } + + lib->GEOS_finish_r(handle); + return toGEOSString(error.data(), error.length()); +} + // // Linear Reference // CR_GEOS_Status CR_GEOS_Interpolate(CR_GEOS* lib, CR_GEOS_Slice a, double distance, CR_GEOS_String* interpolatedPointEWKB) { - std::string error; - auto handle = initHandleWithErrorBuffer(lib, &error); - auto geom = CR_GEOS_GeometryFromSlice(lib, handle, a); - *interpolatedPointEWKB = {.data = NULL, .len = 0}; - if (geom != nullptr) { - auto interpolatedPoint = lib->GEOSInterpolate_r(handle, geom, distance); - if (interpolatedPoint != nullptr) { - auto srid = lib->GEOSGetSRID_r(handle, geom); - CR_GEOS_writeGeomToEWKB(lib, handle, interpolatedPoint, interpolatedPointEWKB, srid); - lib->GEOSGeom_destroy_r(handle, interpolatedPoint); - } - lib->GEOSGeom_destroy_r(handle, geom); - } - lib->GEOS_finish_r(handle); - return toGEOSString(error.data(), error.length()); + std::string error; + auto handle = initHandleWithErrorBuffer(lib, &error); + auto geom = CR_GEOS_GeometryFromSlice(lib, handle, a); + *interpolatedPointEWKB = {.data = NULL, .len = 0}; + if (geom != nullptr) { + auto interpolatedPoint = lib->GEOSInterpolate_r(handle, geom, distance); + if (interpolatedPoint != nullptr) { + auto srid = lib->GEOSGetSRID_r(handle, geom); + CR_GEOS_writeGeomToEWKB(lib, handle, interpolatedPoint, interpolatedPointEWKB, srid); + lib->GEOSGeom_destroy_r(handle, interpolatedPoint); + } + lib->GEOSGeom_destroy_r(handle, geom); + } + lib->GEOS_finish_r(handle); + return toGEOSString(error.data(), error.length()); } // // Binary operators // -CR_GEOS_Status CR_GEOS_Distance(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, double *ret) { +CR_GEOS_Status CR_GEOS_Distance(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, double* ret) { return CR_GEOS_BinaryOperator(lib, lib->GEOSDistance_r, a, b, ret); } diff --git a/pkg/geo/geos/geos.go b/pkg/geo/geos/geos.go index 2b0d97a53d97..a48ccd0812d9 100644 --- a/pkg/geo/geos/geos.go +++ b/pkg/geo/geos/geos.go @@ -320,6 +320,45 @@ func Centroid(ewkb geopb.EWKB) (geopb.EWKB, error) { return cStringToSafeGoBytes(cEWKB), nil } +// PointOnSurface returns an EWKB with a point that is on the surface of the given EWKB. +func PointOnSurface(ewkb geopb.EWKB) (geopb.EWKB, error) { + g, err := ensureInitInternal() + if err != nil { + return nil, err + } + var cEWKB C.CR_GEOS_String + if err := statusToError(C.CR_GEOS_PointOnSurface(g, goToCSlice(ewkb), &cEWKB)); err != nil { + return nil, err + } + return cStringToSafeGoBytes(cEWKB), nil +} + +// Intersection returns an EWKB which contains the geometries of intersection between A and B. +func Intersection(a geopb.EWKB, b geopb.EWKB) (geopb.EWKB, error) { + g, err := ensureInitInternal() + if err != nil { + return nil, err + } + var cEWKB C.CR_GEOS_String + if err := statusToError(C.CR_GEOS_Intersection(g, goToCSlice(a), goToCSlice(b), &cEWKB)); err != nil { + return nil, err + } + return cStringToSafeGoBytes(cEWKB), nil +} + +// Union returns an EWKB which is a union of shapes A and B. +func Union(a geopb.EWKB, b geopb.EWKB) (geopb.EWKB, error) { + g, err := ensureInitInternal() + if err != nil { + return nil, err + } + var cEWKB C.CR_GEOS_String + if err := statusToError(C.CR_GEOS_Union(g, goToCSlice(a), goToCSlice(b), &cEWKB)); err != nil { + return nil, err + } + return cStringToSafeGoBytes(cEWKB), nil +} + // InterpolateLine returns the point along the given LineString which is at // a given distance from starting point. // Note: For distance less than 0 it returns start point similarly for distance diff --git a/pkg/geo/geos/geos.h b/pkg/geo/geos/geos.h index 4ee3f5fa1634..f971d722a2dc 100644 --- a/pkg/geo/geos/geos.h +++ b/pkg/geo/geos/geos.h @@ -75,13 +75,24 @@ CR_GEOS_Status CR_GEOS_Buffer(CR_GEOS* lib, CR_GEOS_Slice ewkb, CR_GEOS_BufferPa CR_GEOS_Status CR_GEOS_Area(CR_GEOS* lib, CR_GEOS_Slice a, double* ret); CR_GEOS_Status CR_GEOS_Length(CR_GEOS* lib, CR_GEOS_Slice a, double* ret); + +// +// Topology operators. +// + CR_GEOS_Status CR_GEOS_Centroid(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_String* centroidEWKB); +CR_GEOS_Status CR_GEOS_Union(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, + CR_GEOS_String* unionEWKB); +CR_GEOS_Status CR_GEOS_PointOnSurface(CR_GEOS* lib, CR_GEOS_Slice a, + CR_GEOS_String* pointOnSurfaceEWKB); +CR_GEOS_Status CR_GEOS_Intersection(CR_GEOS* lib, CR_GEOS_Slice a, CR_GEOS_Slice b, + CR_GEOS_String* intersectionEWKB); // // Linear reference. // -CR_GEOS_Status CR_GEOS_Interpolate(CR_GEOS *lib, CR_GEOS_Slice a, double distance, - CR_GEOS_String *interpolatedPoint); +CR_GEOS_Status CR_GEOS_Interpolate(CR_GEOS* lib, CR_GEOS_Slice a, double distance, + CR_GEOS_String* interpolatedPoint); // // Binary operators. diff --git a/pkg/sql/logictest/testdata/logic_test/geospatial b/pkg/sql/logictest/testdata/logic_test/geospatial index 49f43e750064..62af0b82da95 100644 --- a/pkg/sql/logictest/testdata/logic_test/geospatial +++ b/pkg/sql/logictest/testdata/logic_test/geospatial @@ -376,23 +376,113 @@ Square (left) 1 1 0 4 Square (right) 1 1 0 4 Square overlapping left and right square 1.1 1.1 0 4.2 -query TT +# Topology operators. + +query TTT +SELECT + a.dsc, + b.dsc, + ST_AsEWKT(ST_Union(a.geom, b.geom)) +FROM geom_operators_test a +JOIN geom_operators_test b ON (1=1) +WHERE a.dsc NOT LIKE '%Empty%' and b.dsc NOT LIKE '%Empty%' -- Remove after #49209 is resolved. +ORDER BY a.dsc, b.dsc +---- +Faraway point Faraway point POINT (5 5) +Faraway point Line going through left and right square GEOMETRYCOLLECTION (POINT (5 5), LINESTRING (-0.5 0.5, 0.5 0.5)) +Faraway point NULL NULL +Faraway point Point middle of Left Square MULTIPOINT (5 5, -0.5 0.5) +Faraway point Point middle of Right Square MULTIPOINT (5 5, 0.5 0.5) +Faraway point Square (left) GEOMETRYCOLLECTION (POINT (5 5), POLYGON ((-1 0, 0 0, 0 1, -1 1, -1 0))) +Faraway point Square (right) GEOMETRYCOLLECTION (POINT (5 5), POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))) +Faraway point Square overlapping left and right square GEOMETRYCOLLECTION (POINT (5 5), POLYGON ((-0.1 0, 1 0, 1 1, -0.1 1, -0.1 0))) +Line going through left and right square Faraway point GEOMETRYCOLLECTION (LINESTRING (-0.5 0.5, 0.5 0.5), POINT (5 5)) +Line going through left and right square Line going through left and right square LINESTRING (-0.5 0.5, 0.5 0.5) +Line going through left and right square NULL NULL +Line going through left and right square Point middle of Left Square LINESTRING (-0.5 0.5, 0.5 0.5) +Line going through left and right square Point middle of Right Square LINESTRING (-0.5 0.5, 0.5 0.5) +Line going through left and right square Square (left) GEOMETRYCOLLECTION (LINESTRING (0 0.5, 0.5 0.5), POLYGON ((0 0.5, 0 0, -1 0, -1 1, 0 1, 0 0.5))) +Line going through left and right square Square (right) GEOMETRYCOLLECTION (LINESTRING (-0.5 0.5, 0 0.5), POLYGON ((0 0.5, 0 1, 1 1, 1 0, 0 0, 0 0.5))) +Line going through left and right square Square overlapping left and right square GEOMETRYCOLLECTION (LINESTRING (-0.5 0.5, -0.1 0.5), POLYGON ((-0.1 0.5, -0.1 1, 1 1, 1 0, -0.1 0, -0.1 0.5))) +NULL Faraway point NULL +NULL Line going through left and right square NULL +NULL NULL NULL +NULL Point middle of Left Square NULL +NULL Point middle of Right Square NULL +NULL Square (left) NULL +NULL Square (right) NULL +NULL Square overlapping left and right square NULL +Point middle of Left Square Faraway point MULTIPOINT (-0.5 0.5, 5 5) +Point middle of Left Square Line going through left and right square LINESTRING (-0.5 0.5, 0.5 0.5) +Point middle of Left Square NULL NULL +Point middle of Left Square Point middle of Left Square POINT (-0.5 0.5) +Point middle of Left Square Point middle of Right Square MULTIPOINT (-0.5 0.5, 0.5 0.5) +Point middle of Left Square Square (left) POLYGON ((-1 0, -1 1, 0 1, 0 0, -1 0)) +Point middle of Left Square Square (right) GEOMETRYCOLLECTION (POINT (-0.5 0.5), POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))) +Point middle of Left Square Square overlapping left and right square GEOMETRYCOLLECTION (POINT (-0.5 0.5), POLYGON ((-0.1 0, 1 0, 1 1, -0.1 1, -0.1 0))) +Point middle of Right Square Faraway point MULTIPOINT (0.5 0.5, 5 5) +Point middle of Right Square Line going through left and right square LINESTRING (-0.5 0.5, 0.5 0.5) +Point middle of Right Square NULL NULL +Point middle of Right Square Point middle of Left Square MULTIPOINT (0.5 0.5, -0.5 0.5) +Point middle of Right Square Point middle of Right Square POINT (0.5 0.5) +Point middle of Right Square Square (left) GEOMETRYCOLLECTION (POINT (0.5 0.5), POLYGON ((-1 0, 0 0, 0 1, -1 1, -1 0))) +Point middle of Right Square Square (right) POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0)) +Point middle of Right Square Square overlapping left and right square POLYGON ((-0.1 0, -0.1 1, 1 1, 1 0, -0.1 0)) +Square (left) Faraway point GEOMETRYCOLLECTION (POLYGON ((-1 0, 0 0, 0 1, -1 1, -1 0)), POINT (5 5)) +Square (left) Line going through left and right square GEOMETRYCOLLECTION (LINESTRING (0 0.5, 0.5 0.5), POLYGON ((0 0.5, 0 0, -1 0, -1 1, 0 1, 0 0.5))) +Square (left) NULL NULL +Square (left) Point middle of Left Square POLYGON ((-1 0, -1 1, 0 1, 0 0, -1 0)) +Square (left) Point middle of Right Square GEOMETRYCOLLECTION (POLYGON ((-1 0, 0 0, 0 1, -1 1, -1 0)), POINT (0.5 0.5)) +Square (left) Square (left) POLYGON ((0 0, -1 0, -1 1, 0 1, 0 0)) +Square (left) Square (right) POLYGON ((0 0, -1 0, -1 1, 0 1, 1 1, 1 0, 0 0)) +Square (left) Square overlapping left and right square POLYGON ((-0.1 0, -1 0, -1 1, -0.1 1, 0 1, 1 1, 1 0, 0 0, -0.1 0)) +Square (right) Faraway point GEOMETRYCOLLECTION (POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0)), POINT (5 5)) +Square (right) Line going through left and right square GEOMETRYCOLLECTION (LINESTRING (-0.5 0.5, 0 0.5), POLYGON ((0 0.5, 0 1, 1 1, 1 0, 0 0, 0 0.5))) +Square (right) NULL NULL +Square (right) Point middle of Left Square GEOMETRYCOLLECTION (POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0)), POINT (-0.5 0.5)) +Square (right) Point middle of Right Square POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0)) +Square (right) Square (left) POLYGON ((0 1, 1 1, 1 0, 0 0, -1 0, -1 1, 0 1)) +Square (right) Square (right) POLYGON ((1 0, 0 0, 0 1, 1 1, 1 0)) +Square (right) Square overlapping left and right square POLYGON ((1 0, 0 0, -0.1 0, -0.1 1, 0 1, 1 1, 1 0)) +Square overlapping left and right square Faraway point GEOMETRYCOLLECTION (POLYGON ((-0.1 0, 1 0, 1 1, -0.1 1, -0.1 0)), POINT (5 5)) +Square overlapping left and right square Line going through left and right square GEOMETRYCOLLECTION (LINESTRING (-0.5 0.5, -0.1 0.5), POLYGON ((-0.1 0.5, -0.1 1, 1 1, 1 0, -0.1 0, -0.1 0.5))) +Square overlapping left and right square NULL NULL +Square overlapping left and right square Point middle of Left Square GEOMETRYCOLLECTION (POLYGON ((-0.1 0, 1 0, 1 1, -0.1 1, -0.1 0)), POINT (-0.5 0.5)) +Square overlapping left and right square Point middle of Right Square POLYGON ((-0.1 0, -0.1 1, 1 1, 1 0, -0.1 0)) +Square overlapping left and right square Square (left) POLYGON ((0 0, -0.1 0, -1 0, -1 1, -0.1 1, 0 1, 1 1, 1 0, 0 0)) +Square overlapping left and right square Square (right) POLYGON ((0 0, -0.1 0, -0.1 1, 0 1, 1 1, 1 0, 0 0)) +Square overlapping left and right square Square overlapping left and right square POLYGON ((1 0, -0.1 0, -0.1 1, 1 1, 1 0)) + +query TTT SELECT a.dsc, - ST_AsEWKT(ST_Centroid(a.geom)) + b.dsc, + ST_AsEWKT(ST_Intersection(a.geom, b.geom)) +FROM geom_operators_test a +JOIN geom_operators_test b ON (1=1) +WHERE a.dsc = 'Square (left)' and b.dsc = 'Square (right)' -- Remove after #49209 is resolved. +ORDER BY a.dsc, b.dsc +---- +Square (left) Square (right) LINESTRING (0 0, 0 1) + +query TTT +SELECT + a.dsc, + ST_AsEWKT(ST_Centroid(a.geom)), + ST_AsEWKT(ST_PointOnSurface(a.geom)) FROM geom_operators_test a ORDER BY a.dsc ---- -Empty GeometryCollection POINT EMPTY -Empty LineString POINT EMPTY -Faraway point POINT (5 5) -Line going through left and right square POINT (0 0.5) -NULL NULL -Point middle of Left Square POINT (-0.5 0.5) -Point middle of Right Square POINT (0.5 0.5) -Square (left) POINT (-0.5 0.5) -Square (right) POINT (0.5 0.5) -Square overlapping left and right square POINT (0.4499999999999999 0.5) +Empty GeometryCollection POINT EMPTY POINT EMPTY +Empty LineString POINT EMPTY POINT EMPTY +Faraway point POINT (5 5) POINT (5 5) +Line going through left and right square POINT (0 0.5) POINT (-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) POINT (-0.5 0.5) POINT (-0.5 0.5) +Square (right) POINT (0.5 0.5) POINT (0.5 0.5) +Square overlapping left and right square POINT (0.4499999999999999 0.5) POINT (0.45 0.5) # Functions which take in strings as input as well. query TT diff --git a/pkg/sql/sem/builtins/geo_builtins.go b/pkg/sql/sem/builtins/geo_builtins.go index eed347031d6e..e98a019b4663 100644 --- a/pkg/sql/sem/builtins/geo_builtins.go +++ b/pkg/sql/sem/builtins/geo_builtins.go @@ -1461,43 +1461,6 @@ Note ST_Perimeter is only valid for Polygon - use ST_Length for LineString.`, tree.VolatilityImmutable, ), ), - "st_centroid": makeBuiltin( - defProps(), - geometryOverload1( - func(ctx *tree.EvalContext, g *tree.DGeometry) (tree.Datum, error) { - centroid, err := geomfn.Centroid(g.Geometry) - if err != nil { - return nil, err - } - return tree.NewDGeometry(centroid), err - }, - types.Geometry, - infoBuilder{ - info: "Returns the centroid of the given geometry.", - libraryUsage: usesGEOS, - }, - tree.VolatilityImmutable, - ), - stringOverload1( - func(ctx *tree.EvalContext, s string) (tree.Datum, error) { - g, err := geo.ParseGeometry(s) - if err != nil { - return nil, err - } - centroid, err := geomfn.Centroid(g) - if err != nil { - return nil, err - } - return tree.NewDGeometry(centroid), err - }, - types.Geometry, - infoBuilder{ - info: "Returns the centroid of the given string, which will be parsed as a geometry object.", - libraryUsage: usesGEOS, - }.String(), - tree.VolatilityImmutable, - ), - ), "st_geometrytype": makeBuiltin( defProps(), geometryOverload1( @@ -1944,6 +1907,99 @@ Note If the result has zero or one points, it will be returned as a POINT. If it }, ), + // Topology operations + "st_centroid": makeBuiltin( + defProps(), + geometryOverload1( + func(ctx *tree.EvalContext, g *tree.DGeometry) (tree.Datum, error) { + centroid, err := geomfn.Centroid(g.Geometry) + if err != nil { + return nil, err + } + return tree.NewDGeometry(centroid), err + }, + types.Geometry, + infoBuilder{ + info: "Returns the centroid of the given geometry.", + libraryUsage: usesGEOS, + }, + tree.VolatilityImmutable, + ), + stringOverload1( + func(ctx *tree.EvalContext, s string) (tree.Datum, error) { + g, err := geo.ParseGeometry(s) + if err != nil { + return nil, err + } + centroid, err := geomfn.Centroid(g) + if err != nil { + return nil, err + } + return tree.NewDGeometry(centroid), err + }, + types.Geometry, + infoBuilder{ + info: "Returns the centroid of the given string, which will be parsed as a geometry object.", + libraryUsage: usesGEOS, + }.String(), + tree.VolatilityImmutable, + ), + ), + "st_pointonsurface": makeBuiltin( + defProps(), + geometryOverload1( + func(ctx *tree.EvalContext, g *tree.DGeometry) (tree.Datum, error) { + pointOnSurface, err := geomfn.PointOnSurface(g.Geometry) + if err != nil { + return nil, err + } + return tree.NewDGeometry(pointOnSurface), err + }, + types.Geometry, + infoBuilder{ + info: "Returns a point that intersects with the given Geometry.", + libraryUsage: usesGEOS, + }, + tree.VolatilityImmutable, + ), + ), + "st_intersection": makeBuiltin( + defProps(), + geometryOverload2( + func(ctx *tree.EvalContext, a *tree.DGeometry, b *tree.DGeometry) (tree.Datum, error) { + intersection, err := geomfn.Intersection(a.Geometry, b.Geometry) + if err != nil { + return nil, err + } + return tree.NewDGeometry(intersection), err + }, + types.Geometry, + infoBuilder{ + info: "Returns the point intersections of the given geometries.", + libraryUsage: usesGEOS, + }, + tree.VolatilityImmutable, + ), + ), + "st_union": makeBuiltin( + defProps(), + geometryOverload2( + func(ctx *tree.EvalContext, a *tree.DGeometry, b *tree.DGeometry) (tree.Datum, error) { + union, err := geomfn.Union(a.Geometry, b.Geometry) + if err != nil { + return nil, err + } + return tree.NewDGeometry(union), err + }, + types.Geometry, + infoBuilder{ + info: "Returns the union of the given geometries as a single Geometry object.", + libraryUsage: usesGEOS, + }, + tree.VolatilityImmutable, + ), + ), + // // Transformations //