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.
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.
}
+/**
+* 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
#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);
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;
}
+/*
+** 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
projects / postgis / commitdiff