self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-23 16:42:35 +00:00
Code Issues Packages Projects Releases Wiki Activity

Improve support for ST_Area(geography) over dateline and poles (#2006, #2039)

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@10407 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Paul Ramsey 2012-10-11 22:48:11 +00:00
parent c6e5b63756
commit e9bf8eeb82
5 changed files with 121 additions and 84 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
NEWS
View file

@ -45,6 +45,7 @@ PostGIS 2.1.0
- #1978, wrong answer when calculating length of a closed circular arc (circle)
- #1780, support ST_GeoHash for geography
- #2021, Added multi-band support to ST_Union(raster, ...) aggregate function
- #2006, better support of ST_Area(geography) over poles and dateline
* Fixes *

23
liblwgeom/cunit/cu_geodetic.c
View file

@ -990,7 +990,7 @@ static void test_spheroid_area(void)
a1 = lwgeom_area_sphere(lwg, &s);
a2 = lwgeom_area_spheroid(lwg, &s);
//printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
CU_ASSERT_DOUBLE_EQUAL(a1, 89.7211470368, 0.0001); /* sphere */
CU_ASSERT_DOUBLE_EQUAL(a1, 89.7127703297, 0.0001); /* sphere */
CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.0001); /* spheroid */
lwgeom_free(lwg);
@ -1009,7 +1009,7 @@ static void test_spheroid_area(void)
lwgeom_calculate_gbox_geodetic(lwg, &gbox);
a1 = lwgeom_area_sphere(lwg, &s);
a2 = lwgeom_area_spheroid(lwg, &s);
//printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2);
CU_ASSERT_DOUBLE_EQUAL(a1, 12360265021.1, 10.0); /* sphere */
CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */
lwgeom_free(lwg);
@ -1157,6 +1157,24 @@ static void test_lwgeom_segmentize_sphere(void)
return;
}
static void test_lwgeom_area_sphere(void)
{
LWGEOM *lwg;
double area;
SPHEROID s;
/* Init to WGS84 */
spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
/* Simple case */
lwg = lwgeom_from_wkt("POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", LW_PARSER_CHECK_NONE);
area = lwgeom_area_sphere(lwg, &s);
CU_ASSERT_DOUBLE_EQUAL(area, 12360265021.3561, 0.01);
lwgeom_free(lwg);
return;
}
/*
** Used by test harness to register the tests in this file.
*/
@ -1164,6 +1182,7 @@ CU_TestInfo geodetic_tests[] =
{
PG_TEST(test_sphere_direction),
PG_TEST(test_sphere_project),
PG_TEST(test_lwgeom_area_sphere),
PG_TEST(test_signum),
PG_TEST(test_gbox_from_spherical_coordinates),
PG_TEST(test_gserialized_get_gbox_geocentric),

178
liblwgeom/lwgeodetic.c
View file

@ -564,12 +564,12 @@ int crosses_dateline(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e)
}
/**
* Returns true if the point p is on the great circle plane.
* Forms the scalar triple product of A,B,p and if the volume of the
* resulting parallelepiped is near zero the point p is on the
* great circle plane.
* Returns -1 if the point is to the left of the plane formed
* by the edge, 1 if the point is to the right, and 0 if the
* point is on the plane.
*/
int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
static int
edge_point_side(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
{
POINT3D normal, pt;
double w;
@ -583,9 +583,79 @@ int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
if ( FP_IS_ZERO(w) )
{
LWDEBUG(4, "point is on plane (dot product is zero)");
return LW_TRUE;
return 0;
}
LWDEBUG(4, "point is not on plane");
if ( w < 0 )
return -1;
else
return 1;
}
/**
* Returns the angle in radians at point B of the triangle formed by A-B-C
*/
static double
sphere_angle(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
{
POINT3D normal1, normal2;
robust_cross_product(b, a, &normal1);
robust_cross_product(b, c, &normal2);
normalize(&normal1);
normalize(&normal2);
return sphere_distance_cartesian(&normal1, &normal2);
}
/**
* Computes the spherical area of a triangle. If C is to the left of A/B,
* the area is negative. If C is to the right of A/B, the area is positive.
*
* @param a The first triangle vertex.
* @param b The second triangle vertex.
* @param c The last triangle vertex.
* @return the signed area in radians.
*/
static double
sphere_signed_area(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, const GEOGRAPHIC_POINT *c)
{
double angle_a, angle_b, angle_c;
double area_radians = 0.0;
int side;
GEOGRAPHIC_EDGE e;
angle_a = sphere_angle(b,a,c);
angle_b = sphere_angle(a,b,c);
angle_c = sphere_angle(b,c,a);
area_radians = angle_a + angle_b + angle_c - M_PI;
/* What's the direction of the B/C edge? */
e.start = *a;
e.end = *b;
side = edge_point_side(&e, c);
/* Co-linear points implies no area */
if ( side == 0 )
return 0.0;
/* Add the sign to the area */
return side * area_radians;
}
/**
* Returns true if the point p is on the great circle plane.
* Forms the scalar triple product of A,B,p and if the volume of the
* resulting parallelepiped is near zero the point p is on the
* great circle plane.
*/
int edge_point_on_plane(const GEOGRAPHIC_EDGE *e, const GEOGRAPHIC_POINT *p)
{
int side = edge_point_side(e, p);
if ( side == 0 )
return LW_TRUE;
return LW_FALSE;
}
@ -806,7 +876,6 @@ double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, do
return heading;
}
/**
* Computes the spherical excess of a spherical triangle defined by
* the three vectices A, B, C. Computes on the unit sphere (i.e., divides
@ -832,7 +901,6 @@ static double sphere_excess(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b
}
/**
* Returns true if the point p is on the minor edge defined by the
* end points of e.
@ -1749,72 +1817,31 @@ lwgeom_segmentize_sphere(const LWGEOM *lwg_in, double max_seg_length)
* Returns the area of the ring (ring must be closed) in square radians (surface of
* the sphere is 4*PI).
*/
double ptarray_area_sphere(const POINTARRAY *pa, const POINT2D *pt_outside)
double
ptarray_area_sphere(const POINTARRAY *pa)
{
GEOGRAPHIC_POINT a, b, c;
POINT2D p;
int i;
const POINT2D *p;
GEOGRAPHIC_POINT a, b, c;
double area = 0.0;
/* Return zero on non-sensical inputs */
/* Return zero on nonsensical inputs */
if ( ! pa || pa->npoints < 4 )
return 0.0;
geographic_point_init(pt_outside->x, pt_outside->y, &c);
/* Initialize first point */
getPoint2d_p(pa, 0, &p);
geographic_point_init(p.x, p.y, &a);
for ( i = 1; i < pa->npoints; i++ )
p = getPoint2d_cp(pa, 0);
geographic_point_init(p->x, p->y, &a);
p = getPoint2d_cp(pa, 1);
geographic_point_init(p->x, p->y, &b);
for ( i = 2; i < pa->npoints-1; i++ )
{
double excess = 0.0;
getPoint2d_p(pa, i, &p);
geographic_point_init(p.x, p.y, &b);
if ( crosses_dateline(&a, &b) )
{
GEOGRAPHIC_POINT a1 = a, b1 = b, c1 = c;
double shift;
if ( a.lon > 0.0 )
shift = (M_PI - a.lon) + 0.088; /* About 5deg more */
else
shift = (M_PI - b.lon) + 0.088; /* About 5deg more */
LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g) c1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon, c1.lat, c1.lon);
point_shift(&a1, shift);
point_shift(&b1, shift);
point_shift(&c1, shift);
LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g) c1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon, c1.lat, c1.lon);
excess = sphere_excess(&a1, &b1, &c1);
}
else if ( crosses_dateline(&a, &c) )
{
GEOGRAPHIC_POINT a1 = a, b1 = b, c1 = c;
double shift;
if ( a.lon > 0.0 )
shift = (M_PI - a.lon) + 0.088; /* About 5deg more */
else
shift = (M_PI - c.lon) + 0.088; /* About 5deg more */
point_shift(&a1, shift);
point_shift(&b1, shift);
point_shift(&c1, shift);
excess = sphere_excess(&a1, &b1, &c1);
}
else
{
excess = sphere_excess(&a, &b, &c);
}
area += excess;
/* B gets incremented in the next loop, so we save the value here */
a = b;
p = getPoint2d_cp(pa, i);
geographic_point_init(p->x, p->y, &c);
area += sphere_signed_area(&a, &b, &c);
b = c;
}
return fabs(area);
}
@ -2106,7 +2133,6 @@ static double ptarray_distance_spheroid(const POINTARRAY *pa1, const POINTARRAY
double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
{
int type;
POINT2D pt_outside;
double radius2 = spheroid->radius * spheroid->radius;
GBOX gbox;
gbox.flags = 0;
@ -2124,16 +2150,6 @@ double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
if ( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
return 0.0;
/* Make sure we have boxes */
if ( lwgeom->bbox )
gbox = *(lwgeom->bbox);
else
lwgeom_calculate_gbox_geodetic(lwgeom, &gbox);
gbox_pt_outside(&gbox, &pt_outside);
LWDEBUGF(2, "pt_outside = POINT(%.20g %.20g)", pt_outside.x, pt_outside.y);
/* Actually calculate area */
if ( type == POLYGONTYPE )
{
@ -2146,12 +2162,12 @@ double lwgeom_area_sphere(const LWGEOM *lwgeom, const SPHEROID *spheroid)
return 0.0;
/* First, the area of the outer ring */
area += radius2 * ptarray_area_sphere(poly->rings[0], &pt_outside);
area += radius2 * ptarray_area_sphere(poly->rings[0]);
/* Subtract areas of inner rings */
for ( i = 1; i < poly->nrings; i++ )
{
area -= radius2 * ptarray_area_sphere(poly->rings[i], &pt_outside);
area -= radius2 * ptarray_area_sphere(poly->rings[i]);
}
return area;
}

2
regress/tickets.sql
View file

@ -272,7 +272,7 @@ SELECT ST_AsText(the_geog) as the_pt,
ST_Area(ST_Buffer(the_geog,10)) As the_area,
ST_Area(geography(ST_Transform(ST_Buffer(ST_Transform(geometry(the_geog),utm_srid),10),4326))) As geog_utm_area
FROM utm_dots
WHERE ST_Area(ST_Buffer(the_geog,10)) NOT between 307 and 314
WHERE ST_Area(ST_Buffer(the_geog,10)) NOT between 307 and 315
LIMIT 10;
SELECT '#304.a', Count(*) FROM utm_dots WHERE ST_DWithin(the_geog, 'POINT(0 0)'::geography, 3000000);

1
regress/tickets_expected
View file

@ -64,6 +64,7 @@ NOTICE: No points or linestrings in input array
#277|<gml:Point><gml:coordinates>1,1e+308</gml:coordinates></gml:Point>
#299|2
#304
ERROR: GetProj4StringSPI: Cannot find SRID (32704) in spatial_ref_sys
#304.a|21
#304.b|1
#408|Too few points in geometry component[