ST_Length(geography) per #266

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

diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h
index 32935df466046bde49529bbec7c1f1ce09d03323..5714ab92c3327af1f56b1d8825a0538731e525df 100644 (file)
--- 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 the geodetic length of a lwgeom on the unit sphere. The result
+* will have to by multiplied by the real radius to get the real length.
+*/
+extern double lwgeom_length_sphere(LWGEOM *geom);
+
 /**
 * Calculate covers predicate for two lwgeoms on the sphere. Currently
 * only handles point-in-polygon.

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 9ebdbc6862eb2e8e5ad29bb9e87ea71cb5068b59..678394354000b688428ba47b8fa9a23a1296dc5f 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -2270,10 +2270,75 @@ int lwgeom_check_geodetic(const LWGEOM *geom)
        return LW_FALSE;
 }
 
+double ptarray_length_sphere(POINTARRAY *pa)
+{
+       GEOGRAPHIC_POINT a, b;
+       POINT2D p;
+       int i;
+       double length = 0.0;
+       
+       /* Return zero on non-sensical inputs */
+       if( ! pa || pa->npoints < 2 )
+               return 0.0;
+       
+       /* Initialize first point */
+       getPoint2d_p(pa, 0, &p);
+       geographic_point_init(p.x, p.y, &a);
 
+       for( i = 1; i < pa->npoints; i++ )
+       {               
+               getPoint2d_p(pa, i, &p);
+               geographic_point_init(p.x, p.y, &b);
+               
+               length += sphere_distance(a, b);
+                               
+               /* B gets incremented in the next loop, so we save the value here */
+               a = b;
+       }
+       return length;
+}
 
+double lwgeom_length_sphere(LWGEOM *geom)
+{
+       int type;
+       int i = 0;
+       double length = 0.0;
+       
+       assert(geom);
+       
+       /* No area in nothing */
+       if( lwgeom_is_empty(geom) )
+               return 0.0;
 
+       type = TYPE_GETTYPE(geom->type);
 
+       if ( type == POINTTYPE || type == MULTIPOINTTYPE )
+               return 0.0;
+               
+       if ( type == LINETYPE )
+               return ptarray_length_sphere(((LWLINE*)geom)->points);
+               
+       if ( type == POLYGONTYPE )
+       {
+               LWPOLY *poly = (LWPOLY*)geom;
+               for( i = 0; i < poly->nrings; i++ )
+               {
+                       length += ptarray_length_sphere(poly->rings[i]);
+               }
+               return length;
+       }
 
-
-
+       if( lwgeom_is_collection( type ) )
+       {
+               LWCOLLECTION *col = (LWCOLLECTION*)geom;
+               
+               for( i = 0; i < col->ngeoms; i++ )
+               {
+                       length += lwgeom_length_sphere(col->geoms[i]);
+               }
+               return length;
+       }
+       
+       lwerror("unsupported type passed to lwgeom_length_sphere");
+       return 0.0;
+}      

diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 808d8679e2f61f6c9caae020ad684e8377994daf..83fee0daa74e3451d51f729629289fbd8efc9739 100644 (file)
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -77,3 +77,4 @@ int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test
 double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside);
 double latitude_degrees_normalize(double lat);
 double longitude_degrees_normalize(double lon);
+double ptarray_length_sphere(POINTARRAY *pa);

diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c
index 7108980f10b06adc14f0469cf114dbe908b96668..70df8311925ce4293f2f1d21f4a1bbd6cb20d214 100644 (file)
--- a/postgis/geography.sql.in.c
+++ b/postgis/geography.sql.in.c
@@ -465,6 +465,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_Length(geography)
+       RETURNS float8
+       AS 'MODULE_PATHNAME','geography_length_sphere'
+       LANGUAGE 'C' IMMUTABLE STRICT;
+
 -- Availability: 1.5.0
 CREATE OR REPLACE FUNCTION ST_PointOutside(geography)
        RETURNS geography

diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index ceb691e0e493c2fdd2efca502ed6f3900619487e..2a2d2e80671ae04c89b14e1bd1ea28b81a5556bb 100644 (file)
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -25,11 +25,11 @@
 
 Datum geography_distance_sphere(PG_FUNCTION_ARGS);
 Datum geography_area_sphere(PG_FUNCTION_ARGS);
+Datum geography_length_sphere(PG_FUNCTION_ARGS);
 Datum geography_expand(PG_FUNCTION_ARGS);
 Datum geography_point_outside(PG_FUNCTION_ARGS);
 Datum geography_covers(PG_FUNCTION_ARGS);
 
-
 /*
 ** geography_distance_sphere(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance) 
 ** returns double distance in meters
@@ -161,11 +161,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 */
@@ -178,6 +178,42 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS)
 
 }
 
+/*
+** geography_length_sphere(GSERIALIZED *g) 
+** returns double length in meters
+*/
+PG_FUNCTION_INFO_V1(geography_length_sphere);
+Datum geography_length_sphere(PG_FUNCTION_ARGS)
+{
+       LWGEOM *lwgeom = NULL;
+       GSERIALIZED *g = NULL;
+       double length;
+       
+       /* Get our geometry object loaded into memory. */
+       g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+               
+       lwgeom = lwgeom_from_gserialized(g);
+       
+       /* Calculate the length */
+       length = lwgeom_length_sphere(lwgeom);
+
+       /* Something went wrong... */
+       if( length < 0.0 )
+       {
+               elog(ERROR, "geography_length_sphere returned length < 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 */
+       length = length * WGS84_RADIUS;
+
+       /* Clean up, but not all the way to the point arrays */
+       lwgeom_release(lwgeom);
+       
+       PG_RETURN_FLOAT8(length);
+}
+
 /*
 ** geography_point_outside(GSERIALIZED *g) 
 ** returns point outside the object