Add ST_PointOutside() function for testing purposes.

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

diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h
index f234e589b2909d47bc84c2ceec9da2ca6c3e768b..2a0bc632156cb8fa21440c86911c7468880750d5 100644 (file)
--- a/liblwgeom/libgeom.h
+++ b/liblwgeom/libgeom.h
@@ -393,6 +393,11 @@ extern int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point);
 */
 extern int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox);
 
+/**
+* Calculate a spherical point that falls outside the geocentric gbox
+*/
+void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside);
+
 /**
 * Create a new gbox with the dimensionality indicated by the flags. Caller
 * is responsible for freeing.

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index a0f60575b0a3c420f435ae392f651ce5c2e9a99e..3ac393b255ae338e20735da59bdd06e1fbbdf7c1 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1198,16 +1198,18 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
 }
 
 /**
-* Given a gbox, return a cartesian unit vector to a point that is 
+* 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).
 */
-static void gbox_pt_outside(GBOX gbox, POINT3D *pt)
+void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside)
 {
        static double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
        double d;
        int i;
        GBOX ge;
        POINT3D corners[8];
+       POINT3D pt;
+       GEOGRAPHIC_POINT g;
        
        /* Assign our box and expand it slightly. */
        ge = gbox;
@@ -1256,40 +1258,49 @@ static void gbox_pt_outside(GBOX gbox, POINT3D *pt)
                normalize(&(corners[i]));
                if( ! gbox_contains_point3d(gbox, corners[i]) )
                {
-                       *pt = corners[i];
+                       pt = corners[i];
+                       normalize(&pt);
+                       cart2geog(pt, &g);
+                       pt_outside->x = rad2deg(g.lon);
+                       pt_outside->y = rad2deg(g.lat);
                        return;
                }       
        }
 
-       pt->x = 1.0;
-       pt->y = 0.0;
-       pt->z = 0.0;
+       pt.x = 1.0;
+       pt.y = 0.0;
+       pt.z = 0.0;
        
        if((1.0 - gbox.xmax) > 0.1) 
        {
-               pt->x = gbox.xmax + (1.0 - gbox.xmax) * 0.01;
-               d = sqrt((1.0 - pow(pt->x, 2.0))/2.0);
-               pt->y = d;
-               pt->z = d;
+               pt.x = gbox.xmax + (1.0 - gbox.xmax) * 0.01;
+               d = sqrt((1.0 - pow(pt.x, 2.0))/2.0);
+               pt.y = d;
+               pt.z = d;
        } 
        else if((1.0 - gbox.ymax) > 0.1) 
        {
-               pt->y = gbox.ymax + (1.0 - gbox.ymax) * 0.01;
-               d = sqrt((1.0 - pow(pt->y, 2.0))/2.0);
-               pt->x = d;
-               pt->z = d;
+               pt.y = gbox.ymax + (1.0 - gbox.ymax) * 0.01;
+               d = sqrt((1.0 - pow(pt.y, 2.0))/2.0);
+               pt.x = d;
+               pt.z = d;
        } 
        else if((1.0 - gbox.zmax) > 0.1) 
        {
-               pt->z = gbox.zmax + (1.0 - gbox.zmax) * 0.01;
-               d = sqrt((1.0 - pow(pt->z, 2.0))/2.0);
-               pt->x = d;
-               pt->y = d;
+               pt.z = gbox.zmax + (1.0 - gbox.zmax) * 0.01;
+               d = sqrt((1.0 - pow(pt.z, 2.0))/2.0);
+               pt.x = d;
+               pt.y = d;
        }
-       normalize(pt);
+       normalize(&pt);
+       cart2geog(pt, &g);
+       pt_outside->x = rad2deg(g.lon);
+       pt_outside->y = rad2deg(g.lat);
        return;
 }
 
+
+
 /**
 * Returns the area of the ring (ring must be closed) in square radians (surface of
 * the sphere is 4*PI).
@@ -1514,8 +1525,6 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t
 double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox)
 {
        int type;
-       POINT3D p;
-       GEOGRAPHIC_POINT g;
        POINT2D pt_outside;
        
        assert(lwgeom);
@@ -1531,10 +1540,9 @@ double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox)
        if( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
                return 0.0;
 
-       gbox_pt_outside(gbox, &p);
-       cart2geog(p, &g);
-       pt_outside.x = rad2deg(g.lon);
-       pt_outside.y = rad2deg(g.lat);
+       gbox_pt_outside(gbox, &pt_outside);
+       
+       LWDEBUGF(2, "pt_outside = POINT(%.20g %.20g)", pt_outside.x, pt_outside.y);
        
        /* Actually calculate area */
        if( type == POLYGONTYPE )
