github-actions[bot] commented on code in PR #69032:
URL: https://github.com/apache/doris/pull/69032#discussion_r4378237155
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -382,16 +1613,25 @@ struct StAreaSquareMeters {
auto return_type = block.get_data_type(result);
auto col =
block.get_by_position(arguments[0]).column->convert_to_full_column_if_const();
+ const auto& input_type = block.get_data_type(arguments[0]);
+ RETURN_IF_ERROR(validate_geography_semantics(input_type, NAME));
const auto size = col->size();
auto res = ColumnFloat64::create();
res->reserve(size);
auto null_map = ColumnUInt8::create(size, 0);
auto& null_map_data = null_map->get_data();
std::unique_ptr<GeoShape> shape;
- for (int row = 0; row < size; ++row) {
+ for (size_t row = 0; row < size; ++row) {
+ // Propagate NULL before using nested placeholders; a NULL
geometry has an
+ // empty payload, not a WKB shape.
+ if (any_geo_argument_is_null(block, arguments, row)) {
+ null_map_data[row] = 1;
+ res->insert_default();
+ continue;
+ }
auto shape_value = col->get_data_at(row);
- shape = GeoShape::from_encoded(shape_value.data, shape_value.size);
+ RETURN_IF_ERROR(decode_calculation_shape(shape_value, input_type,
&shape, NAME));
Review Comment:
[P2] Stabilize typed polygon rings before computing area
For default GEOGRAPHY WKB, let A be `POLYGON((-44.99999999999997
-45,134.99999999999997 45,45 45,-44.99999999999997 -45))` and B the same ring
with `(45 0)` inserted on the first edge. These are valid rings already
admitted by the new topology tests, and that inserted point is the exact
raw-arc midpoint, so their areas must agree. This decode still builds
independently rounded S2 vertices via `to_s2loop`; neither area unit uses the
stable raw-arc representation. Modeling the pinned [S2 area
arithmetic](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2measures.cc)
and loop integration gives about 4.69446e13 versus 4.99644e13 square meters, a
roughly 6% difference solely from subdivision (source/arithmetic evidence, not
a Doris runtime test). This is a separate measurement consumer from existing
polygon-interior finding 4377327193 and distance finding 4376498668. Use stable
raw-coordinate rings for both typed area functions, preserving hole roles and
legacy STRING behavior
, and cover subdivision/reversal/holes in both units.
##########
fe/fe-core/src/main/java/org/apache/doris/nereids/rules/analysis/BindSink.java:
##########
@@ -1005,6 +1006,8 @@ private Plan
bindConnectorTableSink(MatchingContext<UnboundConnectorTableSink<Pl
throw new AnalysisException("insert into cols should be
corresponding to the query output. "
+ "Expected " + boundSink.getCols().size() + " columns but
got " + child.getOutput().size());
}
+ ConnectorSpatialWriteValidator.validate(bindColumns,
child.getOutput().stream()
Review Comment:
[P2] Preserve typed NULL provenance for INSERT SELECT
For the default `geometry_values(id INT, geom GEOMETRY, geog GEOGRAPHY)`
fixture, `INSERT INTO geometry_values SELECT 5, CAST(NULL AS
GEOMETRY('EPSG:3857')), NULL` is rejected as a CRS mismatch although the
equivalent VALUES row is explicitly supported. A reduced plan inferred from
this path is:
```text
ConnectorSink[geom target=GEOMETRY('OGC:CRS84')]
OneRowRelation[CAST(NULL AS GEOMETRY('EPSG:3857')) AS geom]
```
`InsertUtils.normalizePlanWithoutLock` only normalizes `UnboundInlineTable`;
SELECT reaches `BindSink` before later constant rewrites.
`LogicalOneRowRelation.computeOutput`/`Alias.toSlot` expose the producer as a
slot retaining EPSG:3857 but losing its NULL value. Passing only those slots
here prevents `ConnectorSpatialWriteValidator`'s `isNullLiteral` exemption from
matching, so `validateSpatialType` throws. `CAST(NULL AS INT)` and nested
typed-NULL leaves have the same problem. Resolve proven NULL output lineage and
normalize it to the target before sink coercion, while retaining rejection of
non-NULL identity mismatches; cover SELECT/project/CTE alongside the VALUES
controls. Existing 4290141173 repaired actual VALUES literals; this is the
separate SELECT producer-to-slot validation boundary in BindSink.
##########
be/src/exprs/function/geo/functions_geo.cpp:
##########
@@ -18,32 +18,1163 @@
#include "exprs/function/geo/functions_geo.h"
#include <glog/logging.h>
+#include <s2/mutable_s2shape_index.h>
+#include <s2/s2boolean_operation.h>
+#include <s2/s2builderutil_s2polyline_vector_layer.h>
+#include <s2/s2closest_edge_query.h>
+#include <s2/s2closest_point_query.h>
+#include <s2/s2contains_point_query.h>
+#include <s2/s2earth.h>
+#include <s2/s2lax_polygon_shape.h>
+#include <s2/s2loop.h>
+#include <s2/s2loop_measures.h>
#include <s2/s2point.h>
+#include <s2/s2point_index.h>
+#include <s2/s2point_vector_shape.h>
+#include <s2/s2pointutil.h>
+#include <s2/s2polygon.h>
+#include <s2/s2polyline.h>
+#include <s2/s2predicates.h>
#include <algorithm>
+#include <array>
+#include <bit>
#include <boost/iterator/iterator_facade.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <numbers>
+#include <optional>
+#include <set>
#include <utility>
+#include "common/cast_set.h"
#include "common/compiler_util.h"
#include "core/assert_cast.h"
#include "core/block/block.h"
#include "core/block/column_with_type_and_name.h"
#include "core/column/column.h"
#include "core/column/column_array.h"
#include "core/column/column_execute_util.h"
+#include "core/column/column_map.h"
#include "core/column/column_nullable.h"
+#include "core/column/column_struct.h"
+#include "core/column/column_varbinary.h"
+#include "core/custom_allocator.h"
#include "core/data_type/data_type_array.h"
+#include "core/data_type/data_type_map.h"
#include "core/data_type/data_type_nullable.h"
#include "core/data_type/data_type_number.h"
+#include "core/data_type/data_type_spatial.h"
#include "core/data_type/data_type_string.h"
+#include "core/data_type/data_type_struct.h"
+#include "core/data_type/data_type_varbinary.h"
#include "core/data_type/define_primitive_type.h"
#include "core/string_ref.h"
#include "exprs/function/geo/geo_common.h"
#include "exprs/function/geo/geo_types.h"
+#include "exprs/function/geo/wkb_parse.h"
#include "exprs/function/simple_function_factory.h"
+#include "exprs/function/string_hex_util.h"
namespace doris {
+static bool is_spatial_type(const DataTypePtr& type) {
+ const auto primitive_type = remove_nullable(type)->get_primitive_type();
+ return primitive_type == TYPE_GEOMETRY || primitive_type == TYPE_GEOGRAPHY;
+}
+
+static Status validate_geography_semantics(const DataTypePtr& type, const
char* function_name) {
+ const auto& nested_type = remove_nullable(type);
+ if (!is_spatial_type(nested_type)) {
+ return Status::OK();
+ }
+
+ const auto* spatial_type = dynamic_cast<const
DataTypeSpatial*>(nested_type.get());
+ DCHECK(spatial_type != nullptr);
+ if (spatial_type->get_primitive_type() == TYPE_GEOGRAPHY &&
+ spatial_type->crs() == "OGC:CRS84" && spatial_type->algorithm() ==
"spherical") {
+ return Status::OK();
+ }
+ return Status::NotSupported(
+ "Function {} requires GEOGRAPHY(OGC:CRS84, spherical) for spatial
inputs",
+ function_name);
+}
+
+static std::unique_ptr<GeoShape> decode_geo_shape(StringRef value, const
DataTypePtr& type) {
+ if (!is_spatial_type(type)) {
+ return GeoShape::from_encoded(value.data, value.size);
+ }
+
+ GeoParseStatus status;
+ auto shape = GeoShape::from_wkb_bytes(value.data, value.size, status);
+ return status == GeoParseStatus::GEO_PARSE_OK ? std::move(shape) : nullptr;
+}
+
+static bool any_geo_argument_is_null(const Block& block, const ColumnNumbers&
arguments,
+ size_t row) {
+ return std::ranges::any_of(arguments, [&](auto argument) {
+ return block.get_by_position(argument).column->is_null_at(row);
+ });
+}
+
+static Status decode_calculation_shape(StringRef value, const DataTypePtr&
type,
+ std::unique_ptr<GeoShape>* shape,
+ const char* function_name) {
+ *shape = decode_geo_shape(value, type);
+ // POLYGON EMPTY decodes successfully, but distance helpers retain
infinity/max sentinels.
+ // Reject it for typed calculations just like empty points/lines, leaving
raw WKB untouched.
+ if (is_spatial_type(type) &&
+ (*shape == nullptr ||
+ ((*shape)->type() == GeoShapeType::GEO_SHAPE_POLYGON &&
+ assert_cast<const GeoPolygon&>(**shape).polygon()->num_loops() ==
0))) {
+ return Status::NotSupported("Function {} does not support this
Geography WKB shape",
+ function_name);
+ }
+ return Status::OK();
+}
+
+// Calculations use canonical positions; raw spatial accessors still return
the original WKB.
+static S2Point canonical_geography_point(const S2Point& point) {
+ const auto coordinate = S2LatLng(point);
+ const double latitude = coordinate.lat().radians();
+ // All longitudes at a pole and +/-180 at the date line denote the same
position.
+ if (std::abs(latitude) == M_PI_2) {
+ return {0, 0, latitude > 0 ? 1.0 : -1.0};
+ }
+ if (coordinate.lng().radians() == M_PI) {
+ // Flip only sin(longitude): rebuilding latitude can change its last
bit.
+ return {point.x(), -point.y(), point.z()};
+ }
+ return point;
+}
+
+// Error-free TwoSum preserves tiny differences from a half-turn (e.g. -1e-15
+ 180).
+static std::pair<double, double> geography_degree_sum(double a, double b) {
+ const double high = a + b;
+ const double offset = high - a;
+ return {high, (a - (high - offset)) + (b - offset)};
+}
+
+// Reduce the exact degree sum before converting to radians. In particular,
+// 180 + (-1e-300) must retain its small sine instead of becoming sin(pi).
+static std::pair<long double, long double> geography_degree_sincos(double a,
double b = 0,
+ bool half =
false) {
+ const auto [high, low] = geography_degree_sum(a, b);
+ int quadrant;
+ const long double reduced =
+ std::remquo(static_cast<long double>(high), half ? 180.0L : 90.0L,
&quadrant) +
+ static_cast<long double>(low);
+ const auto radians = reduced * (std::numbers::pi_v<long double> / 180) /
(half ? 2 : 1);
+ const auto sine = std::sin(radians);
+ const auto cosine = std::cos(radians);
+ switch ((quadrant % 4 + 4) % 4) {
+ case 0:
+ return {sine, cosine};
+ case 1:
+ return {cosine, -sine};
+ case 2:
+ return {-sine, -cosine};
+ default:
+ return {-cosine, sine};
+ }
+}
+
+// Unlike cross(ToPoint(a), ToPoint(b)), these product identities do not
subtract
+// rounded nearly opposite vectors. Longitude-midpoint coordinates give
+// (-sin(dLon/2)*sin(latA+latB), cos(dLon/2)*sin(latA-latB),
+// cos(latA)*cos(latB)*sin(dLon)) for the oriented normal.
+struct GeographyRawArc {
+ S2Point normal;
+ S2Point midpoint;
+
+ GeographyRawArc(const GeoCoordinate& a, const GeoCoordinate& b) {
+ // Wide exponents retain products of the smallest binary64 degree
residuals.
+ static_assert(std::numeric_limits<long double>::digits >= 64 &&
+ std::numeric_limits<long double>::min_exponent <= -4096);
+ const auto [sd, cd] = geography_degree_sincos(b.x, -a.x, true);
+ const auto [sm, cm] = geography_degree_sincos(a.x, b.x, true);
+ const auto rotate_unit = [&](long double x, long double y, long double
z) {
+ const auto norm = std::sqrt(x * x + y * y + z * z);
+ DCHECK_GT(norm, 0);
+ return S2Point(static_cast<double>((cm * x - sm * y) / norm),
+ static_cast<double>((sm * x + cm * y) / norm),
+ static_cast<double>(z / norm));
+ };
+ normal = rotate_unit(-sd * geography_degree_sincos(a.y, b.y).first,
+ cd * geography_degree_sincos(a.y, -b.y).first,
+ geography_degree_sincos(a.y).second *
+ geography_degree_sincos(b.y).second *
+ geography_degree_sincos(b.x, -a.x).first);
+ const auto [sl, cl] = geography_degree_sincos(a.y, b.y, true);
+ const auto [st, ct] = geography_degree_sincos(b.y, -a.y, true);
+ midpoint = rotate_unit(cd * cl * ct, -sd * sl * st, sl * ct);
+ }
+};
+
+#ifdef BE_TEST
+static thread_local size_t geography_exact_proofs = 0;
+static thread_local size_t geography_candidate_visits = 0;
+#endif
+
+// Test coplanarity of the original degree ordinates, not independently
rounded S2 vectors.
+// For example (45,-45), (90,0), (135,45) share an arc although ToPoint's
determinant is nonzero.
+static bool geography_common_great_circle(const GeoCoordinate& a, const
GeoCoordinate& b,
+ const GeoCoordinate& c) {
+#ifdef BE_TEST
+ ++geography_exact_proofs;
+#endif
+ using boost::multiprecision::cpp_int;
+ const std::array<double, 6> ordinates {a.y, a.x, b.y, b.x, c.y, c.x};
+ std::array<uint64_t, 6> significands;
+ std::array<int, 6> exponents;
+ int denominator_bits = 0;
+ for (size_t i = 0; i < ordinates.size(); ++i) {
+ const auto bits = std::bit_cast<uint64_t>(ordinates[i]);
+ const int biased = (bits >> 52) & 0x7ff;
+ uint64_t significand = bits & ((uint64_t {1} << 52) - 1);
+ int exponent = biased ? biased - 1075 : -1074;
+ if (biased) {
+ significand |= uint64_t {1} << 52;
+ }
+ if (significand != 0) {
+ const int zeros = std::countr_zero(significand);
+ significand >>= zeros;
+ exponent += zeros;
+ denominator_bits = std::max(denominator_bits, -exponent);
+ } else {
+ exponent = 0;
+ }
+ significands[i] = significand;
+ exponents[i] = exponent;
+ }
+ std::array<cpp_int, 6> angles;
+ for (size_t i = 0; i < angles.size(); ++i) {
+ angles[i] = cpp_int(significands[i]) << (denominator_bits +
exponents[i]);
+ if (std::signbit(ordinates[i])) {
+ angles[i] = -angles[i];
+ }
+ }
+ // Binary-rational degrees have denominator 2^q. Exp(i*angle) is a power
of a
+ // root of unity of order 360*2^q. Its minimal polynomial is Phi_90(t),
with
+ // t = z^(2^(q+2)), and Phi_90(t) = t^24+t^21-t^15-t^12-t^9+t^3+1.
+ // Expand cos(lat_i)cos(lat_j)sin(lat_k)sin(lon_j-lon_i), cyclically, into
+ // 48 integer-coefficient powers. Equal residues modulo 2^(q+2) share one
+ // polynomial; reduction by Phi_90 proves zero exactly, including tiny
nonzero controls.
+ const int shift = denominator_bits + 2;
+ const cpp_int basis = cpp_int(1) << shift;
+ const cpp_int period = 90 * basis;
+ std::map<cpp_int, std::array<int, 90>> polynomials;
+ for (int i = 0; i < 3; ++i) {
+ const int j = (i + 1) % 3;
+ const int k = (i + 2) % 3;
+ for (int mask = 0; mask < 16; ++mask) {
+ std::array<int, 4> signs;
+ for (int bit = 0; bit < 4; ++bit) {
+ signs[bit] = (mask & (1 << bit)) ? 1 : -1;
+ }
+ cpp_int power = signs[0] * angles[2 * i] + signs[1] * angles[2 *
j] +
+ signs[2] * angles[2 * k] +
+ signs[3] * (angles[2 * j + 1] - angles[2 * i + 1]);
+ power %= period;
+ if (power < 0) {
+ power += period;
+ }
+ const int quotient = (power >> shift).convert_to<int>();
+ const cpp_int residue = power - (cpp_int(quotient) << shift);
+ polynomials[residue][quotient] += signs[2] * signs[3];
+ }
+ }
+ for (auto& [residue, coefficients] : polynomials) {
+ for (int i = 89; i >= 24; --i) {
+ const int coefficient = coefficients[i];
+ coefficients[i] = 0;
+ coefficients[i - 3] -= coefficient;
+ coefficients[i - 9] += coefficient;
+ coefficients[i - 12] += coefficient;
+ coefficients[i - 15] += coefficient;
+ coefficients[i - 21] -= coefficient;
+ coefficients[i - 24] -= coefficient;
+ }
+ if (std::ranges::any_of(coefficients, [](int value) { return value !=
0; })) {
+ return false;
+ }
+ }
+ return true;
+}
+
+// Keep each atomic segment in its own shape. Pinned S2 ProcessEdge1 uses a
chain's
+// global start as a local offset in one endpoint path, so multiple chains per
line
+// shape are unsafe. Closed AA shapes restore internal vertices in the OPEN
model.
+class GeographyLineShape final : public S2Shape {
+public:
+ GeographyLineShape(const S2Point& a, const S2Point& b) {
_edge.emplace_back(a, b); }
+ int num_edges() const override { return 1; }
+ Edge edge(int i) const override {
+ DCHECK_EQ(i, 0);
+ return _edge.front();
+ }
+ int dimension() const override { return 1; }
+ ReferencePoint GetReferencePoint() const override { return
ReferencePoint::Contained(false); }
+ int num_chains() const override { return 1; }
+ Chain chain(int i) const override {
+ DCHECK_EQ(i, 0);
+ return {0, 1};
+ }
+ Edge chain_edge(int i, int j) const override {
+ DCHECK_EQ(i, 0);
+ return edge(j);
+ }
+ ChainPosition chain_position(int i) const override {
+ DCHECK_EQ(i, 0);
+ return {0, 0};
+ }
+
+private:
+ // Owning storage is checked/accounted by the current query's allocator.
+ DorisVector<Edge> _edge;
+};
+
+class GeographyEdgeNodes {
+public:
+ GeographyEdgeNodes(StringRef left, StringRef right) {
+ GeoCoordinateList coordinates;
+ const auto left_status = WkbParse::xy_coordinates(left.data,
left.size, &coordinates);
+ const auto right_status = WkbParse::xy_coordinates(right.data,
right.size, &coordinates);
+ DCHECK_EQ(left_status, GeoParseStatus::GEO_PARSE_OK);
+ DCHECK_EQ(right_status, GeoParseStatus::GEO_PARSE_OK);
+ DorisMap<double, DorisVector<S2Point>> longitude_vertices;
+ for (auto coordinate : coordinates.list) {
+ if (coordinate.x == 180) {
+ coordinate.x = -180;
+ }
+ const auto point = canonical_geography_point(
+ S2LatLng::FromDegrees(coordinate.y,
coordinate.x).ToPoint());
+ if (!_coordinates.emplace(point, coordinate).second) {
+ continue;
+ }
+ longitude_vertices[coordinate.x].push_back(point);
+ if (std::abs(coordinate.y) != 90) {
+ _meridians[meridian(coordinate.x)][coordinate.x < 0].push_back(
+ {coordinate.y, point});
+ }
+ }
+ // Keep dense longitude groups separately so a circle can skip its
known
+ // endpoint longitudes before enumerating nearby points. At sqrt(N)
points
+ // per group there are at most sqrt(N) extra indexes, even with
unrelated
+ // vertices padding a near-antipodal spoke input.
+ const auto dense_size = std::max(
+ size_t {2},
static_cast<size_t>(std::ceil(std::sqrt(_coordinates.size()))));
+ for (const auto& [longitude, points] : longitude_vertices) {
+ auto& index =
+ points.size() >= dense_size ? _dense_longitudes[longitude]
: _vertex_index;
+ for (const auto& point : points) {
+ index.Add(point, 0);
+ }
+ }
+ for (auto& entry : _meridians) {
+ for (auto& side : entry.second) {
+ side.push_back({-90, {0, 0, -1}});
+ side.push_back({90, {0, 0, 1}});
+ std::ranges::sort(side, {}, &Vertex::latitude);
+ }
+ }
+ }
+
+ bool contains_point(const S2Point& raw_start, const S2Point& raw_end,
+ const S2Point& raw_point) const {
+ const auto start = canonical_geography_point(raw_start);
+ const auto end = canonical_geography_point(raw_end);
+ const auto point = canonical_geography_point(raw_point);
+ // Preserve exact endpoints, including date-line and pole aliases.
+ if (point == start || point == end) {
+ return true;
+ }
+ if (start == end) {
+ return false;
+ }
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const auto& c = _coordinates.at(point);
+ const auto key = meridian(std::abs(a.y) == 90 ? b.x : a.x);
+ bool member = false;
+ if (meridian_edge(start, end, [&](const auto& vertices, double from,
double to) {
+ // Opposite-longitude edges can pass through a pole; the pole
belongs
+ // to both legs, whereas other positions must be on the
selected side.
+ const bool same_side =
+ std::abs(c.y) == 90 ||
+ (meridian(c.x) == key && &_meridians.at(key)[c.x < 0]
== &vertices);
+ member |= same_side && c.y >= std::min(from, to) && c.y <=
std::max(from, to);
+ })) {
+ return member;
+ }
+ const bool wraps = geography_degree_sum(std::max(a.x, b.x),
-std::min(a.x, b.x)) >
+ std::pair {180.0, 0.0};
+ const bool within = wraps ? (c.x < std::min(a.x, b.x) || c.x >
std::max(a.x, b.x))
+ : (c.x > std::min(a.x, b.x) && c.x <
std::max(a.x, b.x));
+ // Near-antipodal endpoint conversion can rotate the S2 plane
substantially.
+ // POINT(45 0) is on (-44.99999999999997 -45,134.99999999999997 45);
+ // decide from raw degrees rather than accepting a larger distance
tolerance.
+ return within && geography_common_great_circle(a, b, c);
+ }
+
+ void append_edge(const S2Point& raw_start, const S2Point& raw_end,
+ DorisVector<S2Point>* output) const {
+ const auto start = canonical_geography_point(raw_start);
+ const auto end = canonical_geography_point(raw_end);
+ output->push_back(start);
+ if (start == end) {
+ return;
+ }
+ if (meridian_edge(start, end, [&](const auto& vertices, double from,
double to) {
+ append_range(vertices, from, to, output);
+ if (std::abs(to) == 90) {
+ const S2Point pole(0, 0, to > 0 ? 1 : -1);
+ if (end != pole) {
+ output->push_back(pole);
+ }
+ }
+ })) {
+ return;
+ }
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const bool wraps = geography_degree_sum(std::max(a.x, b.x),
-std::min(a.x, b.x)) >
+ std::pair {180.0, 0.0};
+ const bool east = (a.x < b.x) != wraps;
+ // A near-antipodal rounded edge can miss true raw-coordinate nodes.
Two
+ // stable half-arcs recall those nodes without a hemispherical search.
+ const auto midpoint = GeographyRawArc(a, b).midpoint;
+ DorisVector<S2Point> on_edge;
+ for (const auto& endpoint : {start, end}) {
+ visit_edge_candidates(
+ endpoint, midpoint, [](double) { return false; },
+ [&](const S2Point& point) {
+ const auto& c = _coordinates.at(point);
+ // Longitude is strictly monotone along every
non-meridian minor arc.
+ // Check the original interval (including date-line
wrap), not a distance
+ // tolerance that could append a point just beyond an
endpoint.
+ const bool within =
+ wraps ? (c.x < std::min(a.x, b.x) || c.x >
std::max(a.x, b.x))
+ : (c.x > std::min(a.x, b.x) && c.x <
std::max(a.x, b.x));
+ if (within && geography_common_great_circle(a, b, c)) {
+ on_edge.push_back(point);
+ }
+ });
+ }
+ std::ranges::sort(on_edge, [&](const S2Point& x, const S2Point& y) {
+ const double x_lon = _coordinates.at(x).x;
+ const double y_lon = _coordinates.at(y).x;
+ if ((x_lon < a.x) != (y_lon < a.x)) {
+ return east ? x_lon >= a.x : x_lon < a.x;
+ }
+ return east ? x_lon < y_lon : x_lon > y_lon;
+ });
+ on_edge.erase(std::ranges::unique(on_edge).begin(), on_edge.end());
+ output->insert(output->end(), on_edge.begin(), on_edge.end());
+ }
+
+ // Calculation-only points never enter the raw-degree identity/noding map.
+ std::optional<S2Point> arc_midpoint(const S2Point& a, const S2Point& b)
const {
+ // For (-45+e,-45)->(135-e,45), independently rounded endpoints rotate
+ // the plane. A raw midpoint makes both subedges shorter than pi/2.
+ if (a.DotProd(b) < 0) {
+ const auto coordinate = [&](const S2Point& point) -> GeoCoordinate
{
+ // Meridian noding can introduce a pole absent from the
original WKB.
+ if (std::abs(point.z()) == 1 && point.x() == 0 && point.y() ==
0) {
+ return {.x = 0, .y = point.z() > 0 ? 90.0 : -90.0};
+ }
+ return _coordinates.at(point);
+ };
+ return GeographyRawArc(coordinate(a), coordinate(b)).midpoint;
+ }
+ return std::nullopt;
+ }
+
+ void add_line(const GeoLine& shape, MutableS2ShapeIndex* index) const {
+ using Edge = std::pair<S2Point, S2Point>;
+ using EdgeSet = std::set<Edge, std::less<>, CustomStdAllocator<Edge>>;
+ using Interval = std::pair<double, double>;
+ DorisMap<const DorisVector<Vertex>*, DorisVector<Interval>> intervals;
+ EdgeSet raw_edges;
+ const auto* line = shape.polyline();
+ const auto first = canonical_geography_point(line->vertex(0));
+ const auto last =
canonical_geography_point(line->vertex(line->num_vertices() - 1));
+ for (int i = 0; i + 1 < line->num_vertices(); ++i) {
+ const auto a = canonical_geography_point(line->vertex(i));
+ const auto b = canonical_geography_point(line->vertex(i + 1));
+ if (a == b) {
+ continue;
+ }
+ if (!meridian_edge(a, b, [&](const auto& vertices, double from,
double to) {
+ intervals[&vertices].emplace_back(std::min(from, to),
std::max(from, to));
+ })) {
+ raw_edges.emplace(std::min(a, b), std::max(a, b));
+ }
+ }
+ EdgeSet edges;
+ auto add_edge = [&](const S2Point& a, const S2Point& b) {
+ edges.emplace(std::min(a, b), std::max(a, b));
+ };
+ // Merge latitude coverage before splitting: N repeated or overlapping
edges
+ // and M collinear vertices produce O(N+M) scratch, never an N*M
vertex walk.
+ for (auto& [vertices, ranges] : intervals) {
+ std::ranges::sort(ranges);
+ auto emit = [&](double from, double to) {
+ auto begin = std::ranges::lower_bound(*vertices, from, {},
&Vertex::latitude);
+ const auto end = std::ranges::upper_bound(*vertices, to, {},
&Vertex::latitude);
+ for (auto next = begin + 1; next != end; ++begin, ++next) {
+ add_edge(begin->point, next->point);
+ }
+ };
+ double from = ranges.front().first;
+ double to = ranges.front().second;
+ for (const auto& range : ranges) {
+ if (range.first > to) {
+ emit(from, to);
+ from = range.first;
+ }
+ to = std::max(to, range.second);
+ }
+ emit(from, to);
+ }
+ add_circle_edges(raw_edges, add_edge);
+ auto add = [&](const S2Point& a, const S2Point& b) {
+ index->Add(std::make_unique<GeographyLineShape>(a, b));
+ };
+ if (edges.empty()) {
+ // A raw pole-alias line may collapse to one canonical calculation
position.
+ add(first, first);
+ }
+ std::set<S2Point, std::less<>, CustomStdAllocator<S2Point>>
interior_vertices;
+ for (const auto& [a, b] : edges) {
+ const auto midpoint = arc_midpoint(a, b);
+ if (midpoint) {
+ add(a, *midpoint);
+ add(*midpoint, b);
+ // The synthetic split is interior, not a new OPEN line
boundary:
+ // a crossing at this point must not become ST_Touches=true.
+ interior_vertices.insert(*midpoint);
+ } else {
+ add(a, b);
+ }
+ interior_vertices.insert(a);
+ interior_vertices.insert(b);
+ }
+ for (const auto& point : interior_vertices) {
+ if (first == last || (point != first && point != last)) {
+ add(point, point);
+ }
+ }
+ }
+
+private:
+ struct Vertex {
+ double latitude;
+ S2Point point;
+ };
+ struct Circle {
+ S2Point a;
+ S2Point b;
+ S2Point normal;
+ DorisVector<std::pair<double, double>> ranges;
+ DorisMap<double, S2Point> vertices;
+ };
+
+ template <typename Edges, typename AddEdge>
+ void add_circle_edges(const Edges& edges, AddEdge&& add_edge) const {
+ DorisVector<Circle> circles;
+ S2PointIndex<size_t> normal_index;
+ const double rounding = 64 * std::numeric_limits<double>::epsilon();
+ for (const auto& [a, b] : edges) {
+ const auto normal = GeographyRawArc(_coordinates.at(a),
_coordinates.at(b)).normal;
+ size_t group = circles.size();
+ S2ClosestPointQuery<size_t> query(&normal_index);
+ // Raw trigonometric products avoid the epsilon/sin(length) error
+ // amplification. Normals only recall candidates; exact proof
decides.
+
query.mutable_options()->set_conservative_max_distance(S1Angle::Radians(rounding));
+ S2ClosestPointQuery<size_t>::PointTarget target(normal);
+ for (const auto& candidate : query.FindClosestPoints(&target)) {
+#ifdef BE_TEST
+ ++geography_candidate_visits;
+#endif
+ const auto& circle = circles[candidate.data()];
+ const auto& x = _coordinates.at(circle.a);
+ const auto& y = _coordinates.at(circle.b);
+ const auto on_circle = [&](const S2Point& point) {
+ return point == circle.a || point == circle.b ||
+ geography_common_great_circle(x, y,
_coordinates.at(point));
+ };
+ if (on_circle(a) && on_circle(b)) {
+ group = candidate.data();
+ break;
+ }
+ }
+ if (group == circles.size()) {
+ normal_index.Add(normal, group);
+ normal_index.Add(-normal, group);
+ circles.push_back({a, b, normal, {}, {}});
+ }
+ auto& circle = circles[group];
+ const double lo = std::min(_coordinates.at(a).x,
_coordinates.at(b).x);
+ const double hi = std::max(_coordinates.at(a).x,
_coordinates.at(b).x);
+ // Longitude cut points are interval sentinels, not fabricated
spatial
+ // vertices. For 170->-170 the real wrap edge is emitted below.
+ if (geography_degree_sum(hi, -lo) > std::pair {180.0, 0.0}) {
+ circle.ranges.emplace_back(-180, lo);
+ circle.ranges.emplace_back(hi, 180);
+ } else {
+ circle.ranges.emplace_back(lo, hi);
+ }
+ // Keep every raw endpoint, especially original open-line endpoints
+ // inside a larger merged interval: OPEN must still exclude those
points.
+ circle.vertices.emplace(_coordinates.at(a).x, a);
+ circle.vertices.emplace(_coordinates.at(b).x, b);
+ }
+ for (auto& circle : circles) {
+ collect_circle_vertices(circle);
+ emit_circle_edges(circle, add_edge);
+ }
+ }
+
+ template <typename SkipLongitude, typename Visit>
+ void visit_edge_candidates(const S2Point& a, const S2Point& b,
SkipLongitude&& skip_longitude,
+ Visit&& visit) const {
+ const auto query_index = [&](const S2PointIndex<int>& index) {
+ S2ClosestPointQuery<int> query(&index);
+ query.mutable_options()->set_conservative_max_distance(
+ S1Angle::Radians(64 *
std::numeric_limits<double>::epsilon()));
+ S2ClosestPointQuery<int>::EdgeTarget target(a, b);
+ for (const auto& candidate : query.FindClosestPoints(&target)) {
+#ifdef BE_TEST
+ ++geography_candidate_visits;
+#endif
+ visit(candidate.point());
+ }
+ };
+ query_index(_vertex_index);
+ for (const auto& [longitude, index] : _dense_longitudes) {
+#ifdef BE_TEST
+ ++geography_candidate_visits;
+#endif
+ // A non-meridian circle has only one nonpolar point at each
longitude.
+ // Dense endpoint groups must be skipped before the S2 query,
including
+ // when unrelated vertices make them a minority of the complete
input.
+ if (!skip_longitude(longitude)) {
+ query_index(index);
+ }
+ }
+ }
+
+ void collect_circle_vertices(Circle& circle) const {
+ // Query a whole circle once, using four non-antipodal quarter arcs.
+ const auto u = S2::Ortho(circle.normal);
+ const auto v = circle.normal.CrossProd(u).Normalize();
+ const std::array<S2Point, 4> quarter {u, v, -u, -v};
+ std::set<S2Point, std::less<>, CustomStdAllocator<S2Point>> seen;
+ for (size_t i = 0; i < quarter.size(); ++i) {
+ visit_edge_candidates(
+ quarter[i], quarter[(i + 1) % 4],
+ [&](double longitude) { return
circle.vertices.contains(longitude); },
+ [&](const S2Point& point) {
+ const auto& coordinate = _coordinates.at(point);
+ if (seen.insert(point).second &&
!circle.vertices.contains(coordinate.x) &&
+
geography_common_great_circle(_coordinates.at(circle.a),
+
_coordinates.at(circle.b), coordinate)) {
+ circle.vertices.emplace(coordinate.x, point);
+ }
+ });
+ }
+ }
+
+ template <typename AddEdge>
+ static void emit_circle_edges(Circle& circle, AddEdge&& add_edge) {
+ std::ranges::sort(circle.ranges);
+ DorisVector<std::pair<double, double>> merged;
+ for (const auto& range : circle.ranges) {
+ if (merged.empty() || range.first > merged.back().second) {
+ merged.push_back(range);
+ } else {
+ merged.back().second = std::max(merged.back().second,
range.second);
+ }
+ }
+ for (const auto& [lo, hi] : merged) {
+ const auto begin = circle.vertices.lower_bound(lo);
+ const auto end = circle.vertices.upper_bound(hi);
+ // Each original interval retains an endpoint, including either
side
+ // of the date-line cut. Merging and node collection never remove
it.
+ DCHECK(begin != end);
+ auto previous = begin;
+ for (auto next = std::next(begin); next != end; ++previous,
++next) {
+ add_edge(previous->second, next->second);
+ }
+ }
+ if (merged.front().first == -180 && merged.back().second == 180) {
+ add_edge(circle.vertices.rbegin()->second,
circle.vertices.begin()->second);
+ }
+ }
+
+ static std::pair<double, double> meridian(double longitude) {
+ return longitude >= 0 ? std::pair {longitude, 0.0} :
geography_degree_sum(longitude, 180);
+ }
+ template <typename Append>
+ bool meridian_edge(const S2Point& start, const S2Point& end, Append&&
append) const {
+ const auto& a = _coordinates.at(start);
+ const auto& b = _coordinates.at(end);
+ const bool a_pole = std::abs(a.y) == 90;
+ const bool b_pole = std::abs(b.y) == 90;
+ const auto key = meridian(a_pole ? b.x : a.x);
+ if (!a_pole && !b_pole && key != meridian(b.x)) {
+ return false;
+ }
+ // Exact opposite endpoints (e.g. (0,10),(180,-10)) have no unique
minor arc.
+ if (a.y == -b.y && ((a_pole && b_pole) || (a.x < 0) != (b.x < 0))) {
+ throw Exception(ErrorCode::NOT_IMPLEMENTED_ERROR,
+ "Geography calculations do not support antipodal
edges");
+ }
+ const bool a_side = (a_pole ? b.x : a.x) < 0;
+ const bool b_side = (b_pole ? a.x : b.x) < 0;
+ const auto& group = _meridians.at(key);
+ if (a_side == b_side) {
+ append(group[a_side], a.y, b.y);
+ } else {
+ // (10,80)->(-170,80) goes via the north pole; retain raw latitude
on each side.
+ const double latitude = a.y > -b.y ? 90 : -90;
+ append(group[a_side], a.y, latitude);
+ append(group[b_side], latitude, b.y);
+ }
+ return true;
+ }
+ static void append_range(const DorisVector<Vertex>& vertices, double from,
double to,
+ DorisVector<S2Point>* output) {
+ const auto first =
+ std::ranges::upper_bound(vertices, std::min(from, to), {},
&Vertex::latitude);
+ const auto last =
+ std::ranges::lower_bound(vertices, std::max(from, to), {},
&Vertex::latitude);
+ if (from < to) {
+ for (auto it = first; it != last; ++it) {
+ output->push_back(it->point);
+ }
+ } else {
+ auto current = last;
+ while (current != first) {
+ output->push_back((--current)->point);
+ }
+ }
+ }
+ DorisMap<S2Point, GeoCoordinate> _coordinates;
+ DorisMap<double, S2PointIndex<int>> _dense_longitudes;
+ S2PointIndex<int> _vertex_index;
+ DorisMap<std::pair<double, double>, std::array<DorisVector<Vertex>, 2>>
_meridians;
+};
+
+static void geography_shape_index(const GeoShape& shape, const
GeographyEdgeNodes& nodes,
+ MutableS2ShapeIndex* index) {
+ switch (shape.type()) {
+ case GeoShapeType::GEO_SHAPE_POINT:
+ index->Add(std::make_unique<S2PointVectorShape>(std::vector<S2Point> {
+ canonical_geography_point(*assert_cast<const
GeoPoint&>(shape).point())}));
+ break;
+ case GeoShapeType::GEO_SHAPE_LINE_STRING:
+ nodes.add_line(assert_cast<const GeoLine&>(shape), index);
+ break;
+ case GeoShapeType::GEO_SHAPE_POLYGON: {
+ const auto* polygon = assert_cast<const GeoPolygon&>(shape).polygon();
+ std::vector<S2LaxPolygonShape::Loop> loops;
+ for (int i = 0; i < polygon->num_loops(); ++i) {
+ auto& vertices = loops.emplace_back();
+ const auto* loop = polygon->loop(i);
+ DorisVector<S2Point> edge;
+ for (int j = 0; j < loop->num_vertices(); ++j) {
+ edge.clear();
+ nodes.append_edge(loop->vertex(j), loop->vertex(j + 1), &edge);
+ edge.push_back(canonical_geography_point(loop->vertex(j + 1)));
+ for (size_t k = 0; k + 1 < edge.size(); ++k) {
+ vertices.push_back(edge[k]);
+ // Use the same stable long arcs for interiors and
distances.
+ // The raw (-45+e,-45)->(135-e,45) edge otherwise changes
sides
+ // around POINT(45 0.000001) when its midpoint is inserted.
+ if (const auto midpoint = nodes.arc_midpoint(edge[k],
edge[k + 1])) {
+ vertices.push_back(*midpoint);
+ }
+ }
+ }
+ // A thin near-antipodal ring may have been reversed by
normalization
+ // of the rounded edge. Restore the smaller-region shell
orientation
+ // after stabilization while preserving the admitted hole depths.
+ // Span normalization also retains canonical pole-alias
degeneracies.
+ if (S2::IsNormalized(vertices) == loop->is_hole()) {
+ std::ranges::reverse(vertices);
+ }
+ }
+ index->Add(std::make_unique<S2LaxPolygonShape>(loops));
+ break;
+ }
+ default:
+ throw Exception(ErrorCode::INTERNAL_ERROR, "Unexpected Geography
calculation shape");
+ }
+}
+
+#ifdef BE_TEST
+size_t geography_calculation_edges_for_test(const std::string& left, const
std::string& right,
+ size_t* exact_proofs, size_t*
candidate_visits) {
+ geography_exact_proofs = 0;
+ geography_candidate_visits = 0;
+ GeoParseStatus status;
+ const auto a = GeoShape::from_wkb_bytes(left.data(), left.size(), status);
+ const auto b = GeoShape::from_wkb_bytes(right.data(), right.size(),
status);
+ const GeographyEdgeNodes nodes {StringRef(left), StringRef(right)};
+ MutableS2ShapeIndex a_index;
+ MutableS2ShapeIndex b_index;
+ geography_shape_index(*a, nodes, &a_index);
+ geography_shape_index(*b, nodes, &b_index);
+ size_t edges = 0;
+ for (const auto* shape : a_index) {
+ edges += shape->num_edges();
+ }
+ for (const auto* shape : b_index) {
+ edges += shape->num_edges();
+ }
+ if (exact_proofs != nullptr) {
+ *exact_proofs = geography_exact_proofs;
+ }
+ if (candidate_visits != nullptr) {
+ *candidate_visits = geography_candidate_visits;
+ }
+ return edges;
+}
+#endif
+
+static S1ChordAngle geography_conversion_error() {
+ // S2LatLng::ToPoint and degrees-to-radians rounding can separate
positions on
+ // the same meridian. This is about 1e-8 metres, not a user snap radius.
+ return S1ChordAngle(S1Angle::Radians(8 *
std::numeric_limits<double>::epsilon()));
+}
+
+static bool geography_point_relation(std::string_view name, const GeoShape&
lhs,
+ const GeoShape& rhs, const
GeographyEdgeNodes& nodes) {
+ const auto& point_shape = lhs.type() == GeoShapeType::GEO_SHAPE_POINT ?
lhs : rhs;
+ const auto& other = lhs.type() == GeoShapeType::GEO_SHAPE_POINT ? rhs :
lhs;
+ const auto point =
+ canonical_geography_point(*assert_cast<const
GeoPoint&>(point_shape).point());
+ bool interior = false;
+ bool boundary = false;
+ if (other.type() == GeoShapeType::GEO_SHAPE_POINT) {
+ interior = point == canonical_geography_point(*assert_cast<const
GeoPoint&>(other).point());
+ } else if (other.type() == GeoShapeType::GEO_SHAPE_LINE_STRING) {
+ const auto* line = assert_cast<const GeoLine&>(other).polyline();
+ bool member = false;
+ for (int i = 0; i + 1 < line->num_vertices(); ++i) {
+ member = member || nodes.contains_point(line->vertex(i),
line->vertex(i + 1), point);
+ }
+ const auto first = canonical_geography_point(line->vertex(0));
+ const auto last =
canonical_geography_point(line->vertex(line->num_vertices() - 1));
+ boundary = first != last && (point == first || point == last);
+ interior = member && !boundary;
+ } else {
+ DCHECK(other.type() == GeoShapeType::GEO_SHAPE_POLYGON);
+ const auto* polygon = assert_cast<const GeoPolygon&>(other).polygon();
+ // Polygon rings have closing edges; open polylines above do not.
+ for (int i = 0; i < polygon->num_loops(); ++i) {
+ const auto* loop = polygon->loop(i);
+ for (int j = 0; j < loop->num_vertices(); ++j) {
+ boundary = boundary ||
+ nodes.contains_point(loop->vertex(j),
loop->vertex(j + 1), point);
+ }
+ }
+ // Raw boundary membership above stays authoritative. Off-boundary
points
+ // need the stable polygon too, not the original rounded
near-antipodal loop.
+ if (!boundary) {
+ MutableS2ShapeIndex index;
+ geography_shape_index(other, nodes, &index);
+ interior = MakeS2ContainsPointQuery(&index).Contains(point);
+ }
+ }
+ if (name == "st_intersects") {
+ return interior || boundary;
+ }
+ if (name == "st_disjoint") {
+ return !interior && !boundary;
+ }
+ if (name == "st_touches") {
+ return boundary;
+ }
+ return interior;
+}
+
+static Status geography_relation(std::string_view name, const GeoShape& lhs,
const GeoShape& rhs,
+ StringRef left_wkb, StringRef right_wkb,
bool* result) {
+ const auto& container = name == "st_within" ? rhs : lhs;
+ const auto& value = name == "st_within" ? lhs : rhs;
+ const bool containment = name == "st_contains" || name == "st_within";
+ // A point cannot contain a line or polygon; this is a defined false
result.
+ if (containment && container.type() == GeoShapeType::GEO_SHAPE_POINT &&
+ value.type() != GeoShapeType::GEO_SHAPE_POINT) {
+ *result = false;
+ return Status::OK();
+ }
+ if (name == "st_equals" && lhs.type() != rhs.type()) {
+ *result = false;
+ return Status::OK();
+ }
+ if (name == "st_overlaps" &&
+ (lhs.type() != rhs.type() || lhs.type() ==
GeoShapeType::GEO_SHAPE_POINT)) {
+ *result = false;
+ return Status::OK();
+ }
+ if (lhs.type() == GeoShapeType::GEO_SHAPE_POINT ||
+ rhs.type() == GeoShapeType::GEO_SHAPE_POINT) {
+ // Point equality has the same canonical position test as point
containment.
+ const GeographyEdgeNodes nodes(left_wkb, right_wkb);
+ *result = geography_point_relation(name, lhs, rhs, nodes);
+ return Status::OK();
+ }
+ MutableS2ShapeIndex left;
+ MutableS2ShapeIndex right;
+ const GeographyEdgeNodes nodes(left_wkb, right_wkb);
+ geography_shape_index(container, nodes, &left);
+ geography_shape_index(value, nodes, &right);
+ S2BooleanOperation::Options options;
+ options.set_polygon_model(S2BooleanOperation::PolygonModel::CLOSED);
+ options.set_polyline_model(S2BooleanOperation::PolylineModel::CLOSED);
+ options.set_polyline_loops_have_boundaries(false);
+ if (name == "st_equals") {
+ *result = S2BooleanOperation::Equals(left, right, options);
+ return Status::OK();
+ }
+ if (containment && container.type() == GeoShapeType::GEO_SHAPE_LINE_STRING
&&
+ value.type() == GeoShapeType::GEO_SHAPE_POLYGON) {
+ *result = false;
+ return Status::OK();
+ }
+ if (containment && !S2BooleanOperation::Contains(left, right, options)) {
+ *result = false;
+ return Status::OK();
+ }
+ auto interior_options = options;
+ interior_options.set_polygon_model(S2BooleanOperation::PolygonModel::OPEN);
+
interior_options.set_polyline_model(S2BooleanOperation::PolylineModel::OPEN);
+ if (containment) {
+ // CLOSED Contains is covers: a polygon must not contain a line solely
on its boundary.
+ *result = S2BooleanOperation::Intersects(left, right,
interior_options);
+ return Status::OK();
+ }
+ if (name == "st_overlaps") {
+ // Neither equality nor containment is a partial overlap.
+ if (S2BooleanOperation::Contains(left, right, options) ||
+ S2BooleanOperation::Contains(right, left, options)) {
+ *result = false;
+ return Status::OK();
+ }
+ if (lhs.type() == GeoShapeType::GEO_SHAPE_POLYGON) {
+ *result = S2BooleanOperation::Intersects(left, right,
interior_options);
+ return Status::OK();
+ }
+ // S2 emits a crossing point as a degenerate polyline: only a real
shared edge overlaps.
+ std::vector<std::unique_ptr<S2Polyline>> intersection;
+ S2BooleanOperation operation(
+ S2BooleanOperation::OpType::INTERSECTION,
+
std::make_unique<s2builderutil::S2PolylineVectorLayer>(&intersection),
+ interior_options);
+ S2Error error;
+ if (!operation.Build(left, right, &error)) {
+ return Status::InternalError("Geography intersection failed: {}",
error.text());
+ }
+ *result = std::ranges::any_of(intersection, [](const auto& line) {
+ for (int i = 1; i < line->num_vertices(); ++i) {
+ if (line->vertex(i - 1) != line->vertex(i)) {
+ return true;
+ }
+ }
+ return false;
+ });
+ return Status::OK();
+ }
+ bool intersects = S2BooleanOperation::Intersects(left, right, options);
+ // A nearby disjoint curve must not become a boundary contact through the
legacy
+ // conversion-distance fallback used below by Intersects/Disjoint.
+ if (name == "st_touches") {
+ *result = intersects && !S2BooleanOperation::Intersects(left, right,
interior_options);
+ return Status::OK();
+ }
+ // Independently converted LINESTRING(10 45,10 46) and (10 45.3,10 45.9)
+ // are not exactly coplanar either. Apply the same conversion bound to all
+ // admitted nonpoint pairs, retaining polygon interiors and real open-line
edges.
+ if (!intersects) {
+ S2ClosestEdgeQuery query(&left);
+ S2ClosestEdgeQuery::ShapeIndexTarget target(&right);
+ intersects = query.IsConservativeDistanceLessOrEqual(&target,
geography_conversion_error());
+ }
+ DCHECK(name == "st_intersects" || name == "st_disjoint");
+ *result = name == "st_intersects" ? intersects : !intersects;
+ return Status::OK();
+}
+
+static std::pair<double, double> geography_point_coordinate(StringRef wkb) {
+ double longitude;
+ double latitude;
+ const auto status = WkbParse::point_coordinates(wkb.data, wkb.size,
&longitude, &latitude);
+ DCHECK_EQ(status, GeoParseStatus::GEO_PARSE_OK);
+ // Use original coordinates: reconstructing a shared longitude from two S2
vectors
+ // can turn due north into a tiny negative bearing that wraps to almost
2*pi.
+ if (std::abs(latitude) == 90) {
+ longitude = 0;
+ } else if (longitude == 180) {
+ longitude = -180;
+ }
+ return {longitude, latitude};
+}
+
+static bool geography_bearing(StringRef from, StringRef to, double* bearing) {
+ const auto start = geography_point_coordinate(from);
+ const auto end = geography_point_coordinate(to);
+ // Coincident or antipodal positions have no unique initial bearing.
Compare original degrees:
+ // haversine distance rounds (0, 0) -> (179.999999, 0) to pi despite its
unique east bearing.
+ const bool antipodal =
+ start.second == -end.second &&
+ (std::abs(start.second) == 90 || std::abs(start.first - end.first)
== 180);
+ if (start == end || antipodal) {
+ return false;
+ }
+ *bearing = S2Earth::GetInitialBearing(S2LatLng::FromDegrees(start.second,
start.first),
Review Comment:
[P2] Stabilize the near-antipodal bearing calculation
For default GEOGRAPHY A=`POINT(0 45)` and B=`POINT(180 -44.99999999999999)`
(B's latitude is `nextafter(-45,+infinity)`), the unique shorter great-circle
arc goes due north through the north pole, so `ST_Azimuth(A,B)` must be 0. The
antipodal guard correctly admits this pair because the latitudes are not
opposites. However, after +180 is canonicalized to -180, the pinned [S2 bearing
formula](https://github.com/google/s2geometry/blob/v0.10.0/src/s2/s2earth.cc#L43-L54)
cancels its north component to `x=0` while retaining `y=sin(-pi)*cos(lat2)<0`;
a source/arithmetic model gives `3*pi/2` (west). With N=`POINT(0 46)`, shared
`ST_Angle(N,A,B)` likewise becomes `3*pi/2` instead of 0. No Doris runtime test
was run. This is the numeric computation after admission, separate from
existing 4290141166's antipode guard/false-NULL issue and 4269360833's repaired
legacy axis/order defect. Use stable original-degree tangent/arc-normal
arithmetic here and cover near-antipodal meridian directions a
nd the shared angle path.
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: [email protected]
For queries about this service, please contact Infrastructure at:
[email protected]
---------------------------------------------------------------------
To unsubscribe, e-mail: [email protected]
For additional commands, e-mail: [email protected]
