#2534, st_distance is returning incorrect results for large geographies

git-svn-id: http://svn.osgeo.org/postgis/trunk@12228 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/NEWS b/NEWS
index e97a1244cc870075cab071be62a29b26165a8053..7ffde6541482153c6adba05cf26f75ec75684a66 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -43,6 +43,7 @@ PostGIS 2.2.0
   - #2512, Support for foreign tables and materialized views in
            raster_columns and raster_overviews
   - #2532, Add missing raster/geometry commutator operators
+  - #2534, st_distance is returning incorrect results for large geographies
   - #2540, Change GUC name for GDAL_DATA to postgis.gdal_datapath
   - #2543, invalid join selectivity error from simple query
   - #2552, Fix handling of NULL raster in ST_AsTIFF, ST_AsPNG

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 2a22b1ff4a67f18c28205732d47179d266c98ebf..f9b41d54260537ee30e2bf5de2210fb753c3362f 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -674,6 +674,31 @@ static void test_edge_intersects(void)
        line2pts("LINESTRING(90.0 80.0, -90.0 90.0)", &B1, &B2);
        rv = edge_intersects(&A1, &A2, &B1, &B2);
        CU_ASSERT(rv == (PIR_INTERSECTS|PIR_B_TOUCH_LEFT|PIR_A_TOUCH_RIGHT) );
+
+       /* Antipodal straddles. Great circles cross but at opposite */
+       /* sides of the globe */
+       /* #2534 */
+       /* http://www.gcmap.com/mapui?P=60N+90E-20S+90E%0D%0A0N+0E-90.04868865037885W+57.44011727050777S%0D%0A&MS=wls&DU=mi */
+       line2pts("LINESTRING(90.0 60.0, 90.0 -20.0)", &A1, &A2);
+       line2pts("LINESTRING(0.0 0.0, -90.04868865037885 -57.44011727050777)", &B1, &B2);
+       rv = edge_intersects(&A1, &A2, &B1, &B2);
+       CU_ASSERT(rv == 0);
+
+       line2pts("LINESTRING(-5 0, 5 0)", &A1, &A2);
+       line2pts("LINESTRING(179 -5, 179 5)", &B1, &B2);
+       rv = edge_intersects(&A1, &A2, &B1, &B2);
+       CU_ASSERT(rv == 0);
+
+       line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
+       line2pts("LINESTRING(65 0, -105 0)", &B1, &B2);
+       rv = edge_intersects(&A1, &A2, &B1, &B2);
+       CU_ASSERT(rv == 0);
+       
+       line2pts("LINESTRING(175 -85, 175 85)", &A1, &A2);
+       line2pts("LINESTRING(45 0, -125 0)", &B1, &B2);
+       rv = edge_intersects(&A1, &A2, &B1, &B2);
+       CU_ASSERT(rv == 0);
+       
 }
 
 static void test_edge_distance_to_point(void)

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index a542a870e28e24a50eadd511f87395ddc4ec638a..bc8bf05948435cd7b2a2eb3ee493009a6d245e6d 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -3078,7 +3078,7 @@ dot_product_side(const POINT3D *p, const POINT3D *q)
 int 
 edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const POINT3D *B2)
 {
-       POINT3D AN, BN;  /* Normals to plane A and plane B */
+       POINT3D AN, BN, VN;  /* Normals to plane A and plane B */
        double ab_dot;
        int a1_side, a2_side, b1_side, b2_side;
        int rv = PIR_NO_INTERACT;
@@ -3126,8 +3126,21 @@ edge_intersects(const POINT3D *A1, const POINT3D *A2, const POINT3D *B1, const P
        if ( a1_side != a2_side && (a1_side + a2_side) == 0 &&
             b1_side != b2_side && (b1_side + b2_side) == 0 )
        {
-               /* Mid-point intersection! */
-               return PIR_INTERSECTS;
+               /* Have to check if intersection point is inside both arcs */
+               unit_normal(&AN, &BN, &VN);
+               if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
+               {
+                       return PIR_INTERSECTS;
+               }
+
+               /* Have to check if intersection point is inside both arcs */
+               vector_scale(&VN, -1);
+               if ( point_in_cone(A1, A2, &VN) && point_in_cone(B1, B2, &VN) )
+               {
+                       return PIR_INTERSECTS;
+               }
+               
+               return PIR_NO_INTERACT;
        }
 
        /* The rest are all intersects variants... */