@@ -1800,8 +1808,8 @@ int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test)
 {
        int i;
        int in_hole_count = 0;
-       POINT3D p, q;
-       GEOGRAPHIC_POINT g, gpt_to_test;
+       POINT3D p;
+       GEOGRAPHIC_POINT gpt_to_test;
        POINT2D pt_outside;
        
        /* Nulls and empties don't contain anything! */
@@ -1813,15 +1821,12 @@ int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test)
 
        /* Point not in box? Done! */
        geographic_point_init(pt_to_test.x, pt_to_test.y, &gpt_to_test);
-       geog2cart(gpt_to_test, &q);
-       if( ! gbox_contains_point3d(gbox, q) )
+       geog2cart(gpt_to_test, &p);
+       if( ! gbox_contains_point3d(gbox, p) )
                return LW_FALSE;
        
        /* Calculate our outside point from the gbox */
-       gbox_pt_outside(gbox, &p);
-       cart2geog(p, &g);
-       pt_outside.x = rad2deg(g.lon);
-       pt_outside.y = rad2deg(g.lat);
+       gbox_pt_outside(gbox, &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);

diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c
index eb2eb1ac9cd0a0676ad8f82dbe8b140f098fe1d2..7207909439ec8695a67eec98951fe37d62aa6d21 100644 (file)
--- a/postgis/geography.sql.in.c
+++ b/postgis/geography.sql.in.c
@@ -456,6 +456,12 @@ CREATE OR REPLACE FUNCTION ST_Area(geography)
        AS 'MODULE_PATHNAME','geography_area_sphere'
        LANGUAGE 'C' IMMUTABLE STRICT;
 
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION ST_PointOutside(geography)
+       RETURNS geography
+       AS 'MODULE_PATHNAME','geography_point_outside'
+       LANGUAGE 'C' IMMUTABLE STRICT;
+
 -- ---------- ---------- ---------- ---------- ---------- ---------- ----------
 
 COMMIT;

diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 6c02ad9f774195915c7812f5ef802432ff0ce712..45efdd5f1efa3c7018cd493b60b5866f77a12d57 100644 (file)
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -159,11 +159,11 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS)
        area = lwgeom_area_sphere(lwgeom, gbox);
 
        /* Something went wrong... */
-       if( area < 0.0 )
-       {
-               elog(ERROR, "lwgeom_area_sphere returned area < 0.0");
-               PG_RETURN_NULL();
-       }
+//     if( area < 0.0 )
+//     {
+//             elog(ERROR, "lwgeom_area_sphere returned area < 0.0");
+//             PG_RETURN_NULL();
+//     }
 
        /* Currently normalizing with a fixed WGS84 radius, in future this
           should be the average radius of the SRID in play */
@@ -174,4 +174,40 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS)
        
        PG_RETURN_FLOAT8(area);
 
+}
+
+/*
+** geography_point_outside(GSERIALIZED *g) 
+** returns point outside the object
+*/
+PG_FUNCTION_INFO_V1(geography_point_outside);
+Datum geography_point_outside(PG_FUNCTION_ARGS)
+{
+       GBOX gbox;
+       GSERIALIZED *g = NULL;
+       GSERIALIZED *g_out = NULL;
+       size_t g_out_size;
+       LWPOINT *lwpoint = NULL;
+       POINT2D pt;
+       
+       /* Get our geometry object loaded into memory. */
+       g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+       
+       /* We need the bounding box to get an outside point for area algorithm */
+       if( ! gbox_from_gserialized(g, &gbox) )
+       {
+               elog(ERROR, "Error in gbox_from_gserialized calculation.");
+               PG_RETURN_NULL();
+       }
+       
+       /* Get an exterior point, based on this gbox */
+       gbox_pt_outside(gbox, &pt);
+       
+       lwpoint = make_lwpoint2d(4326, pt.x, pt.y);
+       
+       g_out = gserialized_from_lwgeom((LWGEOM*)lwpoint, 1, &g_out_size);
+       SET_VARSIZE(g_out, g_out_size);
+       
+       PG_RETURN_POINTER(g_out);
+
 }
\ No newline at end of file