Skip to content

Commit

Permalink
Merge #49618 #49649
Browse files Browse the repository at this point in the history 
49618: geo: audit functions for duplicate points in Geometry r=sumeerbhola a=otan

Audited all the existing functions to check they work for duplicate
point cases (e.g. `LINESTRING(1.0 1.0, 1.0 1.0, 2.0 2.0)`. The ones I
found that were handled incorrect were distance. In particular, distance
is completely wrong for s2's Distance calculator, so I've opted to
replace that with our own in house instead.

Added a few small tests. If more are incorrect, it will be up to
randomized testing to handle it.

Release note: None

49649: builtins: implement more PostGIS tutorial available builtins r=sumeerbhola a=otan

Release note (sql change): Implemented the following builtins available
on Geometry types:
* ST_Area2D (resolves #48869)
* ST_NRings (resolves #48996)
* ST_PointN (resolves #49008)
* ST_GeometryType (resolves #48946)

Co-authored-by: Oliver Tan <otan@cockroachlabs.com>
  • Loading branch information
@otan
craig[bot] and otan committed May 28, 2020
3 parents ca57e95 + 195eff7 + f4aa72e commit 3f08089
Show file tree
Hide file tree
Showing 9 changed files with 253 additions and 275 deletions.
10 changes: 10 additions & 0 deletions docs/generated/sql/functions.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -701,6 +701,9 @@ has no relationship with the commit order of concurrent transactions.</p>
<tr><td><a name="st_area"></a><code>st_area(geometry: geometry) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the area of the given geometry.</p>
<p>This function utilizes the GEOS module.</p>
</span></td></tr>
<tr><td><a name="st_area2d"></a><code>st_area2d(geometry: geometry) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the area of the given geometry.</p>
<p>This function utilizes the GEOS module.</p>
</span></td></tr>
<tr><td><a name="st_asbinary"></a><code>st_asbinary(geography: geography) &rarr; <a href="bytes.html">bytes</a></code></td><td><span class="funcdesc"><p>Returns the WKB representation of a given Geography.</p>
</span></td></tr>
<tr><td><a name="st_asbinary"></a><code>st_asbinary(geography: geography, xdr_or_ndr: <a href="string.html">string</a>) &rarr; <a href="bytes.html">bytes</a></code></td><td><span class="funcdesc"><p>Returns the WKB representation of a given Geography. This variant has a second argument denoting the encoding - <code>xdr</code> for big endian and <code>ndr</code> for little endian.</p>
Expand Down Expand Up @@ -844,6 +847,9 @@ has no relationship with the commit order of concurrent transactions.</p>
</span></td></tr>
<tr><td><a name="st_geometryn"></a><code>st_geometryn(geometry: geometry, n: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the n-th Geometry (1-indexed). Returns NULL if out of bounds.</p>
</span></td></tr>
<tr><td><a name="st_geometrytype"></a><code>st_geometrytype(geometry: geometry) &rarr; <a href="string.html">string</a></code></td><td><span class="funcdesc"><p>Returns the type of geometry as a string.</p>
<p>This function utilizes the GEOS module.</p>
</span></td></tr>
<tr><td><a name="st_geomfromewkb"></a><code>st_geomfromewkb(val: <a href="bytes.html">bytes</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from an EWKB representation.</p>
</span></td></tr>
<tr><td><a name="st_geomfromewkt"></a><code>st_geomfromewkt(val: <a href="string.html">string</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from an EWKT representation.</p>
Expand Down Expand Up @@ -971,6 +977,8 @@ has no relationship with the commit order of concurrent transactions.</p>
</span></td></tr>
<tr><td><a name="st_npoints"></a><code>st_npoints(geometry: geometry) &rarr; <a href="int.html">int</a></code></td><td><span class="funcdesc"><p>Returns the number of points in a given Geometry. Works for any shape type.</p>
</span></td></tr>
<tr><td><a name="st_nrings"></a><code>st_nrings(geometry: geometry) &rarr; <a href="int.html">int</a></code></td><td><span class="funcdesc"><p>Returns the number of rings in a Polygon Geometry. Returns 0 if the shape is not a Polygon.</p>
</span></td></tr>
<tr><td><a name="st_numgeometries"></a><code>st_numgeometries(geometry: geometry) &rarr; <a href="int.html">int</a></code></td><td><span class="funcdesc"><p>Returns the number of shapes inside a given Geometry.</p>
</span></td></tr>
<tr><td><a name="st_numinteriorring"></a><code>st_numinteriorring(geometry: geometry) &rarr; <a href="int.html">int</a></code></td><td><span class="funcdesc"><p>Returns the number of interior rings in a Polygon Geometry. Returns NULL if the shape is not a Polygon.</p>
Expand Down Expand Up @@ -1004,6 +1012,8 @@ has no relationship with the commit order of concurrent transactions.</p>
</span></td></tr>
<tr><td><a name="st_pointfromwkb"></a><code>st_pointfromwkb(wkb: <a href="bytes.html">bytes</a>, srid: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from a WKB representation with an SRID. If the shape underneath is not Point, NULL is returned.</p>
</span></td></tr>
<tr><td><a name="st_pointn"></a><code>st_pointn(geometry: geometry, n: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the n-th Point of a LineString (1-indexed). Returns NULL if out of bounds or not a LineString.</p>
</span></td></tr>
<tr><td><a name="st_polyfromtext"></a><code>st_polyfromtext(str: <a href="string.html">string</a>, srid: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns 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.</p>
</span></td></tr>
<tr><td><a name="st_polyfromtext"></a><code>st_polyfromtext(val: <a href="string.html">string</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the Geometry from a WKT or EWKT representation. If the shape underneath is not Polygon, NULL is returned.</p>
Expand Down
176 changes: 49 additions & 127 deletions pkg/geo/geogfn/distance.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,7 @@ func Distance(
return 0, err
}
spheroid := geographiclib.WGS84Spheroid
if useSphereOrSpheroid == UseSpheroid {
return distanceSpheroidRegions(spheroid, aRegions, bRegions, 0.0)
}
return distanceSphereRegions(spheroid, aRegions, bRegions)
return distanceGeographyRegions(spheroid, useSphereOrSpheroid, aRegions, bRegions, 0)
}

//
Expand Down Expand Up @@ -163,11 +160,11 @@ func (c *s2GeodistEdgeCrosser) ChainCrossing(p geodist.Point) bool {
return c.EdgeCrosser.ChainCrossingSign(p.(*s2GeodistPoint).Point) != s2.DoNotCross
}

// distanceSpheroidRegions calculates the distance between two sets of regions on a spheroid.
// distanceGeographyRegions calculates the distance between two sets of regions.
// It will quit if it finds a distance that is less than stopAfterLE.
// It is not guaranteed to find the absolute minimum distance if stopAfterLE > 0.
//
// !!! SURPRISING BEHAVIOR WARNING !!!
// !!! SURPRISING BEHAVIOR WARNING FOR SPHEROIDS !!!
// PostGIS evaluates the distance between spheroid regions by computing the min of
// the pair-wise distance between the cross-product of the regions in A and the regions
// in B, where the pair-wise distance is computed as:
Expand All @@ -181,8 +178,12 @@ func (c *s2GeodistEdgeCrosser) ChainCrossing(p geodist.Point) bool {
// Since we aim to be compatible with PostGIS, we adopt the same approach.
//
// TODO(otan): accelerate checks with bounding boxes.
func distanceSpheroidRegions(
spheroid geographiclib.Spheroid, aRegions []s2.Region, bRegions []s2.Region, stopAfterLE float64,
func distanceGeographyRegions(
spheroid geographiclib.Spheroid,
useSphereOrSpheroid UseSphereOrSpheroid,
aRegions []s2.Region,
bRegions []s2.Region,
stopAfterLE float64,
) (float64, error) {
minDistance := math.MaxFloat64
for _, aRegion := range aRegions {
Expand All @@ -191,13 +192,13 @@ func distanceSpheroidRegions(
return 0, err
}
for _, bRegion := range bRegions {
minDistanceUpdater := newSpheroidMinDistanceUpdater(spheroid, stopAfterLE)
minDistanceUpdater := newGeographyMinDistanceUpdater(spheroid, useSphereOrSpheroid, stopAfterLE)
bGeodist, err := regionToGeodistShape(bRegion)
if err != nil {
return 0, err
}
earlyExit, err := geodist.ShapeDistance(
&spheroidDistanceCalculator{updater: minDistanceUpdater},
&geographyDistanceCalculator{updater: minDistanceUpdater},
aGeodist,
bGeodist,
)
Expand All @@ -213,45 +214,55 @@ func distanceSpheroidRegions(
return minDistance, nil
}

// spheroidMinDistanceUpdater finds the minimum distance using spheroid calculations.
// geographyMinDistanceUpdater finds the minimum distance using a sphere.
// Methods will return early if it finds a minimum distance <= stopAfterLE.
type spheroidMinDistanceUpdater struct {
spheroid geographiclib.Spheroid
minEdge s2.Edge
minD s1.ChordAngle
stopAfterLE s1.ChordAngle
type geographyMinDistanceUpdater struct {
spheroid geographiclib.Spheroid
useSphereOrSpheroid UseSphereOrSpheroid
minEdge s2.Edge
minD s1.ChordAngle
stopAfterLE s1.ChordAngle
}

var _ geodist.DistanceUpdater = (*spheroidMinDistanceUpdater)(nil)
var _ geodist.DistanceUpdater = (*geographyMinDistanceUpdater)(nil)

// newSpheroidMinDistanceUpdater returns a new spheroidMinDistanceUpdater with the
// newGeographyMinDistanceUpdater returns a new geographyMinDistanceUpdater with the
// correct arguments set up.
func newSpheroidMinDistanceUpdater(
spheroid geographiclib.Spheroid, stopAfterLE float64,
) *spheroidMinDistanceUpdater {
return &spheroidMinDistanceUpdater{
spheroid: spheroid,
minD: math.MaxFloat64,
func newGeographyMinDistanceUpdater(
spheroid geographiclib.Spheroid, useSphereOrSpheroid UseSphereOrSpheroid, stopAfterLE float64,
) *geographyMinDistanceUpdater {
multiplier := 1.0
if useSphereOrSpheroid == UseSpheroid {
// Modify the stopAfterLE distance to be 95% of the proposed stopAfterLE, since
// we use the sphere to calculate the distance and we want to leave a buffer
// for spheroid distances being slightly off.
// Distances should differ by a maximum of (2 * spheroid.Flattening)%,
// so the 5% margin is pretty safe.
stopAfterLE: s1.ChordAngleFromAngle(s1.Angle(stopAfterLE/spheroid.SphereRadius)) * 0.95,
multiplier = 0.95
}
stopAfterLEChordAngle := s1.ChordAngleFromAngle(s1.Angle(stopAfterLE * multiplier / spheroid.SphereRadius))
return &geographyMinDistanceUpdater{
spheroid: spheroid,
minD: math.MaxFloat64,
useSphereOrSpheroid: useSphereOrSpheroid,
stopAfterLE: stopAfterLEChordAngle,
}
}

// Distance implements the DistanceUpdater interface.
func (u *spheroidMinDistanceUpdater) Distance() float64 {
func (u *geographyMinDistanceUpdater) Distance() float64 {
// If the distance is zero, avoid the call to spheroidDistance and return early.
if u.minD == 0 {
return 0
}
return spheroidDistance(u.spheroid, u.minEdge.V0, u.minEdge.V1)
if u.useSphereOrSpheroid == UseSpheroid {
return spheroidDistance(u.spheroid, u.minEdge.V0, u.minEdge.V1)
}
return u.minD.Angle().Radians() * u.spheroid.SphereRadius
}

// Update implements the geodist.DistanceUpdater interface.
func (u *spheroidMinDistanceUpdater) Update(
func (u *geographyMinDistanceUpdater) Update(
aInterface geodist.Point, bInterface geodist.Point,
) bool {
a := aInterface.(*s2GeodistPoint).Point
Expand All @@ -272,30 +283,30 @@ func (u *spheroidMinDistanceUpdater) Update(
}

// OnIntersects implements the geodist.DistanceUpdater interface.
func (u *spheroidMinDistanceUpdater) OnIntersects() bool {
func (u *geographyMinDistanceUpdater) OnIntersects() bool {
u.minD = 0
return true
}

// IsMaxDistance implements the geodist.DistanceUpdater interface.
func (u *spheroidMinDistanceUpdater) IsMaxDistance() bool {
func (u *geographyMinDistanceUpdater) IsMaxDistance() bool {
return false
}

// spheroidDistanceCalculator implements geodist.DistanceCalculator
type spheroidDistanceCalculator struct {
updater *spheroidMinDistanceUpdater
// geographyDistanceCalculator implements geodist.DistanceCalculator
type geographyDistanceCalculator struct {
updater *geographyMinDistanceUpdater
}

var _ geodist.DistanceCalculator = (*spheroidDistanceCalculator)(nil)
var _ geodist.DistanceCalculator = (*geographyDistanceCalculator)(nil)

// DistanceUpdater implements geodist.DistanceCalculator.
func (c *spheroidDistanceCalculator) DistanceUpdater() geodist.DistanceUpdater {
func (c *geographyDistanceCalculator) DistanceUpdater() geodist.DistanceUpdater {
return c.updater
}

// NewEdgeCrosser implements geodist.DistanceCalculator.
func (c *spheroidDistanceCalculator) NewEdgeCrosser(
func (c *geographyDistanceCalculator) NewEdgeCrosser(
edge geodist.Edge, startPoint geodist.Point,
) geodist.EdgeCrosser {
return &s2GeodistEdgeCrosser{
Expand All @@ -308,7 +319,7 @@ func (c *spheroidDistanceCalculator) NewEdgeCrosser(
}

// PointInLinearRing implements geodist.DistanceCalculator.
func (c *spheroidDistanceCalculator) PointInLinearRing(
func (c *geographyDistanceCalculator) PointInLinearRing(
point geodist.Point, polygon geodist.LinearRing,
) bool {
return polygon.(*s2GeodistLinearRing).ContainsPoint(point.(*s2GeodistPoint).Point)
Expand All @@ -323,7 +334,7 @@ func (c *spheroidDistanceCalculator) PointInLinearRing(
//
// For visualization and more, see: Section 6 / Figure 4 of
// "Projective configuration theorems: old wine into new wineskins", Tabachnikov, Serge, 2016/07/16
func (c *spheroidDistanceCalculator) ClosestPointToEdge(
func (c *geographyDistanceCalculator) ClosestPointToEdge(
edge geodist.Edge, pointInterface geodist.Point,
) (geodist.Point, bool) {
eV0 := edge.V0.(*s2GeodistPoint).Point
Expand Down Expand Up @@ -359,92 +370,3 @@ func regionToGeodistShape(r s2.Region) (geodist.Shape, error) {
}
return nil, errors.Newf("unknown region: %T", r)
}

//
// Spheres
//

// distanceSphereRegions calculates the distance between two objects on a sphere.
func distanceSphereRegions(
spheroid geographiclib.Spheroid, aRegions []s2.Region, bRegions []s2.Region,
) (float64, error) {
aShapeIndex, aPoints, err := s2RegionsToPointsAndShapeIndexes(aRegions)
if err != nil {
return 0, err
}
bShapeIndex, bPoints, err := s2RegionsToPointsAndShapeIndexes(bRegions)
if err != nil {
return 0, err
}

minDistanceUpdater := newSphereMinDistanceUpdater(spheroid)
// Compare aShapeIndex to bShapeIndex as well as all points in B.
if aShapeIndex.Len() > 0 {
if bShapeIndex.Len() > 0 {
if minDistanceUpdater.onShapeIndexToShapeIndex(aShapeIndex, bShapeIndex) {
return minDistanceUpdater.minDistance(), nil
}
}
for _, bPoint := range bPoints {
if minDistanceUpdater.onShapeIndexToPoint(aShapeIndex, bPoint) {
return minDistanceUpdater.minDistance(), nil
}
}
}

// Then try compare all A points against bShapeIndex and bPoints.
for _, aPoint := range aPoints {
if bShapeIndex.Len() > 0 {
if minDistanceUpdater.onShapeIndexToPoint(bShapeIndex, aPoint) {
return minDistanceUpdater.minDistance(), nil
}
}
for _, bPoint := range bPoints {
if minDistanceUpdater.onPointToPoint(aPoint, bPoint) {
return minDistanceUpdater.minDistance(), nil
}
}
}
return minDistanceUpdater.minDistance(), nil
}

// sphereMinDistanceUpdater finds the minimum distance on a sphere.
type sphereMinDistanceUpdater struct {
spheroid geographiclib.Spheroid
minD s1.ChordAngle
}

// newSphereMinDistanceUpdater returns a new sphereMinDistanceUpdater.
func newSphereMinDistanceUpdater(spheroid geographiclib.Spheroid) *sphereMinDistanceUpdater {
return &sphereMinDistanceUpdater{spheroid: spheroid, minD: s1.StraightChordAngle}
}

// onShapeIndexToShapeIndex updates the minimum distance and returns true if distance is 0.
func (u *sphereMinDistanceUpdater) onShapeIndexToShapeIndex(
a *s2.ShapeIndex, b *s2.ShapeIndex,
) bool {
u.minD = minChordAngle(u.minD, s2.NewClosestEdgeQuery(a, nil).Distance(s2.NewMinDistanceToShapeIndexTarget(b)))
return u.minD == 0
}

// onShapeIndexToPoint updates the minimum distance and returns true if distance is 0.
func (u *sphereMinDistanceUpdater) onShapeIndexToPoint(a *s2.ShapeIndex, b s2.Point) bool {
u.minD = minChordAngle(u.minD, s2.NewClosestEdgeQuery(a, nil).Distance(s2.NewMinDistanceToPointTarget(b)))
return u.minD == 0
}

// onPointToPoint updates the minimum distance and return true if the distance is 0.
func (u *sphereMinDistanceUpdater) onPointToPoint(a s2.Point, b s2.Point) bool {
if a == b {
u.minD = 0
return true
}
u.minD = minChordAngle(u.minD, s2.ChordAngleBetweenPoints(a, b))
return u.minD == 0
}

// minDistance returns the minimum distance in meters found in the sphereMinDistanceUpdater
// so far.
func (u *sphereMinDistanceUpdater) minDistance() float64 {
return u.minD.Angle().Radians() * u.spheroid.SphereRadius
}
28 changes: 28 additions & 0 deletions pkg/geo/geogfn/distance_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,20 @@ var distanceTestCases = []struct {
628519.03378753,
627129.50261075,
},
{
"LINESTRING to intersecting LINESTRING, duplicate points",
"LINESTRING(0.0 0.0, 0.0 0.0, 1.0 1.0)",
"LINESTRING(0.0 1.0, 1.0 0.0, 1.0 0.0)",
0,
0,
},
{
"LINESTRING to faraway LINESTRING, duplicate points",
"LINESTRING(0.0 0.0, 1.0 1.0)",
"LINESTRING(5.0 5.0, 6.0 6.0, 6.0 6.0, 6.0 6.0)",
628519.03378753,
627129.50261075,
},
{
"LINESTRING to intersecting POLYGON where LINESTRING is inside the polygon",
"POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 1.0, 0.0 0.0))",
Expand Down Expand Up @@ -173,6 +187,20 @@ var distanceTestCases = []struct {
0,
0,
},
{
"POLYGON to POLYGON completely inside its hole, duplicate points",
"POLYGON((0.0 0.0, 0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 1.0, 0.0 0.0), (0.2 0.2, 0.2 0.4, 0.2 0.4, 0.4 0.4, 0.4 0.2, 0.2 0.2))",
"POLYGON((0.25 0.25, 0.25 0.25, 0.35 0.35, 0.25 0.35, 0.25 0.35, 0.25 0.25))",
5559.65025416,
5565.87138621,
},
{
"POLYGON to POLYGON intersecting through its hole, duplicate points",
"POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 1.0, 0.0 1.0, 0.0 0.0), (0.2 0.2, 0.2 0.4, 0.2 0.4, 0.4 0.4, 0.4 0.2, 0.2 0.2))",
"POLYGON((0.15 0.25, 0.15 0.25, 0.35 0.25, 0.35 0.35, 0.25 0.35, 0.25 0.35, 0.15 0.25))",
0,
0,
},
{
"POLYGON to faraway POLYGON",
"POLYGON((0.0 0.0, 1.0 0.0, 1.0 1.0, 0.0 1.0, 0.0 0.0))",
Expand Down
Loading

0 comments on commit 3f08089

Please sign in to comment.