From: Paul Ramsey Date: Sun, 18 Oct 2009 04:19:04 +0000 (+0000) Subject: ST_Length(geography) per #266 X-Git-Tag: 1.5.0b1~358 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=38adeece1377e135ed86051f79150a970da49816;p=postgis 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 32935df46..5714ab92c 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 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 9ebdbc686..678394354 100644 --- 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 808d8679e..83fee0daa 100644 --- 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 7108980f1..70df83119 100644 --- 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 ceb691e0e..2a2d2e806 100644 --- 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