Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

builtins: implement ST_Intersection and ST_Buffer for geography types #51537

Merged
merged 1 commit into from
Jul 28, 2020
Merged

builtins: implement ST_Intersection and ST_Buffer for geography types #51537

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 44 additions & 0 deletions docs/generated/sql/functions.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -930,6 +930,39 @@ has no relationship with the commit order of concurrent transactions.</p>
<tr><td><a name="st_azimuth"></a><code>st_azimuth(geometry_a: geometry, geometry_b: geometry) &rarr; <a href="float.html">float</a></code></td><td><span class="funcdesc"><p>Returns the azimuth in radians of the segment defined by the given point geometries, or NULL if the two points are coincident.</p>
<p>The azimuth is angle is referenced from north, and is positive clockwise: North = 0; East = π/2; South = π; West = 3π/2.</p>
</span></td></tr>
<tr><td><a name="st_buffer"></a><code>st_buffer(geography: geography, distance: <a href="float.html">float</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>Returns a Geometry that represents all points whose distance is less than or equal to the given distance
from the given Geometry.</p>
<p>This function utilizes the GEOS module.</p>
<p>This operation is done by transforming the object into a Geometry. This occurs by translating
the Geography objects into Geometry objects before applying an LAEA, UTM or Web Mercator
based projection based on the bounding boxes of the given Geography objects. When the result is
calculated, the result is transformed back into a Geography with SRID 4326.</p>
</span></td></tr>
<tr><td><a name="st_buffer"></a><code>st_buffer(geography: geography, distance: <a href="float.html">float</a>, buffer_style_params: <a href="string.html">string</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>Returns a Geometry that represents all points whose distance is less than or equal to the given distance from the
given Geometry.</p>
<p>This variant takes in a space separate parameter string, which will augment the buffer styles. Valid parameters are:</p>
<ul>
<li>quad_segs=&lt;int&gt;, default 8</li>
<li>endcap=&lt;round|flat|butt|square&gt;, default round</li>
<li>join=&lt;round|mitre|miter|bevel&gt;, default round</li>
<li>side=&lt;both|left|right&gt;, default both</li>
<li>mitre_limit=&lt;float&gt;, default 5.0</li>
</ul>
<p>This function utilizes the GEOS module.</p>
<p>This operation is done by transforming the object into a Geometry. This occurs by translating
the Geography objects into Geometry objects before applying an LAEA, UTM or Web Mercator
based projection based on the bounding boxes of the given Geography objects. When the result is
calculated, the result is transformed back into a Geography with SRID 4326.</p>
</span></td></tr>
<tr><td><a name="st_buffer"></a><code>st_buffer(geography: geography, distance: <a href="float.html">float</a>, quad_segs: <a href="int.html">int</a>) &rarr; geography</code></td><td><span class="funcdesc"><p>Returns a Geometry that represents all points whose distance is less than or equal to the given distance from the
given Geometry.</p>
<p>This variant approximates the circle into quad_seg segments per line (the default is 8).</p>
<p>This function utilizes the GEOS module.</p>
<p>This operation is done by transforming the object into a Geometry. This occurs by translating
the Geography objects into Geometry objects before applying an LAEA, UTM or Web Mercator
based projection based on the bounding boxes of the given Geography objects. When the result is
calculated, the result is transformed back into a Geography with SRID 4326.</p>
</span></td></tr>
<tr><td><a name="st_buffer"></a><code>st_buffer(geometry: geometry, distance: <a href="decimal.html">decimal</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns a Geometry that represents all points whose distance is less than or equal to the given distance
from the given Geometry.</p>
<p>This function utilizes the GEOS module.</p>
Expand Down Expand Up @@ -1171,9 +1204,20 @@ Bottom Left.</p>
</span></td></tr>
<tr><td><a name="st_interiorringn"></a><code>st_interiorringn(geometry: geometry, n: <a href="int.html">int</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns 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.</p>
</span></td></tr>
<tr><td><a name="st_intersection"></a><code>st_intersection(geography_a: geography, geography_b: geography) &rarr; geography</code></td><td><span class="funcdesc"><p>Returns the point intersections of the given geographies.</p>
<p>This operation is done by transforming the object into a Geometry. This occurs by translating
the Geography objects into Geometry objects before applying an LAEA, UTM or Web Mercator
based projection based on the bounding boxes of the given Geography objects. When the result is
calculated, the result is transformed back into a Geography with SRID 4326.</p>
<p>This function utilizes the GEOS module.</p>
</span></td></tr>
<tr><td><a name="st_intersection"></a><code>st_intersection(geometry_a: geometry, geometry_b: geometry) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the point intersections of the given geometries.</p>
<p>This function utilizes the GEOS module.</p>
</span></td></tr>
<tr><td><a name="st_intersection"></a><code>st_intersection(geometry_a_str: <a href="string.html">string</a>, geometry_b_str: <a href="string.html">string</a>) &rarr; geometry</code></td><td><span class="funcdesc"><p>Returns the point intersections of the given geometries.</p>
<p>This function utilizes the GEOS module.</p>
<p>This variant will cast all geometry_str arguments into Geometry types.</p>
</span></td></tr>
<tr><td><a name="st_intersects"></a><code>st_intersects(geography_a: geography, geography_b: geography) &rarr; <a href="bool.html">bool</a></code></td><td><span class="funcdesc"><p>Returns true if geography_a shares any portion of space with geography_b.</p>
<p>The calculations performed are have a precision of 1cm.</p>
<p>This function utilizes the S2 library for spherical calculations.</p>
Expand Down
139 changes: 139 additions & 0 deletions pkg/geo/geogfn/best_projection.go
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
// 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 geogfn

import (
"fmt"
"math"

"github.com/cockroachdb/cockroach/pkg/geo/geopb"
"github.com/cockroachdb/cockroach/pkg/geo/geoprojbase"
"github.com/cockroachdb/errors"
"github.com/golang/geo/s1"
"github.com/golang/geo/s2"
)

// BestGeomProjection translates roughly to the ST_BestSRID function in PostGIS.
// It attempts to find the best projection for a bounding box into an accurate
// geometry-type projection.
//
// The algorithm is described by ST_Buffer/ST_Intersection documentation (paraphrased):
// It first determines the best SRID that fits the bounding box of the 2 geography objects (ST_Intersection only).
// It favors a north/south pole projection, then UTM, then LAEA for smaller zones, otherwise falling back
// to web mercator.
// If geography objects are within one half zone UTM but not the same UTM it will pick one of those.
// After the calculation is complete, it will fall back to WGS84 Geography.
func BestGeomProjection(boundingRect s2.Rect) (geoprojbase.Proj4Text, error) {
center := boundingRect.Center()

latWidth := s1.Angle(boundingRect.Lat.Length())
lngWidth := s1.Angle(boundingRect.Lng.Length())

// Check if these fit either the North Pole or South Pole areas.
// If the center has latitude greater than 70 (an arbitrary polar threshold), and it is
// within the polar ranges, return that.
if center.Lat.Degrees() > 70 && boundingRect.Lo().Lat.Degrees() > 45 {
// See: https://epsg.io/3574.
return getGeomProjection(3574)
}
// Same for south pole.
if center.Lat.Degrees() < -70 && boundingRect.Hi().Lat.Degrees() < -45 {
// See: https://epsg.io/3409
return getGeomProjection(3409)
}

// Each UTM zone is 6 degrees wide and distortion is low for geometries that fit within the zone. We use
// UTM if the width is lower than the UTM zone width, even though the geometry may span 2 zones -- using
// the geometry center to pick the UTM zone should result in most of the geometry being in the picked zone.
if lngWidth.Degrees() < 6 {
// Determine the offset of the projection.
// Offset longitude -180 to 0 and divide by 6 to get the zone.
// Note that we treat 180 degree longtitudes as offset 59.
// TODO(#geo): do we care about https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions?
// PostGIS's _ST_BestSRID function doesn't seem to care: .
sridOffset := geopb.SRID(math.Min(math.Floor((center.Lng.Degrees()+180)/6), 59))
if center.Lat.Degrees() >= 0 {
// Start at the north UTM SRID.
return getGeomProjection(32601 + sridOffset)
}
// Start at the south UTM SRID.
// This should make no difference in end result compared to using the north UTMs,
// but for completeness we do it.
return getGeomProjection(32701 + sridOffset)
}

// Attempt to fit into LAEA areas if the width is less than 25 degrees (we can go up to 30
// but want to leave some room for precision issues).
//
// LAEA areas are separated into 3 latitude zones between 0 and 90 and 3 latitude zones
// between -90 and 0. Within each latitude zones, they have different longitude bands:
// * The bands closest to the equator have 12x30 degree longitude zones.
// * The bands in the temperate area 8x45 degree longitude zones.
// * The bands near the poles have 4x90 degree longitude zones.
//
// For each of these bands, we custom define a LAEA area with the center of the LAEA area
// as the lat/lon offset.
//
// See also: https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection.
if latWidth.Degrees() < 25 {
// Convert lat to a known 30 degree zone..
// -3 represents [-90, -60), -2 represents [-60, -30) ... and 2 represents [60, 90].
// (note: 90 is inclusive at the end).
// Treat a 90 degree latitude as band 2.
latZone := math.Min(math.Floor(center.Lat.Degrees()/30), 2)
latZoneCenterDegrees := (latZone * 30) + 15
// Equator bands - 30 degree zones.
if (latZone == 0 || latZone == -1) && lngWidth.Degrees() <= 30 {
lngZone := math.Floor(center.Lng.Degrees() / 30)
return geoprojbase.MakeProj4Text(
fmt.Sprintf(
"+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=%g +lon_0=%g +units=m +no_defs",
latZoneCenterDegrees,
(lngZone*30)+15,
),
), nil
}
// Temperate bands - 45 degree zones.
if (latZone == -2 || latZone == 1) && lngWidth.Degrees() <= 45 {
lngZone := math.Floor(center.Lng.Degrees() / 45)
return geoprojbase.MakeProj4Text(
fmt.Sprintf(
"+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=%g +lon_0=%g +units=m +no_defs",
latZoneCenterDegrees,
(lngZone*45)+22.5,
),
), nil
}
// Polar bands -- 90 degree zones.
if (latZone == -3 || latZone == 2) && lngWidth.Degrees() <= 90 {
lngZone := math.Floor(center.Lng.Degrees() / 90)
return geoprojbase.MakeProj4Text(
fmt.Sprintf(
"+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=%g +lon_0=%g +units=m +no_defs",
latZoneCenterDegrees,
(lngZone*90)+45,
),
), nil
}
}

// Default to Web Mercator.
return getGeomProjection(3857)
}

// getGeomProjection returns the Proj4Text associated with an SRID.
func getGeomProjection(srid geopb.SRID) (geoprojbase.Proj4Text, error) {
proj, ok := geoprojbase.Projection(srid)
if !ok {
return geoprojbase.Proj4Text{}, errors.Newf("unexpected SRID %d", srid)
}
return proj.Proj4Text, nil
}
94 changes: 94 additions & 0 deletions pkg/geo/geogfn/best_projection_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
// 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 geogfn

import (
"testing"

"github.com/cockroachdb/cockroach/pkg/geo/geoprojbase"
"github.com/golang/geo/s2"
"github.com/stretchr/testify/require"
)

func TestBestGeomProjection(t *testing.T) {
testCases := []struct {
desc string
rect s2.Rect
expected geoprojbase.Proj4Text
}{
{
"north pole",
s2.RectFromLatLng(s2.LatLngFromDegrees(75, 75)),
geoprojbase.Projections[3574].Proj4Text,
},
{
"south pole",
s2.RectFromLatLng(s2.LatLngFromDegrees(-75, -75)),
geoprojbase.Projections[3409].Proj4Text},
{
"utm 15 on top hemisphere",
s2.RectFromLatLng(s2.LatLngFromDegrees(15, 93)),
geoprojbase.Projections[32646].Proj4Text},
{
"utm -16 on bottom hemisphere",
s2.RectFromLatLng(s2.LatLngFromDegrees(-15, -111)),
geoprojbase.Projections[32712].Proj4Text,
},
{
"LAEA at equator bands (north half)",
s2.RectFromCenterSize(s2.LatLngFromDegrees(12, 22), s2.LatLngFromDegrees(10, 10)),
geoprojbase.MakeProj4Text("+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=15 +lon_0=15 +units=m +no_defs"),
},
{
"LAEA at equator bands (south half)",
s2.RectFromCenterSize(s2.LatLngFromDegrees(-13, 123), s2.LatLngFromDegrees(10, 10)),
geoprojbase.MakeProj4Text("+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=-15 +lon_0=135 +units=m +no_defs"),
},
{
"LAEA at temperate band (north half)",
s2.RectFromCenterSize(s2.LatLngFromDegrees(33, 87), s2.LatLngFromDegrees(10, 10)),
geoprojbase.MakeProj4Text("+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=45 +lon_0=67.5 +units=m +no_defs"),
},
{
"LAEA at temperate band (south half)",
s2.RectFromCenterSize(s2.LatLngFromDegrees(-53, -120.5), s2.LatLngFromDegrees(10, 10)),
geoprojbase.MakeProj4Text("+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=-45 +lon_0=-112.5 +units=m +no_defs"),
},
{
"LAEA at polar band (north half)",
s2.RectFromCenterSize(s2.LatLngFromDegrees(63, 87), s2.LatLngFromDegrees(10, 10)),
geoprojbase.MakeProj4Text("+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=75 +lon_0=45 +units=m +no_defs"),
},
{
"LAEA at polar band (south half)",
s2.RectFromCenterSize(s2.LatLngFromDegrees(-66, -120.5), s2.LatLngFromDegrees(10, 10)),
geoprojbase.MakeProj4Text("+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=-75 +lon_0=-135 +units=m +no_defs"),
},
{
"UTM which should be 32V, but we return 31V as we do not handle UTM exceptions",
s2.RectFromLatLng(s2.LatLngFromDegrees(59.4136, 5.26)),
geoprojbase.Projections[32631].Proj4Text, // Should be 32632
},
{
"wide example",
s2.RectFromCenterSize(s2.LatLngFromDegrees(0, 0), s2.LatLngFromDegrees(50, 50)),
geoprojbase.Projections[3857].Proj4Text,
},
}

for _, tc := range testCases {
t.Run(tc.desc, func(t *testing.T) {
result, err := BestGeomProjection(tc.rect)
require.NoError(t, err)
require.Equal(t, tc.expected.String(), result.String())
})
}
}
Loading