[postgis-tickets] r17705 - branch 2.5: Geography Distance inconsistent with Intersects

Paul Ramsey pramsey at cleverelephant.ca
Tue Aug 13 09:59:08 PDT 2019


Author: pramsey
Date: 2019-08-13 09:59:08 -0700 (Tue, 13 Aug 2019)
New Revision: 17705

Modified:
   branches/2.5/NEWS
   branches/2.5/liblwgeom/cunit/cu_geodetic.c
   branches/2.5/liblwgeom/g_box.c
   branches/2.5/liblwgeom/liblwgeom.h.in
   branches/2.5/liblwgeom/lwgeodetic.c
   branches/2.5/liblwgeom/lwgeodetic.h
   branches/2.5/liblwgeom/lwgeodetic_tree.c
   branches/2.5/liblwgeom/lwgeodetic_tree.h
   branches/2.5/postgis/geography_measurement_trees.c
Log:
branch 2.5: Geography Distance inconsistent with Intersects 
References #4480



Modified: branches/2.5/NEWS
===================================================================
--- branches/2.5/NEWS	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/NEWS	2019-08-13 16:59:08 UTC (rev 17705)
@@ -20,6 +20,7 @@
   - #4459, Fix ST_Subdivide crash on intermediate EMPTY (Darafei Praliaskouski)
   - #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.5.2

Modified: branches/2.5/liblwgeom/cunit/cu_geodetic.c
===================================================================
--- branches/2.5/liblwgeom/cunit/cu_geodetic.c	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/cunit/cu_geodetic.c	2019-08-13 16:59:08 UTC (rev 17705)
@@ -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.5/liblwgeom/g_box.c
===================================================================
--- branches/2.5/liblwgeom/g_box.c	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/g_box.c	2019-08-13 16:59:08 UTC (rev 17705)
@@ -254,7 +254,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.5/liblwgeom/liblwgeom.h.in
===================================================================
--- branches/2.5/liblwgeom/liblwgeom.h.in	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/liblwgeom.h.in	2019-08-13 16:59:08 UTC (rev 17705)
@@ -1829,7 +1829,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.5/liblwgeom/lwgeodetic.c
===================================================================
--- branches/2.5/liblwgeom/lwgeodetic.c	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/lwgeodetic.c	2019-08-13 16:59:08 UTC (rev 17705)
@@ -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: branches/2.5/liblwgeom/lwgeodetic.h
===================================================================
--- branches/2.5/liblwgeom/lwgeodetic.h	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/lwgeodetic.h	2019-08-13 16:59:08 UTC (rev 17705)
@@ -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.5/liblwgeom/lwgeodetic_tree.c
===================================================================
--- branches/2.5/liblwgeom/lwgeodetic_tree.c	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/lwgeodetic_tree.c	2019-08-13 16:59:08 UTC (rev 17705)
@@ -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), &center3d);
+	vector_scale(&center3d, -1.0);
+	cart2geog(&center3d, &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.5/liblwgeom/lwgeodetic_tree.h
===================================================================
--- branches/2.5/liblwgeom/lwgeodetic_tree.h	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/liblwgeom/lwgeodetic_tree.h	2019-08-13 16:59:08 UTC (rev 17705)
@@ -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: branches/2.5/postgis/geography_measurement_trees.c
===================================================================
--- branches/2.5/postgis/geography_measurement_trees.c	2019-08-13 16:58:58 UTC (rev 17704)
+++ branches/2.5/postgis/geography_measurement_trees.c	2019-08-13 16:59:08 UTC (rev 17705)
@@ -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