From 45dab34e1cbb0fa5c4a973330195cb3f070fa6e5 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 9 Oct 2009 19:23:05 +0000 Subject: [PATCH] Add implementation for ST_Covers(geography, geography) in point-in-polygon case. git-svn-id: http://svn.osgeo.org/postgis/trunk@4635 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/libgeom.h | 6 +++ liblwgeom/liblwgeom.h | 2 +- liblwgeom/lwgeodetic.c | 69 +++++++++++++++++++++++++++++++- liblwgeom/lwgeom.c | 4 +- liblwgeom/measures.c | 6 +-- postgis/geography.sql.in.c | 21 ++++++++++ postgis/geography_measurement.c | 54 ++++++++++++++++++++++++- postgis/lwgeom_dump.c | 6 +-- postgis/lwgeom_functions_basic.c | 4 +- 9 files changed, 158 insertions(+), 14 deletions(-) diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index 2a0bc6321..0d9e6d61e 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -382,6 +382,12 @@ extern double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox */ extern double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox); +/** +* Calculate covers predicate for two lwgeoms on the sphere. Currently +* only handles point-in-polygon. +*/ +extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2); + /** * New function to read doubles directly from the double* coordinate array * of an aligned lwgeom #POINTARRAY (built by de-serializing a #GSERIALIZED). diff --git a/liblwgeom/liblwgeom.h b/liblwgeom/liblwgeom.h index 801e8e746..b570fff7e 100644 --- a/liblwgeom/liblwgeom.h +++ b/liblwgeom/liblwgeom.h @@ -538,7 +538,7 @@ extern void lwgeom_add_bbox(LWGEOM *lwgeom); extern void lwgeom_dropSRID(LWGEOM *lwgeom); /* Determine whether a LWGEOM can contain sub-geometries or not */ -extern int lwgeom_contains_subgeoms(int type); +extern int lwgeom_is_collection(int type); /******************************************************************/ diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 3ac393b25..f00c1e67c 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1758,7 +1758,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX } /* Recurse into collections */ - if( lwgeom_contains_subgeoms(type1) ) + if( lwgeom_is_collection(type1) ) { int i; double distance = MAXFLOAT; @@ -1776,7 +1776,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX } /* Recurse into collections */ - if( lwgeom_contains_subgeoms(type2) ) + if( lwgeom_is_collection(type2) ) { int i; double distance = MAXFLOAT; @@ -1799,6 +1799,71 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX } + +int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2) +{ + int type1, type2; + + assert(lwgeom1); + assert(lwgeom2); + + type1 = TYPE_GETTYPE(lwgeom1->type); + type2 = TYPE_GETTYPE(lwgeom2->type); + + /* Currently a restricted implementation */ + if( ! ( (type1 == POLYGONTYPE || type1 == MULTIPOLYGONTYPE || type1 == COLLECTIONTYPE) && + (type2 == POINTTYPE || type2 == MULTIPOINTTYPE || type2 == COLLECTIONTYPE) ) ) + { + lwerror("lwgeom_covers_lwgeom_sphere: only POLYGON and POINT types are currently supported"); + return LW_FALSE; + } + + /* Handle the polygon/point case */ + if( type1 == POLYGONTYPE && type2 == POINTTYPE ) + { + POINT2D pt_to_test; + getPoint2d_p(((LWPOINT*)lwgeom2)->point, 0, &pt_to_test); + return lwpoly_covers_point2d((LWPOLY*)lwgeom1, gbox1, pt_to_test); + } + + /* If any of the first argument parts covers the second argument, it's true */ + if( lwgeom_is_collection( type1 ) ) + { + int i; + LWCOLLECTION *col = (LWCOLLECTION*)lwgeom1; + + for( i = 0; i < col->ngeoms; i++ ) + { + if( lwgeom_covers_lwgeom_sphere(col->geoms[i], lwgeom2, gbox1, gbox2) ) + { + return LW_TRUE; + } + } + return LW_FALSE; + } + + /* Only if all of the second arguments are covered by the first argument is the condition true */ + if( lwgeom_is_collection( type2 ) ) + { + int i; + LWCOLLECTION *col = (LWCOLLECTION*)lwgeom2; + + for( i = 0; i < col->ngeoms; i++ ) + { + if( ! lwgeom_covers_lwgeom_sphere(lwgeom1, col->geoms[i], gbox1, gbox2) ) + { + return LW_FALSE; + } + } + return LW_TRUE; + } + + /* Don't get here */ + lwerror("lwgeom_covers_lwgeom_sphere: reached end of function without resolution"); + return LW_FALSE; + +} + /** * Given a polygon (lon/lat decimal degrees) and point (lon/lat decimal degrees) and * a guaranteed outside point (lon/lat decimal degrees) (calculate with gbox_pt_outside()) diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c index f5d7aed0c..c65749268 100644 --- a/liblwgeom/lwgeom.c +++ b/liblwgeom/lwgeom.c @@ -404,7 +404,7 @@ lwgeom_as_multi(LWGEOM *lwgeom) ** This funx is a no-op only if a bbox cache is already present ** in input. */ - if ( lwgeom_contains_subgeoms(TYPE_GETTYPE(lwgeom->type)) ) + if ( lwgeom_is_collection(TYPE_GETTYPE(lwgeom->type)) ) { return lwgeom_clone(lwgeom); } @@ -928,7 +928,7 @@ lwgeom_longitude_shift(LWGEOM *lwgeom) /** Return TRUE if the geometry may contain sub-geometries, i.e. it is a MULTI* or COMPOUNDCURVE */ int -lwgeom_contains_subgeoms(int type) +lwgeom_is_collection(int type) { switch (type) diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c index ea6592b81..6f8ff375b 100644 --- a/liblwgeom/measures.c +++ b/liblwgeom/measures.c @@ -667,7 +667,7 @@ lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double toleranc double dist=tolerance; /* Argument 1 is a multitype... recurse */ - if ( lwgeom_contains_subgeoms(t1) ) + if ( lwgeom_is_collection(t1) ) { dist = lwgeom_mindistance2d_recursive_tolerance(g1, lw2, tolerance); if ( dist <= tolerance ) return tolerance; /* can't be closer */ @@ -682,7 +682,7 @@ lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double toleranc int t2 = lwgeom_getType(g2[0]); /* Argument 2 is a multitype... recurse */ - if ( lwgeom_contains_subgeoms(t2) ) + if ( lwgeom_is_collection(t2) ) { dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); if ( dist <= tolerance ) return tolerance; /* can't be closer */ @@ -775,7 +775,7 @@ lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double toleranc lwerror("Unsupported geometry type: %s", lwgeom_typename(t2)); } } -// else if (lwgeom_contains_subgeoms(t1)) /* it's a multitype... recurse */ +// else if (lwgeom_is_collection(t1)) /* it's a multitype... recurse */ // { // dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance); // } diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 720790943..e40f4601f 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -462,6 +462,27 @@ CREATE OR REPLACE FUNCTION ST_PointOutside(geography) AS 'MODULE_PATHNAME','geography_point_outside' LANGUAGE 'C' IMMUTABLE STRICT; +-- Only implemented for polygon-over-point +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION _ST_Covers(geography, geography) + RETURNS boolean + AS 'MODULE_PATHNAME','geography_covers' + LANGUAGE 'C' IMMUTABLE STRICT; + +-- Only implemented for polygon-over-point +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_Covers(geography, geography) + RETURNS boolean + AS 'SELECT $1 && $2 AND _ST_Covers($1, $2)' + LANGUAGE 'SQL' IMMUTABLE STRICT; + +-- Only implemented for polygon-over-point +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_CoveredBy(geography, geography) + RETURNS boolean + AS 'SELECT $1 && $2 AND _ST_Covers($2, $1)' + LANGUAGE 'SQL' IMMUTABLE STRICT; + -- ---------- ---------- ---------- ---------- ---------- ---------- ---------- COMMIT; diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index 45efdd5f1..ceb691e0e 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -26,6 +26,8 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS); Datum geography_area_sphere(PG_FUNCTION_ARGS); Datum geography_expand(PG_FUNCTION_ARGS); +Datum geography_point_outside(PG_FUNCTION_ARGS); +Datum geography_covers(PG_FUNCTION_ARGS); /* @@ -210,4 +212,54 @@ Datum geography_point_outside(PG_FUNCTION_ARGS) PG_RETURN_POINTER(g_out); -} \ No newline at end of file +} + +PG_FUNCTION_INFO_V1(geography_covers); +Datum geography_covers(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom1 = NULL; + LWGEOM *lwgeom2 = NULL; + GBOX gbox1; + GBOX gbox2; + GSERIALIZED *g1 = NULL; + GSERIALIZED *g2 = NULL; + int type1, type2; + int result = LW_FALSE; + + /* Get our geometry objects loaded into memory. */ + g1 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + g2 = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); + + 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(); + } + + /* We need the bounding boxes in case of polygon calculations, + which requires them to generate a stab-line to test point-in-polygon. */ + if( ! gbox_from_gserialized(g1, &gbox1) || + ! gbox_from_gserialized(g2, &gbox2) ) + { + elog(ERROR, "geography_covers: error in gbox_from_gserialized calculation."); + PG_RETURN_NULL(); + } + + /* Construct our working geometries */ + lwgeom1 = lwgeom_from_gserialized(g1); + lwgeom2 = lwgeom_from_gserialized(g2); + + /* Calculate answer */ + result = lwgeom_covers_lwgeom_sphere(lwgeom1, lwgeom2, gbox1, gbox2); + + /* Clean up, but not all the way to the point arrays */ + lwgeom_release(lwgeom1); + lwgeom_release(lwgeom2); + + PG_RETURN_BOOL(result); +} diff --git a/postgis/lwgeom_dump.c b/postgis/lwgeom_dump.c index e439c6fb3..af9717427 100644 --- a/postgis/lwgeom_dump.c +++ b/postgis/lwgeom_dump.c @@ -89,7 +89,7 @@ Datum LWGEOM_dump(PG_FUNCTION_ARGS) state->root = lwgeom; state->stacklen=0; - if ( lwgeom_contains_subgeoms(TYPE_GETTYPE(lwgeom->type)) ) + if ( lwgeom_is_collection(TYPE_GETTYPE(lwgeom->type)) ) { /* * Push a GEOMDUMPNODE on the state stack @@ -127,7 +127,7 @@ Datum LWGEOM_dump(PG_FUNCTION_ARGS) /* Handled simple geometries */ if ( ! state->root ) SRF_RETURN_DONE(funcctx); - if ( ! lwgeom_contains_subgeoms(TYPE_GETTYPE(state->root->type)) ) + if ( ! lwgeom_is_collection(TYPE_GETTYPE(state->root->type)) ) { values[0] = "{}"; values[1] = lwgeom_to_hexwkb(state->root, PARSER_CHECK_NONE, -1); @@ -146,7 +146,7 @@ Datum LWGEOM_dump(PG_FUNCTION_ARGS) if ( node->idx < lwcoll->ngeoms ) { lwgeom = lwcoll->geoms[node->idx]; - if ( ! lwgeom_contains_subgeoms(TYPE_GETTYPE(lwgeom->type)) ) + if ( ! lwgeom_is_collection(TYPE_GETTYPE(lwgeom->type)) ) { /* write address of current geom */ ptr=address; diff --git a/postgis/lwgeom_functions_basic.c b/postgis/lwgeom_functions_basic.c index 129790fd9..bebcdcc04 100644 --- a/postgis/lwgeom_functions_basic.c +++ b/postgis/lwgeom_functions_basic.c @@ -1460,7 +1460,7 @@ Datum LWGEOM_force_collection(PG_FUNCTION_ARGS) lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom)); /* alread a multi*, just make it a collection */ - if ( lwgeom_contains_subgeoms(TYPE_GETTYPE(lwgeom->type)) ) + if ( lwgeom_is_collection(TYPE_GETTYPE(lwgeom->type)) ) { TYPE_SETTYPE(lwgeom->type, COLLECTIONTYPE); } @@ -1503,7 +1503,7 @@ Datum LWGEOM_force_multi(PG_FUNCTION_ARGS) ** in input. If bbox cache is not there we'll need to handle ** automatic bbox addition FOR_COMPLEX_GEOMS. */ - if ( lwgeom_contains_subgeoms(TYPE_GETTYPE(geom->type)) && TYPE_HASBBOX(geom->type) ) + if ( lwgeom_is_collection(TYPE_GETTYPE(geom->type)) && TYPE_HASBBOX(geom->type) ) { PG_RETURN_POINTER(geom); } -- 2.50.1