[postgis-tickets] r17704 - branch 2.4: Geography Distance inconsistent with Intersects
Paul Ramsey
pramsey at cleverelephant.ca
Tue Aug 13 09:58:58 PDT 2019
Author: pramsey
Date: 2019-08-13 09:58:58 -0700 (Tue, 13 Aug 2019)
New Revision: 17704
Modified:
branches/2.4/NEWS
branches/2.4/liblwgeom/cunit/cu_geodetic.c
branches/2.4/liblwgeom/g_box.c
branches/2.4/liblwgeom/liblwgeom.h.in
branches/2.4/liblwgeom/lwgeodetic.c
branches/2.4/liblwgeom/lwgeodetic.h
branches/2.4/liblwgeom/lwgeodetic_tree.c
branches/2.4/liblwgeom/lwgeodetic_tree.h
branches/2.4/postgis/geography_measurement_trees.c
Log:
branch 2.4: Geography Distance inconsistent with Intersects
References #4480
Modified: branches/2.4/NEWS
===================================================================
--- branches/2.4/NEWS 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/NEWS 2019-08-13 16:58:58 UTC (rev 17704)
@@ -17,6 +17,7 @@
- #4461, ST_AsTWKB doesn't always remove duplicate points (Nicklas Avén)
- #4470, ST_GeomFromGeoJSON crash on empty rings (Darafei Praliaskouski)
- #4420, update path does not exists for address_standardizer extension (Regina Obe)
+ - #4480, Geography Distance inconsistent with Intersects (Paul Ramsey)
PostGIS 2.4.7
Modified: branches/2.4/liblwgeom/cunit/cu_geodetic.c
===================================================================
--- branches/2.4/liblwgeom/cunit/cu_geodetic.c 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/cunit/cu_geodetic.c 2019-08-13 16:58:58 UTC (rev 17704)
@@ -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: branches/2.4/liblwgeom/g_box.c
===================================================================
--- branches/2.4/liblwgeom/g_box.c 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/g_box.c 2019-08-13 16:58:58 UTC (rev 17704)
@@ -259,7 +259,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: branches/2.4/liblwgeom/liblwgeom.h.in
===================================================================
--- branches/2.4/liblwgeom/liblwgeom.h.in 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/liblwgeom.h.in 2019-08-13 16:58:58 UTC (rev 17704)
@@ -1837,7 +1837,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: branches/2.4/liblwgeom/lwgeodetic.c
===================================================================
--- branches/2.4/liblwgeom/lwgeodetic.c 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/lwgeodetic.c 2019-08-13 16:58:58 UTC (rev 17704)
@@ -319,34 +319,49 @@
LWDEBUG(4, "checking poles");
LWDEBUGF(4, "gbox %s", gbox_to_string(gbox));
/* 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;
@@ -353,19 +368,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;
}
@@ -458,7 +480,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;
@@ -1446,21 +1468,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;
}
/**
@@ -1467,7 +1545,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;
@@ -1534,7 +1612,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;
}
}
@@ -1543,8 +1621,8 @@
}
/* This should never happen! */
- lwerror("BOOM! Could not generate outside point!");
- return;
+ // lwerror("BOOM! Could not generate outside point!");
+ return LW_FAILURE;
}
@@ -2462,7 +2540,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);
@@ -2791,6 +2869,7 @@
return LW_SUCCESS;
}
+
int ptarray_calculate_gbox_geodetic(const POINTARRAY *pa, GBOX *gbox)
{
int i;
Modified: branches/2.4/liblwgeom/lwgeodetic.h
===================================================================
--- branches/2.4/liblwgeom/lwgeodetic.h 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/lwgeodetic.h 2019-08-13 16:58:58 UTC (rev 17704)
@@ -125,7 +125,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);
@@ -137,6 +137,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: branches/2.4/liblwgeom/lwgeodetic_tree.c
===================================================================
--- branches/2.4/liblwgeom/lwgeodetic_tree.c 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/lwgeodetic_tree.c 2019-08-13 16:58:58 UTC (rev 17704)
@@ -470,7 +470,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: branches/2.4/liblwgeom/lwgeodetic_tree.h
===================================================================
--- branches/2.4/liblwgeom/lwgeodetic_tree.h 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/liblwgeom/lwgeodetic_tree.h 2019-08-13 16:58:58 UTC (rev 17704)
@@ -53,6 +53,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: branches/2.4/postgis/geography_measurement_trees.c
===================================================================
--- branches/2.4/postgis/geography_measurement_trees.c 2019-08-11 06:24:50 UTC (rev 17703)
+++ branches/2.4/postgis/geography_measurement_trees.c 2019-08-13 16:58:58 UTC (rev 17704)
@@ -100,7 +100,6 @@
return (CircTreeGeomCache*)GetGeomCache(fcinfo, &CircTreeCacheMethods, g1, g2);
}
-
static int
CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const POINT4D* in_point)
{
@@ -142,7 +141,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