]> granicus.if.org Git - postgis/commitdiff
Covers support for Polygon-on-polygon, line on line, point on line, point on point...
authorRegina Obe <lr@pcorp.us>
Mon, 28 Aug 2017 01:53:32 +0000 (01:53 +0000)
committerRegina Obe <lr@pcorp.us>
Mon, 28 Aug 2017 01:53:32 +0000 (01:53 +0000)
Patch from Danny Götte
Closes #524

git-svn-id: http://svn.osgeo.org/postgis/trunk@15599 b70326c6-7e19-0410-871a-916f4a2858ee

NEWS
doc/reference_measure.xml
liblwgeom/lwgeodetic.c
liblwgeom/lwgeodetic.h
postgis/geography_measurement.c
regress/geography_covers.sql [new file with mode: 0644]
regress/geography_covers_expected [new file with mode: 0644]

diff --git a/NEWS b/NEWS
index 5f36c701a23032f70981ca67aca8b132d1c77788..83c8f0ad1e599dfc78eee01d10bea127f44fcb5a 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -19,6 +19,8 @@ PostGIS 2.4.0
     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 *
 
index dcde1a9e08baeee5c29d3ad860fa7c892a0fd88a..0a4cb0ebe5b0f27114b8f559564556c774034cc9 100644 (file)
@@ -1753,10 +1753,6 @@ ST_Point      | t          | t              | f           | f
                  <para>Do not call with a <varname>GEOMETRYCOLLECTION</varname> as an argument</para>
                </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 @@ ST_Point      | t          | t              | f           | f
                        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 &gt;= 3.0</para>
index e96b9584074ea5803c79c389e2da17c78e8ee558..54234be0d1016b1ad616b328cdf32abbd187d649 100644 (file)
@@ -2317,18 +2317,19 @@ int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
        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 @@ int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
                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,6 +2487,279 @@ int lwpoly_covers_point2d(const LWPOLY *poly, const POINT2D *pt_to_test)
        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
index 40c2ced5e1d5d9b49497e1dc5712ae9f8a925fdc..a43b3f1b722e5209ac56ce31b80e7f8c286db92c 100644 (file)
@@ -119,6 +119,12 @@ double edge_distance_to_edge(const GEOGRAPHIC_EDGE *e1, const GEOGRAPHIC_EDGE *e
 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);
index be76d0b256b4be33adcd59bc5571ad762ab3d93e..81529464dfa15ac856f38e8c04254c07910e7a4b 100644 (file)
@@ -742,14 +742,6 @@ Datum geography_covers(PG_FUNCTION_ARGS)
        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);
diff --git a/regress/geography_covers.sql b/regress/geography_covers.sql
new file mode 100644 (file)
index 0000000..021eb51
--- /dev/null
@@ -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)
diff --git a/regress/geography_covers_expected b/regress/geography_covers_expected
new file mode 100644 (file)
index 0000000..792ccc5
--- /dev/null
@@ -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