[postgis-tickets] r15599 - Covers support for Polygon-on-polygon, line on line, point on line, point on point for geography
Regina Obe
lr at pcorp.us
Sun Aug 27 18:53:32 PDT 2017
Author: robe
Date: 2017-08-27 18:53:32 -0700 (Sun, 27 Aug 2017)
New Revision: 15599
Added:
trunk/regress/geography_covers.sql
trunk/regress/geography_covers_expected
Modified:
trunk/NEWS
trunk/doc/reference_measure.xml
trunk/liblwgeom/lwgeodetic.c
trunk/liblwgeom/lwgeodetic.h
trunk/postgis/geography_measurement.c
Log:
Covers support for Polygon-on-polygon, line on line, point on line, point on point for geography
Patch from Danny G?\195?\182tte
Closes #524
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2017-08-27 15:17:39 UTC (rev 15598)
+++ trunk/NEWS 2017-08-28 01:53:32 UTC (rev 15599)
@@ -19,6 +19,8 @@
and all stable / immutable (raster and geometry) marked as parallel safe
- #2249, ST_MakeEmptyCoverage for raster (David Zwarg, ainomieli)
- #3709, Allow signed distance for ST_Project (Darafei Praliaskouski)
+ - #524, Covers support for Polygon-on-polygon, line on line, point on line,
+ point on point for geography (Danny Götte)
* Enhancements *
Modified: trunk/doc/reference_measure.xml
===================================================================
--- trunk/doc/reference_measure.xml 2017-08-27 15:17:39 UTC (rev 15598)
+++ trunk/doc/reference_measure.xml 2017-08-28 01:53:32 UTC (rev 15599)
@@ -1754,10 +1754,6 @@
</important>
<important>
- <para>For geography only Polygon covers point is supported.</para>
- </important>
-
- <important>
<para>Do not use this function with invalid geometries. You will get unexpected results.</para>
</important>
@@ -1766,6 +1762,7 @@
the geometries. To avoid index use, use the function
_ST_Covers.</para>
+ <para>Enhanced: 2.4.0 Support for polygon in polygon and line in polygon added for geography type</para>
<para>Enhanced: 2.3.0 Enhancement to PIP short-circuit for geometry extended to support MultiPoints with few points. Prior versions only supported point in polygon.</para>
<para>Availability: 1.5 - support for geography was introduced. </para>
<para>Availability: 1.2.2 - requires GEOS >= 3.0</para>
Modified: trunk/liblwgeom/lwgeodetic.c
===================================================================
--- trunk/liblwgeom/lwgeodetic.c 2017-08-27 15:17:39 UTC (rev 15598)
+++ trunk/liblwgeom/lwgeodetic.c 2017-08-28 01:53:32 UTC (rev 15599)
@@ -2317,18 +2317,19 @@
int type1, type2;
GBOX gbox1, gbox2;
gbox1.flags = gbox2.flags = 0;
-
+
assert(lwgeom1);
assert(lwgeom2);
type1 = lwgeom1->type;
type2 = lwgeom2->type;
- /* Currently a restricted implementation */
- if ( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) &&
- (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) )
+ /* dim(geom2) > dim(geom1) always returns false (because geom2 is bigger) */
+ if ( (type1 == POINTTYPE && type2 == LINETYPE)
+ || (type1 == POINTTYPE && type2 == POLYGONTYPE)
+ || (type1 == LINETYPE && type2 == POLYGONTYPE) )
{
- lwerror("lwgeom_covers_lwgeom_sphere: only POLYGON covers POINT tests are currently supported");
+ LWDEBUG(4, "dimension of geom2 is bigger than geom1");
return LW_FALSE;
}
@@ -2352,6 +2353,26 @@
getPoint2d_p(((LWPOINT*)lwgeom2)->point, 0, &pt_to_test);
return lwpoly_covers_point2d((LWPOLY*)lwgeom1, &pt_to_test);
}
+ else if ( type1 == POLYGONTYPE && type2 == LINETYPE)
+ {
+ return lwpoly_covers_lwline((LWPOLY*)lwgeom1, (LWLINE*)lwgeom2);
+ }
+ else if ( type1 == POLYGONTYPE && type2 == POLYGONTYPE)
+ {
+ return lwpoly_covers_lwpoly((LWPOLY*)lwgeom1, (LWPOLY*)lwgeom2);
+ }
+ else if ( type1 == LINETYPE && type2 == POINTTYPE)
+ {
+ return lwline_covers_lwpoint((LWLINE*)lwgeom1, (LWPOINT*)lwgeom2);
+ }
+ else if ( type1 == LINETYPE && type2 == LINETYPE)
+ {
+ return lwline_covers_lwline((LWLINE*)lwgeom1, (LWLINE*)lwgeom2);
+ }
+ else if ( type1 == POINTTYPE && type2 == POINTTYPE)
+ {
+ return lwpoint_same((LWPOINT*)lwgeom1, (LWPOINT*)lwgeom2);
+ }
/* If any of the first argument parts covers the second argument, it's true */
if ( lwtype_is_collection( type1 ) )
@@ -2466,8 +2487,281 @@
return LW_TRUE;
}
+/**
+ * Given a polygon1 check if all points of polygon2 are inside polygon1 and no
+ * intersections of the polygon edges occur.
+ * return LW_TRUE if polygon is inside or on edge of polygon.
+ */
+int lwpoly_covers_lwpoly(const LWPOLY *poly1, const LWPOLY *poly2)
+{
+ int i;
+ /* Nulls and empties don't contain anything! */
+ if ( ! poly1 || lwgeom_is_empty((LWGEOM*)poly1) )
+ {
+ LWDEBUG(4,"returning false, geometry1 is empty or null");
+ return LW_FALSE;
+ }
+
+ /* Nulls and empties don't contain anything! */
+ if ( ! poly2 || lwgeom_is_empty((LWGEOM*)poly2) )
+ {
+ LWDEBUG(4,"returning false, geometry2 is empty or null");
+ return LW_FALSE;
+ }
+
+ /* check if all vertices of poly2 are inside poly1 */
+ for (i = 0; i < poly2->nrings; i++)
+ {
+
+ /* every other ring is a hole, check if point is inside the actual polygon */
+ if ( i % 2 == 0)
+ {
+ if (LW_FALSE == lwpoly_covers_pointarray(poly1, poly2->rings[i]))
+ {
+ LWDEBUG(4,"returning false, geometry2 has point outside of geometry1");
+ return LW_FALSE;
+ }
+ }
+ else
+ {
+ if (LW_TRUE == lwpoly_covers_pointarray(poly1, poly2->rings[i]))
+ {
+ LWDEBUG(4,"returning false, geometry2 has point inside a hole of geometry1");
+ return LW_FALSE;
+ }
+ }
+ }
+
+ /* check for any edge intersections, so nothing is partially outside of poly1 */
+ for (i = 0; i < poly2->nrings; i++)
+ {
+ if (LW_TRUE == lwpoly_intersects_line(poly1, poly2->rings[i]))
+ {
+ LWDEBUG(4,"returning false, geometry2 is partially outside of geometry1");
+ return LW_FALSE;
+ }
+ }
+
+ /* no abort condition found, so the poly2 should be completly inside poly1 */
+ return LW_TRUE;
+}
+
/**
+ *
+ */
+int lwpoly_covers_lwline(const LWPOLY *poly, const LWLINE *line)
+{
+ /* Nulls and empties don't contain anything! */
+ if ( ! poly || lwgeom_is_empty((LWGEOM*)poly) )
+ {
+ LWDEBUG(4,"returning false, geometry1 is empty or null");
+ return LW_FALSE;
+ }
+
+ /* Nulls and empties don't contain anything! */
+ if ( ! line || lwgeom_is_empty((LWGEOM*)line) )
+ {
+ LWDEBUG(4,"returning false, geometry2 is empty or null");
+ return LW_FALSE;
+ }
+
+ if (LW_FALSE == lwpoly_covers_pointarray(poly, line->points))
+ {
+ LWDEBUG(4,"returning false, geometry2 has point outside of geometry1");
+ return LW_FALSE;
+ }
+
+ /* check for any edge intersections, so nothing is partially outside of poly1 */
+ if (LW_TRUE == lwpoly_intersects_line(poly, line->points))
+ {
+ LWDEBUG(4,"returning false, geometry2 is partially outside of geometry1");
+ return LW_FALSE;
+ }
+
+ /* no abort condition found, so the poly2 should be completly inside poly1 */
+ return LW_TRUE;
+}
+
+/**
+ * return LW_TRUE if all points are inside the polygon
+ */
+int lwpoly_covers_pointarray(const LWPOLY* lwpoly, const POINTARRAY* pta)
+{
+ int i;
+ for (i = 0; i < pta->npoints; i++) {
+ const POINT2D* pt_to_test = getPoint2d_cp(pta, i);
+
+ if ( LW_FALSE == lwpoly_covers_point2d(lwpoly, pt_to_test) ) {
+ LWDEBUG(4,"returning false, geometry2 has point outside of geometry1");
+ return LW_FALSE;
+ }
+ }
+
+ return LW_TRUE;
+}
+
+/**
+ * Checks if any edges of lwpoly intersect with the line formed by the pointarray
+ * return LW_TRUE if any intersection beetween the given polygon and the line
+ */
+int lwpoly_intersects_line(const LWPOLY* lwpoly, const POINTARRAY* line)
+{
+ int i, j, k;
+ POINT3D pa1, pa2, pb1, pb2;
+ for (i = 0; i < lwpoly->nrings; i++)
+ {
+ for (j = 0; j < lwpoly->rings[i]->npoints - 1; j++)
+ {
+ const POINT2D* a1 = getPoint2d_cp(lwpoly->rings[i], j);
+ const POINT2D* a2 = getPoint2d_cp(lwpoly->rings[i], j+1);
+
+ /* Set up our stab line */
+ ll2cart(a1, &pa1);
+ ll2cart(a2, &pa2);
+
+ for (k = 0; k < line->npoints - 1; k++)
+ {
+ const POINT2D* b1 = getPoint2d_cp(line, k);
+ const POINT2D* b2 = getPoint2d_cp(line, k+1);
+
+ /* Set up our stab line */
+ ll2cart(b1, &pb1);
+ ll2cart(b2, &pb2);
+
+ int inter = edge_intersects(&pa1, &pa2, &pb1, &pb2);
+
+ /* ignore same edges */
+ if (inter & PIR_INTERSECTS
+ && !(inter & PIR_B_TOUCH_RIGHT || inter & PIR_COLINEAR) )
+ {
+ return LW_TRUE;
+ }
+ }
+ }
+ }
+
+ return LW_FALSE;
+}
+
+/**
+ * return LW_TRUE if any of the line segments covers the point
+ */
+int lwline_covers_lwpoint(const LWLINE* lwline, const LWPOINT* lwpoint)
+{
+ int i;
+ GEOGRAPHIC_POINT p;
+ GEOGRAPHIC_EDGE e;
+
+ for ( i = 0; i < lwline->points->npoints - 1; i++)
+ {
+ const POINT2D* a1 = getPoint2d_cp(lwline->points, i);
+ const POINT2D* a2 = getPoint2d_cp(lwline->points, i+1);
+
+ geographic_point_init(a1->x, a1->y, &(e.start));
+ geographic_point_init(a2->x, a2->y, &(e.end));
+
+ geographic_point_init(lwpoint_get_x(lwpoint), lwpoint_get_y(lwpoint), &p);
+
+ if ( edge_contains_point(&e, &p) ) {
+ return LW_TRUE;
+ }
+ }
+
+ return LW_FALSE;
+}
+
+/**
+ * Check if first and last point of line2 are covered by line1 and then each
+ * point in between has to be one line1 in the exact same order
+ * return LW_TRUE if all edge points of line2 are on line1
+ */
+int lwline_covers_lwline(const LWLINE* lwline1, const LWLINE* lwline2)
+{
+ int i, j;
+ GEOGRAPHIC_EDGE e1, e2;
+ GEOGRAPHIC_POINT p1, p2;
+ int start = LW_FALSE;
+ int changed = LW_FALSE;
+
+ /* first point on line */
+ if ( ! lwline_covers_lwpoint(lwline1, lwline_get_lwpoint(lwline2, 0)))
+ {
+ LWDEBUG(4,"returning false, first point of line2 is not covered by line1");
+ return LW_FALSE;
+ }
+
+ /* last point on line */
+ if ( ! lwline_covers_lwpoint(lwline1, lwline_get_lwpoint(lwline2, lwline2->points->npoints - 1)))
+ {
+ LWDEBUG(4,"returning false, last point of line2 is not covered by line1");
+ return LW_FALSE;
+ }
+
+ j = 0;
+ i = 0;
+ while (i < lwline1->points->npoints - 1 && j < lwline2->points->npoints - 1)
+ {
+ changed = LW_FALSE;
+ const POINT2D* a1 = getPoint2d_cp(lwline1->points, i);
+ const POINT2D* a2 = getPoint2d_cp(lwline1->points, i+1);
+ const POINT2D* b1 = getPoint2d_cp(lwline2->points, j);
+ const POINT2D* b2 = getPoint2d_cp(lwline2->points, j+1);
+
+ geographic_point_init(a1->x, a1->y, &(e1.start));
+ geographic_point_init(a2->x, a2->y, &(e1.end));
+ geographic_point_init(b1->x, b1->y, &p2);
+
+ /* we already know, that the last point is on line1, so we're done */
+ if ( j == lwline2->points->npoints - 1)
+ {
+ return LW_TRUE;
+ }
+ else if (start == LW_TRUE)
+ {
+ /* point is on current line1 edge, check next point in line2 */
+ if ( edge_contains_point(&e1, &p2)) {
+ j++;
+ changed = LW_TRUE;
+ }
+
+ geographic_point_init(a1->x, a1->y, &(e2.start));
+ geographic_point_init(a2->x, b2->y, &(e2.end));
+ geographic_point_init(a1->x, a1->y, &p1);
+
+ /* point is on current line2 edge, check next point in line1 */
+ if ( edge_contains_point(&e2, &p1)) {
+ i++;
+ changed = LW_TRUE;
+ }
+
+ /* no edge progressed -> point left one line */
+ if ( changed == LW_FALSE )
+ {
+ LWDEBUG(4,"returning false, found point not covered by both lines");
+ return LW_FALSE;
+ }
+ else
+ {
+ continue;
+ }
+ }
+
+ /* find first edge to cover line2 */
+ if (edge_contains_point(&e1, &p2))
+ {
+ start = LW_TRUE;
+ }
+
+ /* next line1 edge */
+ i++;
+ }
+
+ /* no uncovered point found */
+ return LW_TRUE;
+}
+
+/**
* This function can only be used on LWGEOM that is built on top of
* GSERIALIZED, otherwise alignment errors will ensue.
*/
Modified: trunk/liblwgeom/lwgeodetic.h
===================================================================
--- trunk/liblwgeom/lwgeodetic.h 2017-08-27 15:17:39 UTC (rev 15598)
+++ trunk/liblwgeom/lwgeodetic.h 2017-08-28 01:53:32 UTC (rev 15599)
@@ -119,6 +119,12 @@
void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g);
int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test);
int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test);
+int lwpoly_covers_lwpoly(const LWPOLY *lwpoly1, const LWPOLY *lwpoly2);
+int lwpoly_covers_pointarray(const LWPOLY* lwpoly, const POINTARRAY* pta);
+int lwpoly_covers_lwline(const LWPOLY *poly, const LWLINE *line);
+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 ptarray_point_in_ring(const POINTARRAY *pa, const POINT2D *pt_outside, const POINT2D *pt_to_test);
double ptarray_area_sphere(const POINTARRAY *pa);
Modified: trunk/postgis/geography_measurement.c
===================================================================
--- trunk/postgis/geography_measurement.c 2017-08-27 15:17:39 UTC (rev 15598)
+++ trunk/postgis/geography_measurement.c 2017-08-28 01:53:32 UTC (rev 15599)
@@ -742,14 +742,6 @@
type1 = gserialized_get_type(g1);
type2 = gserialized_get_type(g2);
- /* Right now we only handle points and polygons */
- if ( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) &&
- (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) )
- {
- elog(ERROR, "geography_covers: only POLYGON and POINT types are currently supported");
- PG_RETURN_NULL();
- }
-
/* Construct our working geometries */
lwgeom1 = lwgeom_from_gserialized(g1);
lwgeom2 = lwgeom_from_gserialized(g2);
Added: trunk/regress/geography_covers.sql
===================================================================
--- trunk/regress/geography_covers.sql (rev 0)
+++ trunk/regress/geography_covers.sql 2017-08-28 01:53:32 UTC (rev 15599)
@@ -0,0 +1,45 @@
+
+-- poly in poly
+SELECT c, ST_Covers(g1::geography, g2::geography) FROM
+( VALUES
+ ('geog_covers_poly_poly_in', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'POLYGON (( 10 30, 30 30, 30 10, 10 10, 10 30))'),
+ ('geog_covers_poly_poly_part', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'POLYGON (( -10 30, 30 30, 30 10, 10 10, -10 30))'),
+ ('geog_covers_poly_poly_out', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'POLYGON ((0 -40, -40 -40, -40 0, 0 0, 0 -40))'),
+ ('geog_covers_poly_poly_idl', 'POLYGON((-170 -20, 170 -20, 170 20, -170 20, -170 -20))', 'POLYGON((-175 -20, 175 -20, 175 20, -175 20, -175 -20))'),
+ ('geog_covers_poly_poly_pole', 'POLYGON((0 80, 90 80, 180 80, -90 80, 0 80))', 'POLYGON((0 85, 90 85, 180 85, -90 85, 0 85))')
+) AS u(c, g1, g2);
+
+-- line in poly
+SELECT c, ST_Covers(g1::geography, g2::geography) FROM
+( VALUES
+ ('geog_covers_poly_line_in', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'LINESTRING (35 35, 35 15, 15 15)'),
+ ('geog_covers_poly_line_part', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'LINESTRING (-10 30, 30 30, 30 10, 10 10)'),
+ ('geog_covers_poly_line_out', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'LINESTRING (-10 -40, -40 -40, -40 -10, -10 -10)')
+) AS u(c, g1, g2);
+
+
+-- poly in poly (reversed arguments)
+SELECT c, ST_CoveredBy(g1::geography, g2::geography) FROM
+( VALUES
+ ('geog_covered_poly_poly_in', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'POLYGON (( 10 30, 30 30, 30 10, 10 10, 10 30))'),
+ ('geog_covered_poly_poly_part', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'POLYGON (( -10 30, 30 30, 30 10, 10 10, -10 30))'),
+ ('geog_covered_poly_poly_out', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'POLYGON ((0 -40, -40 -40, -40 0, 0 0, 0 -40))'),
+ ('geog_covered_poly_poly_idl', 'POLYGON((-170 -20, 170 -20, 170 20, -170 20, -170 -20))', 'POLYGON((-175 -20, 175 -20, 175 20, -175 20, -175 -20))'),
+ ('geog_covered_poly_poly_pole', 'POLYGON((0 80, 90 80, 180 80, -90 80, 0 80))', 'POLYGON((0 85, 90 85, 180 85, -90 85, 0 85))')
+) AS u(c, g1, g2);
+
+-- line in poly (reversed arguments)
+SELECT c, ST_CoveredBy(g1::geography, g2::geography) FROM
+( VALUES
+ ('geog_covered_poly_line_in', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'LINESTRING (35 35, 35 15, 15 15)'),
+ ('geog_covered_poly_line_part', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'LINESTRING (-10 30, 30 30, 30 10, 10 10)'),
+ ('geog_covered_poly_line_out', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))', 'LINESTRING (-10 -40, -40 -40, -40 -10, -10 -10)')
+) AS u(c, g1, g2);
+
+-- self coverage
+SELECT c, ST_Covers(g::geography, g::geography) FROM
+( VALUES
+ ('geog_covers_self_point', 'POINT (0 20)'),
+ ('geog_covers_self_line', 'LINESTRING (35 35, 35 15, 15 15)'),
+ ('geog_covers_self_polygon', 'POLYGON((0 40, 40 40, 40 0, 0 0, 0 40))')
+) AS u(c, g)
Added: trunk/regress/geography_covers_expected
===================================================================
--- trunk/regress/geography_covers_expected (rev 0)
+++ trunk/regress/geography_covers_expected 2017-08-28 01:53:32 UTC (rev 15599)
@@ -0,0 +1,19 @@
+geog_covers_poly_poly_in|t
+geog_covers_poly_poly_part|f
+geog_covers_poly_poly_out|f
+geog_covers_poly_poly_idl|t
+geog_covers_poly_poly_pole|t
+geog_covers_poly_line_in|t
+geog_covers_poly_line_part|f
+geog_covers_poly_line_out|f
+geog_covered_poly_poly_in|f
+geog_covered_poly_poly_part|f
+geog_covered_poly_poly_out|f
+geog_covered_poly_poly_idl|f
+geog_covered_poly_poly_pole|f
+geog_covered_poly_line_in|f
+geog_covered_poly_line_part|f
+geog_covered_poly_line_out|f
+geog_covers_self_point|t
+geog_covers_self_line|t
+geog_covers_self_polygon|t
More information about the postgis-tickets
mailing list