From: Paul Ramsey Date: Tue, 3 Nov 2009 00:36:02 +0000 (+0000) Subject: First cut of ST_Area(geography) on spheroid. Currently not default, use ST_Area(geog... X-Git-Tag: 1.5.0b1~314 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=fba1c9c339873aa82face5b468d0137897a5e82b;p=postgis First cut of ST_Area(geography) on spheroid. Currently not default, use ST_Area(geog, true) to enable it. Beware of limitations over poles and eequator. git-svn-id: http://svn.osgeo.org/postgis/trunk@4719 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 5b7881450..da3a3981f 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -28,7 +28,7 @@ CU_pSuite register_geodetic_suite(void) } if ( - (NULL == CU_add_test(pSuite, "test_signum()", test_signum)) || +/* (NULL == CU_add_test(pSuite, "test_signum()", test_signum)) || (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) || (NULL == CU_add_test(pSuite, "test_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) || (NULL == CU_add_test(pSuite, "test_clairaut()", test_clairaut)) || @@ -37,7 +37,8 @@ CU_pSuite register_geodetic_suite(void) (NULL == CU_add_test(pSuite, "test_edge_distance_to_edge()", test_edge_distance_to_edge)) || (NULL == CU_add_test(pSuite, "test_lwgeom_distance_sphere()", test_lwgeom_distance_sphere)) || (NULL == CU_add_test(pSuite, "test_ptarray_point_in_ring()", test_ptarray_point_in_ring)) || - (NULL == CU_add_test(pSuite, "test_spheroid_distance()", test_spheroid_distance)) + (NULL == CU_add_test(pSuite, "test_spheroid_distance()", test_spheroid_distance)) || */ + (NULL == CU_add_test(pSuite, "test_spheroid_area()", test_spheroid_area)) ) { CU_cleanup_registry(); @@ -609,4 +610,34 @@ void test_spheroid_distance(void) d = spheroid_distance(g1, g2, s); CU_ASSERT_DOUBLE_EQUAL(d, 10001965.7295318, 0.001); +} + +void test_spheroid_area(void) +{ + LWGEOM *lwg; + GBOX gbox; + double a1, a2; + SPHEROID s; + + /* Init to WGS84 */ + spheroid_init(&s, 6378137.0, 6356752.314245179498); + + gbox.flags = gflags(0, 0, 1); + + /* Medford lot test polygon */ + lwg = lwgeom_from_ewkt("POLYGON((-122.848227067007 42.5007249610493,-122.848309475585 42.5007179884263,-122.848327688675 42.500835880696,-122.848245279942 42.5008428533324,-122.848227067007 42.5007249610493))", PARSER_CHECK_NONE); + lwgeom_calculate_gbox_geodetic(lwg, &gbox); + a1 = lwgeom_area_sphere(lwg, gbox, s); + a2 = lwgeom_area_spheroid(lwg, gbox, s); + //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); + CU_ASSERT_DOUBLE_EQUAL(a1, a2, 0.2); + + /* Big-ass polygon */ + lwg = lwgeom_from_ewkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", PARSER_CHECK_NONE); + lwgeom_calculate_gbox_geodetic(lwg, &gbox); + a1 = lwgeom_area_sphere(lwg, gbox, s); + a2 = lwgeom_area_spheroid(lwg, gbox, s); + //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); + CU_ASSERT_DOUBLE_EQUAL(a1, a2, 100000000.0); + } \ No newline at end of file diff --git a/liblwgeom/cunit/cu_geodetic.h b/liblwgeom/cunit/cu_geodetic.h index 102957a25..594fda1c1 100644 --- a/liblwgeom/cunit/cu_geodetic.h +++ b/liblwgeom/cunit/cu_geodetic.h @@ -36,3 +36,4 @@ void test_ptarray_point_in_ring_winding(void); void test_lwgeom_distance_sphere(void); void test_ptarray_point_in_ring(void); void test_spheroid_distance(void); +void test_spheroid_area(void); diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c index 55813410a..9ba7b2082 100644 --- a/liblwgeom/g_box.c +++ b/liblwgeom/g_box.c @@ -437,7 +437,7 @@ static int lwcircle_calculate_gbox(POINT4D p1, POINT4D p2, POINT4D p3, GBOX *gbo return G_SUCCESS; } -static int ptarray_calculate_gbox( POINTARRAY *pa, GBOX *gbox ) +int ptarray_calculate_gbox( POINTARRAY *pa, GBOX *gbox ) { int i; POINT4D p; diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index 43618d190..3e66933ec 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -383,10 +383,16 @@ extern void spheroid_init(SPHEROID *s, double a, double b); extern double lwgeom_distance_spheroid(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX gbox1, GBOX gbox2, SPHEROID spheroid, 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. +* Calculate the geodetic area of a lwgeom on the sphere. The result +* will be multiplied by the average radius of the supplied spheroid. */ -extern double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID s); +extern double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid); + +/** +* Calculate the geodetic area of a lwgeom on the spheroid. The result +* will have the squared units of the spheroid axes. +*/ +extern double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid); /** * Calculate the geodetic length of a lwgeom on the unit sphere. The result @@ -407,10 +413,15 @@ extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwge extern int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point); /** -* Calculate box and add values to gbox. Return #G_SUCCESS on success. +* Calculate geodetic (x/y/z) box and add values to gbox. Return #G_SUCCESS on success. */ extern int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox); +/** +* Calculate box (x/y) and add values to gbox. Return #G_SUCCESS on success. +*/ +extern int ptarray_calculate_gbox( POINTARRAY *pa, GBOX *gbox ); + /** * Calculate a spherical point that falls outside the geocentric gbox */ diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index d80b6e358..97b8ada09 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -165,7 +165,7 @@ void point_rad2deg(GEOGRAPHIC_POINT *p) /** * Shift a point around by a number of radians */ -static void point_shift(GEOGRAPHIC_POINT *p, double shift) +void point_shift(GEOGRAPHIC_POINT *p, double shift) { double lon = p->lon + shift; if( lon > M_PI ) @@ -365,7 +365,7 @@ void y_to_z(POINT3D *p) } -static inline int crosses_dateline(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e) +int crosses_dateline(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e) { double sign_s = signum(s.lon); double sign_e = signum(e.lon); @@ -1639,7 +1639,7 @@ static double ptarray_distance_spheroid(POINTARRAY *pa1, POINTARRAY *pa2, SPHERO * calculate external ring area and subtract internal ring area. A GBOX is * required to calculate an outside point. */ -double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) +double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) { int type; POINT2D pt_outside; @@ -1693,7 +1693,7 @@ double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) for( i = 0; i < col->ngeoms; i++ ) { - area += lwgeom_area_spheroid(col->geoms[i], gbox, spheroid); + area += lwgeom_area_sphere(col->geoms[i], gbox, spheroid); } return area; } diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index 5738642c1..cd411249f 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -79,7 +79,11 @@ double latitude_degrees_normalize(double lat); double longitude_degrees_normalize(double lon); double ptarray_length_spheroid(POINTARRAY *pa, SPHEROID s); int geographic_point_equals(GEOGRAPHIC_POINT g1, GEOGRAPHIC_POINT g2); +int crosses_dateline(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e); +void point_shift(GEOGRAPHIC_POINT *p, double shift); /* ** Prototypes for spheroid functions. */ double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid); +double spheroid_direction(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, SPHEROID spheroid); +int spheroid_project(GEOGRAPHIC_POINT r, SPHEROID spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g); diff --git a/liblwgeom/lwspheroid.c b/liblwgeom/lwspheroid.c index 8c3eaaf21..e00c278e9 100644 --- a/liblwgeom/lwspheroid.c +++ b/liblwgeom/lwspheroid.c @@ -2,6 +2,9 @@ #define POW2(x) ((x)*(x)) +/** +* Initialize spheroid object based on major and minor axis +*/ void spheroid_init(SPHEROID *s, double a, double b) { s->a = a; @@ -27,6 +30,20 @@ static double spheroid_big_b(double u2) return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2))); } +/** +* Computes the shortest distance along the surface of the spheroid +* between two points. Based on Vincenty's formula for the geodetic +* inverse problem as described in "Geocentric Datum of Australia +* Technical Manual", Chapter 4. Tested against: +* http://mascot.gdbc.gov.bc.ca/mascot/util1a.html +* and +* http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp +* +* @param a - location of first point. +* @param b - location of second point. +* @param s - spheroid to calculate on +* @return spheroidal distance between a and b in spheroid units. +*/ double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid) { double lambda = (b.lon - a.lon); @@ -90,7 +107,7 @@ double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID sphero (cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m)))); i++; } - while ( i < 999 && lambda != 0.0 && fabs((last_lambda - lambda)/lambda) > 1.0e-9 ); + while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) ); u2 = spheroid_mu2(alpha, spheroid); big_a = spheroid_big_a(u2); @@ -110,3 +127,655 @@ double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID sphero return distance; } +/** +* Computes the direction of the geodesic joining two points on +* the spheroid. Based on Vincenty's formula for the geodetic +* inverse problem as described in "Geocentric Datum of Australia +* Technical Manual", Chapter 4. Tested against: +* http://mascot.gdbc.gov.bc.ca/mascot/util1a.html +* and +* http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp +* +* @param r - location of first point +* @param s - location of second point +* @return azimuth of line joining r and s +*/ +double spheroid_direction(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, SPHEROID spheroid) +{ + int i = 0; + double lambda = s.lon - r.lon; + double omf = 1 - spheroid.f; + double u1 = atan(omf * tan(r.lat)); + double cos_u1 = cos(u1); + double sin_u1 = sin(u1); + double u2 = atan(omf * tan(s.lat)); + double cos_u2 = cos(u2); + double sin_u2 = sin(u2); + + double omega = lambda; + double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda; + double sin_alpha, cos_alphasq, C, alphaFD; + do + { + sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) + + POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda))); + sin_sigma = sqrt(sqr_sin_sigma); + cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda); + sigma = atan2(sin_sigma, cos_sigma); + sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma); + + /* Numerical stability issue, ensure asin is not NaN */ + if( sin_alpha > 1.0 ) + alpha = M_PI_2; + else if( sin_alpha < -1.0 ) + alpha = -1.0 * M_PI_2; + else + alpha = asin(sin_alpha); + + alpha = asin(sin_alpha); + cos_alphasq = POW2(cos(alpha)); + cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq); + + /* Numerical stability issue, cos2 is in range */ + if( cos2_sigma_m > 1.0 ) + cos2_sigma_m = 1.0; + if( cos2_sigma_m < -1.0 ) + cos2_sigma_m = -1.0; + + C = (spheroid.f / 16.0) * cos_alphasq * (4.0 + spheroid.f * (4.0 - 3.0 * cos_alphasq)); + last_lambda = lambda; + lambda = omega + (1.0 - C) * spheroid.f * sin(alpha) * (sigma + C * sin(sigma) * + (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m)))); + i++; + } + while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) ); + + alphaFD = atan2((cos_u2 * sin(lambda)), + (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda))); + if (alphaFD < 0.0) + { + alphaFD = alphaFD + 2.0 * M_PI; + } + if (alphaFD > 2.0 * M_PI) + { + alphaFD = alphaFD - 2.0 * M_PI; + } + return alphaFD; +} + + +/** +* Given a location, an azimuth and a distance, computes the +* location of the projected point. Based on Vincenty's formula +* for the geodetic direct problem as described in "Geocentric +* Datum of Australia Technical Manual", Chapter 4. Tested against: +* http://mascot.gdbc.gov.bc.ca/mascot/util1b.html +* and +* http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp +* +* @param r - location of first point. +* @param distance - distance in meters. +* @param azimuth - azimuth in radians. +* @return s - location of projected point. +*/ +int spheroid_project(GEOGRAPHIC_POINT r, SPHEROID spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g) +{ + double omf = 1 - spheroid.f; + double tan_u1 = omf * tan(r.lat); + double u1 = atan(tan_u1); + double sigma, last_sigma, delta_sigma, two_sigma_m; + double sigma1, sin_alpha, alpha, cos_alphasq; + double u2, A, B; + double lat2, lambda, lambda2, C, omega; + int i = 0; + + if (azimuth < 0.0) { + azimuth = azimuth + M_PI * 2.0; + } + if (azimuth > (PI * 2.0)) { + azimuth = azimuth - M_PI * 2.0; + } + + sigma1 = atan2(tan_u1, cos(azimuth)); + sin_alpha = cos(u1) * sin(azimuth); + alpha = asin(sin_alpha); + cos_alphasq = 1.0 - POW2(sin_alpha); + + u2 = spheroid_mu2(alpha, spheroid); + A = spheroid_big_a(u2); + B = spheroid_big_b(u2); + + sigma = (distance / (spheroid.b * A)); + do + { + two_sigma_m = 2.0 * sigma1 + sigma; + delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m)))))); + last_sigma = sigma; + sigma = (distance / (spheroid.b * A)) + delta_sigma; + i++; + } + while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9); + + lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) * + cos(azimuth)), (omf * sqrt(POW2(sin_alpha) + + POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) * + cos(azimuth))))); + lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) - + sin(u1) * sin(sigma) * cos(azimuth))); + C = (spheroid.f / 16.0) * cos_alphasq * (4.0 + spheroid.f * (4.0 - 3.0 * cos_alphasq)); + omega = lambda - (1.0 - C) * spheroid.f * sin_alpha * (sigma + C * sin(sigma) * + (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m))))); + lambda2 = r.lon + omega; + g->lat = lat2; + g->lon = lambda2; + return G_SUCCESS; +} + + +static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, SPHEROID spheroid) +{ + return spheroid.a / (sqrt(1.0 - spheroid.e_sq * POW2(sin(latitude)))); +} + +static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, SPHEROID spheroid) +{ + return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid) + * cos(latitude) + * deltaLongitude; +} + + +/** +* Computes the area on the spheroid of a box bounded by meridians and +* parallels. The box is defined by two points, the South West corner +* and the North East corner. Formula based on Bagratuni 1967. +* +* @param southWestCorner - lower left corner of bounding box. +* @param northEastCorner - upper right corner of bounding box. +* @return area in square meters. +*/ +static double spheroid_boundingbox_area(GEOGRAPHIC_POINT southWestCorner, GEOGRAPHIC_POINT northEastCorner, SPHEROID spheroid) +{ + double z0 = (northEastCorner.lon - southWestCorner.lon) * POW2(spheroid.b) / 2.0; + double e = sqrt(spheroid.e_sq); + double sinPhi1 = sin(southWestCorner.lat); + double sinPhi2 = sin(northEastCorner.lat); + double t1p1 = sinPhi1 / (1.0 - spheroid.e_sq * sinPhi1 * sinPhi1); + double t1p2 = sinPhi2 / (1.0 - spheroid.e_sq * sinPhi2 * sinPhi2); + double oneOver2e = 1.0 / (2.0 * e); + double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1)); + double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2)); + return z0 * (t1p2 + t2p2) - z0 * (t1p1 + t2p1); +} + +/** +* This function doesn't work for edges crossing the dateline or in the southern +* hemisphere. Points are pre-conditioned in ptarray_area_spheroid. +*/ +static double spheroid_striparea(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, double latitude_min, SPHEROID spheroid) +{ + GEOGRAPHIC_POINT A = a; + GEOGRAPHIC_POINT B = b; + GEOGRAPHIC_POINT mL, nR; + double deltaLng, baseArea, topArea; + double bE, tE, ratio, sign; + + mL.lat = latitude_min; + mL.lon = FP_MIN(A.lon, B.lon); + nR.lat = FP_MIN(A.lat, B.lat); + nR.lon = FP_MAX(A.lon, B.lon); + baseArea = spheroid_boundingbox_area(mL, nR, spheroid); + + mL.lat = FP_MIN(A.lat, B.lat); + mL.lon = FP_MIN(A.lon, B.lon); + nR.lat = FP_MAX(A.lat, B.lat); + nR.lon = FP_MAX(A.lon, B.lon); + topArea = spheroid_boundingbox_area(mL, nR, spheroid); + + deltaLng = B.lon - A.lon; + bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid); + tE = spheroid_parallel_arc_length(B.lat, deltaLng, spheroid); + ratio = (bE + tE)/tE; + sign = signum(B.lon - A.lon); + return (baseArea + topArea / ratio) * sign; +} + +static double ptarray_area_spheroid(POINTARRAY *pa, SPHEROID spheroid) +{ + GEOGRAPHIC_POINT a, b; + POINT2D p; + int i; + double area = 0.0; + GBOX gbox2d; + int in_south = LW_FALSE; + double delta_lon_tolerance; + double latitude_min; + + gbox2d.flags = gflags(0, 0, 0); + + /* Return zero on non-sensical inputs */ + if( ! pa || pa->npoints < 4 ) + return 0.0; + + /* Get the raw min/max values for the latitudes */ + ptarray_calculate_gbox(pa, &gbox2d); + + if( signum(gbox2d.ymin) != signum(gbox2d.ymax) ) + lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator"); + + /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */ + if( gbox2d.ymax < 0.0 ) + in_south = LW_TRUE; + + LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax); + + /* Tolerance for strip area calculation */ + if( in_south ) + { + delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0; + latitude_min = deg2rad(fabs(gbox2d.ymax)); + } + else + { + delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0; + latitude_min = deg2rad(gbox2d.ymin); + } + + /* Initialize first point */ + getPoint2d_p(pa, 0, &p); + geographic_point_init(p.x, p.y, &a); + + for( i = 1; i < pa->npoints; i++ ) + { + GEOGRAPHIC_POINT a1, b1; + double strip_area = 0.0; + double delta_lon = 0.0; + LWDEBUGF(4, "edge #%d", i); + + getPoint2d_p(pa, i, &p); + geographic_point_init(p.x, p.y, &b); + + a1 = a; b1 = b; + + /* Flip into north if in south */ + if( in_south ) + { + a1.lat = -1.0 * a1.lat; + b1.lat = -1.0 * b1.lat; + } + + LWDEBUGF(4, "in_south %d", in_south); + + if( crosses_dateline(a, b) ) + { + double shift; + + if( a1.lon > 0.0 ) + shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */ + else + shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */ + + point_shift(&a1, shift); + point_shift(&b1, shift); + } + + LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(a, b) ); + + delta_lon = fabs(b1.lon - a1.lon); + + LWDEBUGF(4, "(%.18g %.18g) (%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon); + LWDEBUGF(4, "delta_lon %.18g", delta_lon); + LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance); + + if( delta_lon > 0.0 ) + { + if( delta_lon < delta_lon_tolerance ) + { + strip_area = spheroid_striparea(a1, b1, latitude_min, spheroid); + LWDEBUGF(4, "strip_area %.12g", strip_area); + area += strip_area; + } + else + { + GEOGRAPHIC_POINT p, q; + double step = floor(delta_lon / delta_lon_tolerance); + double distance = spheroid_distance(a1, b1, spheroid); + double pDistance = 0.0; + int j = 0; + LWDEBUGF(4, "step %.18g", step); + LWDEBUGF(4, "distance %.18g", distance); + step = distance / step; + LWDEBUGF(4, "step %.18g", step); + p = a1; + while (pDistance < (distance - step * 1.01)) + { + double azimuth = spheroid_direction(p, b1, spheroid); + j++; + LWDEBUGF(4, " iteration %d", j); + LWDEBUGF(4, " azimuth %.12g", azimuth); + pDistance = pDistance + step; + LWDEBUGF(4, " pDistance %.12g", pDistance); + spheroid_project(p, spheroid, step, azimuth, &q); + strip_area = spheroid_striparea(p, q, latitude_min, spheroid); + LWDEBUGF(4, " strip_area %.12g", strip_area); + area += strip_area; + LWDEBUGF(4, " area %.12g", area); + p.lat = q.lat; + p.lon = q.lon; + } + strip_area = spheroid_striparea(p, b1, latitude_min, spheroid); + area += strip_area; + } + } + + /* B gets incremented in the next loop, so we save the value here */ + a = b; + } + return fabs(area); +} +/** +* 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 check relationship to equator an outside point. +* WARNING: Does NOT WORK for polygons over equator or pole. +*/ +double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) +{ + int type; + + 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; + + /* 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_spheroid(poly->rings[0], spheroid); + + /* Subtract areas of inner rings */ + for( i = 1; i < poly->nrings; i++ ) + { + area -= ptarray_area_spheroid(poly->rings[i], spheroid); + } + 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_spheroid(col->geoms[i], gbox, spheroid); + } + return area; + } + + /* Shouldn't get here. */ + return 0.0; +} + + + +#if 0 + +private double areaOnSpheroid(KMLOut kml) { + int date = 1900; + boolean inSouthernHemisphere = false; + if(this.getMaxLat() < 0.0) { + inSouthernHemisphere = true; + } + double deltaLngTol = 0.0; + if(inSouthernHemisphere) { + deltaLngTol = (90.0 / ((Math.abs(this.getMinLat()) * 180 / Math.PI) / 8.0) - 2.0) / 10000.0; + } else { + deltaLngTol = (90.0 / ((Math.abs(this.getMaxLat()) * 180 / Math.PI) / 8.0) - 2.0) / 10000.0; + } + double area = 0.0; + double sa = 0.0; + double minLat = this.getMinLat(); + if(inSouthernHemisphere) { + minLat = Math.abs(this.getMaxLat()); + } + for (int ring = 0; ring < this.numberOfRings; ring++) { + EdgeIterator ei = new EdgeIterator(this, ring); + double ringArea = 0.0; + while (ei.hasMoreEdgesOnRing(ring)) { + GreatCircleEdge gce = ei.getNextEdge(); + if (gce.inSH()) { + gce.shiftIntoNH(); + } + double deltaLng = Math.abs(gce.getEnd().getLng() - gce.getStart().getLng()); + if (deltaLng < deltaLngTol) { + sa = spheroid.stripArea(gce, minLat); + ringArea = ringArea + sa; + } else { + double step = Math.floor(deltaLng / deltaLngTol); + double distance = spheroid.distance(gce.getStart(), gce.getEnd()); + step = distance / step; + double pDistance = 0.0; + Geographic p = new Geographic(gce.getStart()); + Geographic q = null; + while (pDistance < (distance - step * 1.01)) { + pDistance = pDistance + step; + double azmuth = spheroid.direction(p, gce.getEnd()); + q = spheroid.project(p, step, azmuth); + if (kml != null) { + date = addPlacemark(kml, p, q, minLat, date); + } + sa = spheroid.stripArea(p, q, minLat); + ringArea = ringArea + sa; + p.setLat(q.getLat()); + p.setLng(q.getLng()); + } + sa = spheroid.stripArea(p, gce.getEnd(), minLat); + if (kml != null) { + date = addPlacemark(kml, p, gce.getEnd(), minLat, date); + } + ringArea = ringArea + sa; + } + } + if (ring == 0) { + area = Math.abs(ringArea); + } else { + area = area - Math.abs(ringArea); + } + } + return area; +} +public double stripArea(Geographic _A, Geographic _B, double minLatitude) { + Geographic A = null; + Geographic B = null; + if(signum(_A.getLng())/signum(_B.getLng()) < 0.0) { + if(_A.getLng() > 0.0) { + double deltaLng = 2.0 * (Math.PI - _A.getLng()); + A = new Geographic(_A.lat, (_A.getLng() + deltaLng) - 2.0*Math.PI); + B = new Geographic(_B.getLat(), _B.getLng() + deltaLng); + } else { + double deltaLng = 2.0 * (Math.PI - _B.getLng()); + B = new Geographic(_B.getLat(), (_B.getLng() + deltaLng) - 2.0*Math.PI); + A = new Geographic(_A.getLat(), _A.getLng() + deltaLng); + } + } else { + A = new Geographic(_A); + B = new Geographic(_B); + } + Geographic mL = new Geographic(minLatitude, Math.min(A.getLng(), B.getLng())); + Geographic nR = new Geographic(Math.min(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng())); + double baseArea = this.boundingBoxArea(mL, nR); + mL = new Geographic(Math.min(A.getLat(), B.getLat()), Math.min(A.getLng(), B.getLng())); + nR = new Geographic(Math.max(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng())); + double topArea = this.boundingBoxArea(mL, nR); + double deltaLongitude = B.getLng() - A.getLng(); + double bE = this.paralleArcLength(A.getLat(), deltaLongitude); + double tE = this.paralleArcLength(B.getLat(), deltaLongitude); + double ratio = (bE + tE)/tE; + double sign = Math.signum(B.getLng() - A.getLng()); + return (baseArea + topArea / ratio) * sign; +} + +/** + * Given a location, an azimuth and a distance, computes the + * location of the projected point. Based on Vincenty's formula + * for the geodetic direct problem as described in "Geocentric + * Datum of Australia Technical Manual", Chapter 4. Tested against: + * http://mascot.gdbc.gov.bc.ca/mascot/util1b.html + * and + * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp + * + * @param r - location of first point. + * @param distance - distance in meters. + * @param azimuth - azimuth in radians. + * @return s - location of projected point. + */ +public Geographic project(Geographic r, double distance, double azimuth) { + if (azimuth < 0.0) { + azimuth = azimuth + PI * 2.0; + } + if (azimuth > (Math.PI * 2.0)) { + azimuth = azimuth - Math.PI * 2.0; + } + + double tanU1 = omf * Math.tan(r.getLat()); + double u1 = Math.atan(tanU1); + double sigma1 = Math.atan2(tanU1, Math.cos(azimuth)); + double sinAlpha = Math.cos(u1) * Math.sin(azimuth); + double alpha = Math.asin(sinAlpha); + double cosAlphaSq = 1.0 - sinAlpha * sinAlpha; + + double u2 = mu2(alpha); + double A = bigA(u2); + double B = bigB(u2); + + double sigma = (distance / (b * A)); + double lastSigma, twoSigmaM; + do { + twoSigmaM = 2.0 * sigma1 + sigma; + double deltaSigma = B * Math.sin(sigma) * (Math.cos(twoSigmaM) + (B / 4.0) * (Math.cos(sigma) * (-1.0 + 2.0 * power(Math.cos(twoSigmaM), 2) - (B / 6.0) * Math.cos(twoSigmaM) * (-3.0 + 4.0 * power(Math.sin(sigma), 2)) * (-3.0 + 4.0 * power(Math.cos(twoSigmaM), 2))))); + lastSigma = sigma; + sigma = (distance / (b * A)) + deltaSigma; + } while (Math.abs((lastSigma - sigma) / sigma) > 1.0e-9); + + double lat2 = Math.atan2((Math.sin(u1) * Math.cos(sigma) + Math.cos(u1) * Math.sin(sigma) * + Math.cos(azimuth)), (omf * Math.sqrt(power(sinAlpha, 2) + + power(Math.sin(u1) * Math.sin(sigma) - Math.cos(u1) * Math.cos(sigma) * + Math.cos(azimuth), 2)))); + double lambda = Math.atan2((Math.sin(sigma) * Math.sin(azimuth)), (Math.cos(u1) * Math.cos(sigma) - + Math.sin(u1) * Math.sin(sigma) * Math.cos(azimuth))); + double C = (f / 16.0) * cosAlphaSq * (4.0 + f * (4.0 - 3.0 * cosAlphaSq)); + double omega = lambda - (1.0 - C) * f * sinAlpha * (sigma + C * Math.sin(sigma) * + (Math.cos(twoSigmaM) + C * Math.cos(sigma) * (-1.0 + 2.0 * power(Math.cos(twoSigmaM), 2)))); + double lambda2 = r.getLng() + omega; + Geographic s = new Geographic(lat2, lambda2); + return s; +} + +/** + * Computes the direction of the geodesic joining two points on + * the spheroid. Based on Vincenty's formula for the geodetic + * inverse problem as described in "Geocentric Datum of Australia + * Technical Manual", Chapter 4. Tested against: + * http://mascot.gdbc.gov.bc.ca/mascot/util1a.html + * and + * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp + * + * @param r - location of first point. + * @param s - location of second point. + * @return azimuth of line joining r and s in meters. + */ +public double direction(Geographic r, Geographic s) { + double lambda = s.getLng() - r.getLng(); + double u1 = Math.atan(omf * Math.tan(r.getLat())); + double cosU1 = Math.cos(u1); + double sinU1 = Math.sin(u1); + double u2 = Math.atan(omf * Math.tan(s.getLat())); + double cosU2 = Math.cos(u2); + double sinU2 = Math.sin(u2); + + double omega = lambda; + double alpha, sigma, sinSigma, cosSigma, cos2SigmaM, sqrSinSigma, lastLambda; + do { + sqrSinSigma = power(cosU2 * Math.sin(lambda), 2) + + power((cosU1 * sinU2 - sinU1 * cosU2 * Math.cos(lambda)), 2); + sinSigma = Math.sqrt(sqrSinSigma); + cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * Math.cos(lambda); + sigma = Math.atan2(sinSigma, cosSigma); + double sinAlpha = cosU1 * cosU2 * Math.sin(lambda) / Math.sin(sigma); + alpha = Math.asin(sinAlpha); + double cosAlphaSq = power(Math.cos(alpha), 2); + cos2SigmaM = Math.cos(sigma) - (2.0 * sinU1 * sinU2 / cosAlphaSq); + double C = (f / 16.0) * cosAlphaSq * (4.0 + f * (4.0 - 3.0 * cosAlphaSq)); + lastLambda = lambda; + lambda = omega + (1.0 - C) * f * Math.sin(alpha) * (sigma + C * Math.sin(sigma) * + (cos2SigmaM + C * Math.cos(sigma) * (-1.0 + 2.0 * power(cos2SigmaM, 2)))); + } while ((lambda != 0) && (Math.abs((lastLambda - lambda) / lambda) > 1.0e-9)); + double alphaFD = Math.atan2((cosU2 * Math.sin(lambda)), + (cosU1 * sinU2 - sinU1 * cosU2 * Math.cos(lambda))); + if (alphaFD < 0.0) { + alphaFD = alphaFD + twoPI; + } + if (alphaFD > twoPI) { + alphaFD = alphaFD - twoPI; + } + return alphaFD; +} + + + +public double stripArea(Geographic _A, Geographic _B, double minLatitude) { + Geographic A = null; + Geographic B = null; + if(Math.signum(_A.getLng())/Math.signum(_B.getLng()) < 0.0) { + if(_A.getLng() > 0.0) { + double deltaLng = 2.0 * (Math.PI - _A.getLng()); + A = new Geographic(_A.getLat(), (_A.getLng() + deltaLng) - 2.0*Math.PI); + B = new Geographic(_B.getLat(), _B.getLng() + deltaLng); + } else { + double deltaLng = 2.0 * (Math.PI - _B.getLng()); + B = new Geographic(_B.getLat(), (_B.getLng() + deltaLng) - 2.0*Math.PI); + A = new Geographic(_A.getLat(), _A.getLng() + deltaLng); + } + } else { + A = new Geographic(_A); + B = new Geographic(_B); + } + Geographic mL = new Geographic(minLatitude, Math.min(A.getLng(), B.getLng())); + Geographic nR = new Geographic(Math.min(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng())); + double baseArea = this.boundingBoxArea(mL, nR); + mL = new Geographic(Math.min(A.getLat(), B.getLat()), Math.min(A.getLng(), B.getLng())); + nR = new Geographic(Math.max(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng())); + double topArea = this.boundingBoxArea(mL, nR); + double deltaLongitude = B.getLng() - A.getLng(); + double bE = this.paralleArcLength(A.getLat(), deltaLongitude); + double tE = this.paralleArcLength(B.getLat(), deltaLongitude); + double ratio = (bE + tE)/tE; + double sign = Math.signum(B.getLng() - A.getLng()); + return (baseArea + topArea / ratio) * sign; +} + + + +#endif + diff --git a/postgis/geography.h b/postgis/geography.h index dc8d6d406..cc0975e1d 100644 --- a/postgis/geography.h +++ b/postgis/geography.h @@ -18,7 +18,8 @@ */ #define WGS84_MAJOR_AXIS 6378137.0 -#define WGS84_MINOR_AXIS 6356752.314245179498 +#define WGS84_INVERSE_FLATTENING 298.257223563 +#define WGS84_MINOR_AXIS (WGS84_MAJOR_AXIS - WGS84_MAJOR_AXIS / WGS84_INVERSE_FLATTENING) #define WGS84_RADIUS ((2.0 * WGS84_MAJOR_AXIS + WGS84_MINOR_AXIS ) / 3.0) diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 233532158..15943dbda 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -479,11 +479,11 @@ CREATE OR REPLACE FUNCTION ST_Area(geography, boolean) AS 'MODULE_PATHNAME','geography_area' LANGUAGE 'C' IMMUTABLE STRICT; --- Currently defaulting to spheroid calculations +-- Currently defaulting to sphere calculations -- Availability: 1.5.0 CREATE OR REPLACE FUNCTION ST_Area(geography) RETURNS float8 - AS 'SELECT ST_Area($1, true)' + AS 'SELECT ST_Area($1, false)' LANGUAGE 'SQL' IMMUTABLE STRICT; -- Availability: 1.5.0 diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index 5f8282036..c26fbb675 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -165,7 +165,7 @@ Datum geography_area(PG_FUNCTION_ARGS) /* User requests spherical calculation, turn our spheroid into a sphere */ if( ! use_spheroid ) - s.a = s.b = s.radius; + s.a = s.b = s.radius; /* We need the bounding box to get an outside point for area algorithm */ if( ! gbox_from_gserialized(g, &gbox) ) @@ -177,12 +177,15 @@ Datum geography_area(PG_FUNCTION_ARGS) lwgeom = lwgeom_from_gserialized(g); /* Calculate the area */ - area = lwgeom_area_spheroid(lwgeom, gbox, s); + if( use_spheroid ) + area = lwgeom_area_spheroid(lwgeom, gbox, s); + else + area = lwgeom_area_sphere(lwgeom, gbox, s); /* Something went wrong... */ if( area < 0.0 ) { - elog(ERROR, "lwgeom_area_sphere returned area < 0.0"); + elog(ERROR, "lwgeom_area_spher(oid) returned area < 0.0"); PG_RETURN_NULL(); } @@ -193,6 +196,7 @@ Datum geography_area(PG_FUNCTION_ARGS) } + /* ** geography_length_sphere(GSERIALIZED *g) ** returns double length in meters