[postgis-tickets] r17706 - Geography Distance inconsistent with Intersects
Paul Ramsey
pramsey at cleverelephant.ca
Tue Aug 13 09:59:20 PDT 2019
Author: pramsey
Date: 2019-08-13 09:59:20 -0700 (Tue, 13 Aug 2019)
New Revision: 17706
Modified:
trunk/liblwgeom/cunit/cu_geodetic.c
trunk/liblwgeom/gbox.c
trunk/liblwgeom/liblwgeom.h.in
trunk/liblwgeom/lwgeodetic.c
trunk/liblwgeom/lwgeodetic.h
trunk/liblwgeom/lwgeodetic_tree.c
trunk/liblwgeom/lwgeodetic_tree.h
trunk/postgis/geography_measurement_trees.c
Log:
Geography Distance inconsistent with Intersects
Closes #4480
Modified: trunk/liblwgeom/cunit/cu_geodetic.c
===================================================================
--- trunk/liblwgeom/cunit/cu_geodetic.c 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/cunit/cu_geodetic.c 2019-08-13 16:59:20 UTC (rev 17706)
@@ -1135,7 +1135,7 @@
lwgeom_free(lwg);
/* Great big ring */
- lwg = lwgeom_from_wkt("POLYGON((-40.0 52.0, 102.0 -6.0, -67.0 -29.0, -40.0 52.0))", LW_PARSER_CHECK_NONE);
+ lwg = lwgeom_from_wkt("POLYGON((-40.0 52.0, -67.0 -29.0, 102.0 -6.0, -40.0 52.0))", LW_PARSER_CHECK_NONE);
poly = (LWPOLY*)lwg;
pt_to_test.x = 4.0;
pt_to_test.y = 11.0;
Modified: trunk/liblwgeom/gbox.c
===================================================================
--- trunk/liblwgeom/gbox.c 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/gbox.c 2019-08-13 16:59:20 UTC (rev 17706)
@@ -247,7 +247,7 @@
int gbox_contains_point3d(const GBOX *gbox, const POINT3D *pt)
{
if ( gbox->xmin > pt->x || gbox->ymin > pt->y || gbox->zmin > pt->z ||
- gbox->xmax < pt->x || gbox->ymax < pt->y || gbox->zmax < pt->z )
+ gbox->xmax < pt->x || gbox->ymax < pt->y || gbox->zmax < pt->z )
{
return LW_FALSE;
}
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/liblwgeom.h.in 2019-08-13 16:59:20 UTC (rev 17706)
@@ -1921,7 +1921,7 @@
/**
* Calculate a spherical point that falls outside the geocentric gbox
*/
-void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside);
+int gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside);
/**
* Create a new gbox with the dimensionality indicated by the flags. Caller
Modified: trunk/liblwgeom/lwgeodetic.c
===================================================================
--- trunk/liblwgeom/lwgeodetic.c 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/lwgeodetic.c 2019-08-13 16:59:20 UTC (rev 17706)
@@ -323,34 +323,49 @@
lwfree(gbox_str);
#endif
/* Z axis */
- if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
- gbox->ymin < 0.0 && gbox->ymax > 0.0 )
+ if (gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
+ gbox->ymin < 0.0 && gbox->ymax > 0.0)
{
- if ( (gbox->zmin + gbox->zmax) > 0.0 )
+ /* Extrema lean positive */
+ if ((gbox->zmin > 0.0) && (gbox->zmax > 0.0))
{
LWDEBUG(4, "enclosed positive z axis");
gbox->zmax = 1.0;
}
- else
+ /* Extrema lean negative */
+ else if ((gbox->zmin < 0.0) && (gbox->zmax < 0.0))
{
LWDEBUG(4, "enclosed negative z axis");
gbox->zmin = -1.0;
}
+ /* Extrema both sides! */
+ else
+ {
+ LWDEBUG(4, "enclosed both z axes");
+ gbox->zmin = -1.0;
+ gbox->zmax = 1.0;
+ }
rv = LW_TRUE;
}
/* Y axis */
- if ( gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
- gbox->zmin < 0.0 && gbox->zmax > 0.0 )
+ if (gbox->xmin < 0.0 && gbox->xmax > 0.0 &&
+ gbox->zmin < 0.0 && gbox->zmax > 0.0)
{
- if ( gbox->ymin + gbox->ymax > 0.0 )
+ if ((gbox->ymin > 0.0) && (gbox->ymax > 0.0))
{
LWDEBUG(4, "enclosed positive y axis");
gbox->ymax = 1.0;
}
+ else if ((gbox->ymin < 0.0) && (gbox->ymax < 0.0))
+ {
+ LWDEBUG(4, "enclosed negative y axis");
+ gbox->ymax = -1.0;
+ }
else
{
- LWDEBUG(4, "enclosed negative y axis");
+ LWDEBUG(4, "enclosed both y axes");
+ gbox->ymax = 1.0;
gbox->ymin = -1.0;
}
rv = LW_TRUE;
@@ -357,19 +372,26 @@
}
/* X axis */
- if ( gbox->ymin < 0.0 && gbox->ymax > 0.0 &&
- gbox->zmin < 0.0 && gbox->zmax > 0.0 )
+ if (gbox->ymin < 0.0 && gbox->ymax > 0.0 &&
+ gbox->zmin < 0.0 && gbox->zmax > 0.0)
{
- if ( gbox->xmin + gbox->xmax > 0.0 )
+ if ((gbox->xmin > 0.0) && (gbox->xmax > 0.0))
{
LWDEBUG(4, "enclosed positive x axis");
gbox->xmax = 1.0;
}
- else
+ else if ((gbox->xmin < 0.0) && (gbox->xmax < 0.0))
{
LWDEBUG(4, "enclosed negative x axis");
gbox->xmin = -1.0;
}
+ else
+ {
+ LWDEBUG(4, "enclosed both x axes");
+ gbox->xmax = 1.0;
+ gbox->xmin = -1.0;
+ }
+
rv = LW_TRUE;
}
@@ -462,7 +484,7 @@
/**
* Scale a vector out by a factor
*/
-static void vector_scale(POINT3D *n, double scale)
+void vector_scale(POINT3D *n, double scale)
{
n->x *= scale;
n->y *= scale;
@@ -1450,21 +1472,77 @@
return LW_SUCCESS;
}
-void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
+/*
+* When we have a globe-covering gbox but we still want an outside
+* point, we do this Very Bad Hack, which is look at the first two points
+* in the ring and then nudge a point to the left of that arc.
+* There is an assumption of convexity built in there, as well as that
+* the shape doesn't have a sharp reversal in it. It's ugly, but
+* it fixes some common cases (large selection polygons) that users
+* are generating. At some point all of geodetic needs a clean-room
+* rewrite.
+* There is also an assumption of CCW exterior ring, which is how the
+* GeoJSON spec defined geographic ring orientation.
+*/
+static int lwpoly_pt_outside_hack(const LWPOLY *poly, POINT2D *pt_outside)
{
+ GEOGRAPHIC_POINT g1, g2, gSum;
+ POINT4D p1, p2;
+ POINT3D q1, q2, qMid, qCross, qSum;
+ POINTARRAY *pa;
+ if (lwgeom_is_empty((LWGEOM*)poly))
+ return LW_FAILURE;
+ if (poly->nrings < 1)
+ return LW_FAILURE;
+ pa = poly->rings[0];
+ if (pa->npoints < 2)
+ return LW_FAILURE;
+
+ /* First two points of ring */
+ getPoint4d_p(pa, 0, &p1);
+ getPoint4d_p(pa, 1, &p2);
+ /* Convert to XYZ unit vectors */
+ geographic_point_init(p1.x, p1.y, &g1);
+ geographic_point_init(p2.x, p2.y, &g2);
+ geog2cart(&g1, &q1);
+ geog2cart(&g2, &q2);
+ /* Mid-point of first two points */
+ vector_sum(&q1, &q2, &qMid);
+ normalize(&qMid);
+ /* Cross product of first two points (perpendicular) */
+ cross_product(&q1, &q2, &qCross);
+ normalize(&qCross);
+ /* Invert it to put it outside, and scale down */
+ vector_scale(&qCross, -0.2);
+ /* Project midpoint to the right */
+ vector_sum(&qMid, &qCross, &qSum);
+ normalize(&qSum);
+ /* Convert back to lon/lat */
+ cart2geog(&qSum, &gSum);
+ pt_outside->x = rad2deg(gSum.lon);
+ pt_outside->y = rad2deg(gSum.lat);
+ return LW_SUCCESS;
+}
+
+int lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside)
+{
+ int rv;
/* Make sure we have boxes */
if ( poly->bbox )
{
- gbox_pt_outside(poly->bbox, pt_outside);
- return;
+ rv = gbox_pt_outside(poly->bbox, pt_outside);
}
else
{
GBOX gbox;
lwgeom_calculate_gbox_geodetic((LWGEOM*)poly, &gbox);
- gbox_pt_outside(&gbox, pt_outside);
- return;
+ rv = gbox_pt_outside(&gbox, pt_outside);
}
+
+ if (rv == LW_FALSE)
+ return lwpoly_pt_outside_hack(poly, pt_outside);
+
+ return rv;
}
/**
@@ -1471,7 +1549,7 @@
* Given a unit geocentric gbox, return a lon/lat (degrees) coordinate point point that is
* guaranteed to be outside the box (and therefore anything it contains).
*/
-void gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
+int gbox_pt_outside(const GBOX *gbox, POINT2D *pt_outside)
{
double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
int i;
@@ -1538,7 +1616,7 @@
pt_outside->x = rad2deg(g.lon);
pt_outside->y = rad2deg(g.lat);
LWDEBUGF(4, "returning POINT(%.8g %.8g) as outside point", pt_outside->x, pt_outside->y);
- return;
+ return LW_SUCCESS;
}
}
@@ -1547,8 +1625,8 @@
}
/* This should never happen! */
- lwerror("BOOM! Could not generate outside point!");
- return;
+ // lwerror("BOOM! Could not generate outside point!");
+ return LW_FAILURE;
}
@@ -2473,7 +2551,7 @@
}
/* Calculate our outside point from the gbox */
- gbox_pt_outside(&gbox, &pt_outside);
+ lwpoly_pt_outside(poly, &pt_outside);
LWDEBUGF(4, "pt_outside POINT(%.18g %.18g)", pt_outside.x, pt_outside.y);
LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test->x, pt_to_test->y);
@@ -2807,6 +2885,7 @@
return LW_SUCCESS;
}
+
int ptarray_calculate_gbox_geodetic(const POINTARRAY *pa, GBOX *gbox)
{
uint32_t i;
Modified: trunk/liblwgeom/lwgeodetic.h
===================================================================
--- trunk/liblwgeom/lwgeodetic.h 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/lwgeodetic.h 2019-08-13 16:59:20 UTC (rev 17706)
@@ -126,7 +126,7 @@
int lwline_covers_lwline(const LWLINE* lwline1, const LWLINE* lwline2);
int lwline_covers_lwpoint(const LWLINE* lwline, const LWPOINT* lwpoint);
int lwpoly_intersects_line(const LWPOLY* lwpoly, const POINTARRAY* line);
-void lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside);
+int lwpoly_pt_outside(const LWPOLY *poly, POINT2D *pt_outside);
int ptarray_point_in_ring(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test);
double ptarray_area_sphere(const POINTARRAY *pa);
double latitude_degrees_normalize(double lat);
@@ -138,6 +138,7 @@
double longitude_radians_normalize(double lon);
double latitude_radians_normalize(double lat);
void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n);
+void vector_scale(POINT3D *a, double s);
double vector_angle(const POINT3D* v1, const POINT3D* v2);
void vector_rotate(const POINT3D* v1, const POINT3D* v2, double angle, POINT3D* n);
void normalize(POINT3D *p);
Modified: trunk/liblwgeom/lwgeodetic_tree.c
===================================================================
--- trunk/liblwgeom/lwgeodetic_tree.c 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/lwgeodetic_tree.c 2019-08-13 16:59:20 UTC (rev 17706)
@@ -472,7 +472,20 @@
}
}
+int circ_tree_get_point_outside(const CIRC_NODE* node, POINT2D* pt)
+{
+ POINT3D center3d;
+ GEOGRAPHIC_POINT g;
+ if (node->radius >= M_PI) return LW_FAILURE;
+ geog2cart(&(node->center), ¢er3d);
+ vector_scale(¢er3d, -1.0);
+ cart2geog(¢er3d, &g);
+ pt->x = rad2deg(g.lon);
+ pt->y = rad2deg(g.lat);
+ return LW_SUCCESS;
+}
+
/**
* Walk the tree and count intersections between the stab line and the edges.
* odd => containment, even => no containment.
Modified: trunk/liblwgeom/lwgeodetic_tree.h
===================================================================
--- trunk/liblwgeom/lwgeodetic_tree.h 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/liblwgeom/lwgeodetic_tree.h 2019-08-13 16:59:20 UTC (rev 17706)
@@ -54,6 +54,7 @@
double circ_tree_distance_tree(const CIRC_NODE* n1, const CIRC_NODE* n2, const SPHEROID *spheroid, double threshold);
CIRC_NODE* lwgeom_calculate_circ_tree(const LWGEOM* lwgeom);
int circ_tree_get_point(const CIRC_NODE* node, POINT2D* pt);
+int circ_tree_get_point_outside(const CIRC_NODE* node, POINT2D* pt);
#endif /* _LWGEODETIC_TREE_H */
Modified: trunk/postgis/geography_measurement_trees.c
===================================================================
--- trunk/postgis/geography_measurement_trees.c 2019-08-13 16:59:08 UTC (rev 17705)
+++ trunk/postgis/geography_measurement_trees.c 2019-08-13 16:59:20 UTC (rev 17706)
@@ -95,7 +95,6 @@
return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2);
}
-
static int
CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const POINT4D* in_point)
{
@@ -137,7 +136,10 @@
pt2d_inside.x = in_point->x;
pt2d_inside.y = in_point->y;
/* Calculate a definitive outside point */
- gbox_pt_outside(&gbox1, &pt2d_outside);
+ if (gbox_pt_outside(&gbox1, &pt2d_outside) == LW_FAILURE)
+ if (circ_tree_get_point_outside(tree1, &pt2d_outside) == LW_FAILURE)
+ lwpgerror("%s: Unable to generate outside point!", __func__);
+
POSTGIS_DEBUGF(3, "p2d_inside=POINT(%g %g) p2d_outside=POINT(%g %g)", pt2d_inside.x, pt2d_inside.y, pt2d_outside.x, pt2d_outside.y);
/* Test the candidate point for strict containment */
POSTGIS_DEBUG(3, "calling circ_tree_contains_point for PiP test");
More information about the postgis-tickets
mailing list