From d5e350a58796e140b66acf0aeaba128fbe6dad3f Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 8 Oct 2009 05:35:59 +0000 Subject: [PATCH] ST_Area(geography) implementation and SQL bindings. git-svn-id: http://svn.osgeo.org/postgis/trunk@4621 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/libgeom.h | 8 +- liblwgeom/lwgeodetic.c | 121 ++++++++++++++++++ liblwgeom/lwgeodetic.h | 1 + postgis/Makefile.in | 2 +- postgis/geography.sql.in.c | 8 +- ...phy_distance.c => geography_measurement.c} | 47 ++++++- 6 files changed, 181 insertions(+), 6 deletions(-) rename postgis/{geography_distance.c => geography_measurement.c} (76%) diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index 030621ced..f234e589b 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -376,14 +376,18 @@ extern int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox); */ extern double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, double tolerance); +/** +* Calculate the geodetic area of a lwgeom on the unit sphere. The result +* will have to by multiplied by the real radius^2 to get the real area. +*/ +extern double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox); + /** * New function to read doubles directly from the double* coordinate array * of an aligned lwgeom #POINTARRAY (built by de-serializing a #GSERIALIZED). */ extern int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point); - - /** * Calculate box and add values to gbox. Return #G_SUCCESS on success. */ diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 5bd031d31..d5fe5f5bc 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -552,6 +552,27 @@ double sphere_distance_cartesian(POINT3D s, POINT3D e) return acos(dot_product(s, e)); } +/** +* Computes the spherical excess of a spherical triangle defined by +* the three vectices A, B, C. Computes on the unit sphere (i.e., divides +* edge lengths by the radius, even if the radius is 1.0). The excess is +* signed based on the sign of the delta longitude of A and B. +* +* @param a The first triangle vertex. +* @param b The second triangle vertex. +* @param c The last triangle vertex. +* @return the signed spherical excess. +*/ +static double sphere_excess(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, GEOGRAPHIC_POINT c) +{ + double a_dist = sphere_distance(b, c); + double b_dist = sphere_distance(c, a); + double c_dist = sphere_distance(a, b); + double sign = (a.lon - b.lon) / fabs(a.lon - b.lon); + double ss = (a_dist + b_dist + c_dist) / 2.0; + double E = tan(ss/2.0) * tan((ss-a_dist)/2.0) * tan((ss-b_dist)/2.0) * tan((ss-c_dist)/2.0); + return 4.0 * atan(sqrt(E)) * sign; +} /** * Given two points on a unit sphere, calculate the direction from s to e. @@ -1269,6 +1290,34 @@ static void gbox_pt_outside(GBOX gbox, POINT3D *pt) return; } +/** +* Returns the area of the ring (ring must be closed) in square radians (surface of +* the sphere is 4*PI). +*/ +double ptarray_area_sphere(POINTARRAY *pa, POINT2D pt_outside) +{ + GEOGRAPHIC_POINT a, b, c; + POINT2D p; + int i; + double area = 0.0; + + /* Return zero on non-sensical inputs */ + if( ! pa || pa->npoints < 4 ) + return 0.0; + + geographic_point_init(pt_outside.x, pt_outside.y, &c); + + for( i = 1; i < pa->npoints; i++ ) + { + getPoint2d_p(pa, i-1, &p); + geographic_point_init(p.x, p.y, &a); + getPoint2d_p(pa, i, &p); + geographic_point_init(p.x, p.y, &b); + area += sphere_excess(a, b, c); + } + return fabs(area); +} + /** * This routine returns LW_TRUE if the point is inside the ring or on the boundary, LW_FALSE otherwise. @@ -1436,6 +1485,78 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t } +/** +* Calculate the area of an LWGEOM. Anything except POLYGON, MULTIPOLYGON +* and GEOMETRYCOLLECTION return zero immediately. Multi's recurse, polygons +* calculate external ring area and subtract internal ring area. A GBOX is +* required to calculate an outside point. +*/ +double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox) +{ + int type; + POINT3D p; + GEOGRAPHIC_POINT g; + POINT2D pt_outside; + + assert(lwgeom); + + /* No area in nothing */ + if( lwgeom_is_empty(lwgeom) ) + return 0.0; + + /* Read the geometry type number */ + type = TYPE_GETTYPE(lwgeom->type); + + /* Anything but polygons and collections returns zero */ + 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); + + /* Actually calculate area */ + if( type == POLYGONTYPE ) + { + LWPOLY *poly = (LWPOLY*)lwgeom; + int i; + double area = 0.0; + + /* Just in case there's no rings */ + if( poly->nrings < 1 ) + return 0.0; + + /* First, the area of the outer ring */ + area += ptarray_area_sphere(poly->rings[0], pt_outside); + + /* Subtract areas of inner rings */ + for( i = 1; i < poly->nrings; i++ ) + { + area -= ptarray_area_sphere(poly->rings[i], pt_outside); + } + return area; + } + + /* Recurse into sub-geometries to get area */ + if( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) + { + LWCOLLECTION *col = (LWCOLLECTION*)lwgeom; + int i; + double area = 0.0; + + for( i = 0; i < col->ngeoms; i++ ) + { + area += lwgeom_area_sphere(col->geoms[i], gbox); + } + return area; + } + + /* Shouldn't get here. */ + return 0.0; +} + + /** * Calculate the distance between two LWGEOMs, using the coordinates are * longitude and latitude. Return immediately when the calulated distance drops diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index 391544c61..64123c90f 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -75,5 +75,6 @@ void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g); int ptarray_point_in_ring_winding(POINTARRAY *pa, POINT2D pt_to_test); int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test); 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); diff --git a/postgis/Makefile.in b/postgis/Makefile.in index 69522c79e..667640ec0 100644 --- a/postgis/Makefile.in +++ b/postgis/Makefile.in @@ -52,7 +52,7 @@ PG_OBJS=lwgeom_pg.o \ geography_inout.o \ geography_gist.o \ geography_estimate.o \ - geography_distance.o + geography_measurement.o # Objects to build using PGXS OBJS=$(PG_OBJS) diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 3adae3214..697f3f9b2 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -449,7 +449,13 @@ CREATE OR REPLACE FUNCTION ST_DWithin(geography, geography, float8) RETURNS boolean AS 'SELECT $1 && _ST_Expand($2,$3) AND $2 && _ST_Expand($1,$3) AND _ST_Distance($1, $2, $3) < $3' LANGUAGE 'SQL' IMMUTABLE; - + +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_Area(geography) + RETURNS float8 + AS 'MODULE_PATHNAME','geography_area_sphere' + LANGUAGE 'C' IMMUTABLE STRICT; + -- ---------- ---------- ---------- ---------- ---------- ---------- ---------- COMMIT; diff --git a/postgis/geography_distance.c b/postgis/geography_measurement.c similarity index 76% rename from postgis/geography_distance.c rename to postgis/geography_measurement.c index b6455741c..656694b35 100644 --- a/postgis/geography_distance.c +++ b/postgis/geography_measurement.c @@ -24,6 +24,7 @@ #include "geography.h" /* For utility functions. */ Datum geography_distance_sphere(PG_FUNCTION_ARGS); +Datum geography_area_sphere(PG_FUNCTION_ARGS); Datum geography_expand(PG_FUNCTION_ARGS); @@ -38,8 +39,6 @@ Datum geography_distance_sphere(PG_FUNCTION_ARGS) LWGEOM *lwgeom2 = NULL; GBOX gbox1; GBOX gbox2; - GIDX *gidx1 = gidx_new(3); - GIDX *gidx2 = gidx_new(3); GSERIALIZED *g1 = NULL; GSERIALIZED *g2 = NULL; double tolerance; @@ -132,3 +131,47 @@ Datum geography_expand(PG_FUNCTION_ARGS) } +/* +** geography_area_sphere(GSERIALIZED *g) +** returns double area in meters square +*/ +PG_FUNCTION_INFO_V1(geography_area_sphere); +Datum geography_area_sphere(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom = NULL; + GBOX gbox; + GSERIALIZED *g = NULL; + double area; + + /* 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(); + } + + lwgeom = lwgeom_from_gserialized(g); + + /* Calculate the area */ + area = lwgeom_area_sphere(lwgeom, gbox); + + /* Something went wrong... */ + if( area < 0.0 ) + { + elog(ERROR, "Error in geography_distance_sphere calculation."); + PG_RETURN_NULL(); + } + + /* Currently normalizing with a fixed WGS84 radius, in future this + should be the average radius of the SRID in play */ + area = area * WGS84_RADIUS * WGS84_RADIUS; + + /* Clean up, but not all the way to the point arrays */ + lwgeom_release(lwgeom); + + PG_RETURN_FLOAT8(area); + +} \ No newline at end of file -- 2.50.1