Fix issue with projecting from the poles, retain the source longitude for more sensib...

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

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 6e35fcf9144f863936cb51e4a2aed52ff5f45a7c..f472c4e545b16229bf15b29f424f1fb7402ccaee 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1053,19 +1053,27 @@ int sphere_project(const GEOGRAPHIC_POINT *r, double distance, double azimuth, G
        double lat2, lon2;
 
        lat2 = asin(sin(lat1)*cos(d) + cos(lat1)*sin(d)*cos(azimuth));
-       lon2 = lon1 + atan2(sin(azimuth)*sin(d)*cos(lat1), cos(d)-sin(lat1)*sin(lat2));
-       n->lat = lat2;
-       n->lon = lon2;
 
+       /* If we're going straight up or straight down, we don't need to calculate the longitude */
+       if ( FP_EQUALS(azimuth, M_PI) || FP_EQUALS(azimuth, 0.0) )
+       {
+               lon2 = r->lon;
+       }
+       else
+       {
+               lon2 = lon1 + atan2(sin(azimuth)*sin(d)*cos(lat1), cos(d)-sin(lat1)*sin(lat2));
+       }
+       
        if ( isnan(lat2) || isnan(lon2) )
                return LW_FAILURE;
 
+       n->lat = lat2;
+       n->lon = lon2;
+
        return LW_SUCCESS;
 }
 
 
-
-
 int edge_calculate_gbox_slow(const GEOGRAPHIC_EDGE *e, GBOX *gbox)
 {
        int steps = 1000000;

diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c
index 17a009a0850a89cb64586e768a7961919583f860..3bab17cc0d838ec8a3feeda44af8e98d2492c7d5 100644 (file)
--- a/liblwgeom/lwgeodetic_tree.c
+++ b/liblwgeom/lwgeodetic_tree.c
@@ -51,6 +51,8 @@ circ_node_leaf_new(const POINTARRAY* pa, int i)
        p2 = (POINT2D*)getPoint_internal(pa, i+1);
        geographic_point_init(p1->x, p1->y, &g1);
        geographic_point_init(p2->x, p2->y, &g2);
+
+       LWDEBUGF(3,"edge #%d (%g %g, %g %g)", i, p1->x, p1->y, p2->x, p2->y);
        
        diameter = sphere_distance(&g1, &g2);
 
@@ -72,6 +74,8 @@ circ_node_leaf_new(const POINTARRAY* pa, int i)
        node->center = gc;
        node->radius = diameter / 2.0;
 
+       LWDEBUGF(3,"edge #%d CENTER(%g %g) RADIUS=%g", i, gc.lon, gc.lat, node->radius);
+
        /* Leaf has no children */
        node->num_nodes = 0;
        node->nodes = NULL;
@@ -181,7 +185,7 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
        double offset1, dist, D, r1, ri;
        int i;
 
-       LWDEBUGF(3, "called with %d nodes", num_nodes);
+       LWDEBUGF(3, "called with %d nodes --", num_nodes);
 
        /* Can't do anything w/ empty input */
        if ( num_nodes < 1 )
@@ -202,33 +206,42 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
 
                LWDEBUGF(3, "distance between new (%g %g) and %i (%g %g) is %g", c1.lon, c1.lat, i, c[i]->center.lon, c[i]->center.lat, dist);
                
-               if ( FP_EQUALS(dist, 0) || dist <= fabs(r1 - ri) )
+               if ( FP_EQUALS(dist, 0) )
+               {
+                       LWDEBUG(3, "  distance between centers is zero");
+                       new_radius = r1 + 2*dist;
+                       new_center = c1;
+               }
+               else if ( dist < fabs(r1 - ri) )
                {
                        /* new contains next */
                        if ( r1 > ri )
                        {
+                               LWDEBUG(3, "  c1 contains ci");
                                new_center = c1;
                                new_radius = r1;
                        }
                        /* next contains new */
                        else
                        {
+                               LWDEBUG(3, "  ci contains c1");
                                new_center = c[i]->center;
                                new_radius = ri;
                        }
                }
                else
                {       
+                       LWDEBUG(3, "  calculating new center");
                        /* New circle diameter */
                        D = dist + r1 + ri;
-                       LWDEBUGF(4,"D is %g", D);
+                       LWDEBUGF(3,"    D is %g", D);
                        
                        /* New radius */
                        new_radius = D / 2.0;
                        
                        /* Distance from cn1 center to the new center */
                        offset1 = ri + (D - (2.0*r1 + 2.0*ri)) / 2.0;
-                       LWDEBUGF(4,"offset1 is %g", offset1);
+                       LWDEBUGF(3,"    offset1 is %g", offset1);
                        
                        /* Sometimes the sphere_direction function fails... this causes the center calculation */
                        /* to fail too. In that case, we're going to fall back ot a cartesian calculation, which */
@@ -240,7 +253,7 @@ circ_node_internal_new(CIRC_NODE** c, int num_nodes)
                                new_radius *= 1.1;
                        }
                }
-               LWDEBUGF(3, "new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);      
+               LWDEBUGF(3, " new center is (%g %g) new radius is %g", new_center.lon, new_center.lat, new_radius);     
        }
        
        node = lwalloc(sizeof(CIRC_NODE));
@@ -292,7 +305,6 @@ circ_tree_new(const POINTARRAY* pa)
        for ( i = 0; i < num_edges; i++ )
        {
                node = circ_node_leaf_new(pa, i);
-               LWDEBUGF(4,"making new leaf node %d", i);
                if ( node ) /* Not zero length? */
                        nodes[j++] = node;
        }
@@ -339,7 +351,7 @@ circ_nodes_merge(CIRC_NODE** nodes, int num_nodes)
                        CIRC_NODE **inodes = lwalloc(sizeof(CIRC_NODE*)*node_size);
                        inodes[0] = nodes[node_size*j];
                        inodes[1] = nodes[node_size*j + 1];
-                       LWDEBUGF(5,"merging nodes %d and %d", node_size*j, node_size*j + 1);                    
+                       LWDEBUGF(3,"merging nodes %d and %d", node_size*j, node_size*j + 1);                    
                        nodes[j++] = circ_node_internal_new(inodes, node_size);
                }
                /* Odd number of children, just copy the last node up a level */

diff --git a/postgis/geography_measurement_trees.c b/postgis/geography_measurement_trees.c
index 1b03e83cfd94bd656b55b1f07f4eac2182479ad6..8c1fcd95a640edeac3c3dacbb49ee233e17c8ea3 100644 (file)
--- a/postgis/geography_measurement_trees.c
+++ b/postgis/geography_measurement_trees.c
@@ -94,8 +94,8 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
                /* Need a gbox to calculate an outside point */
                if ( LW_FAILURE == gserialized_get_gbox_p(g1, &gbox1) )
                {
-                       POSTGIS_DEBUG(3, "unable to read gbox from gserialized, calculating from scratch");
                        LWGEOM* lwgeom1 = lwgeom_from_gserialized(g1);
+                       POSTGIS_DEBUG(3, "unable to read gbox from gserialized, calculating from scratch");
                        lwgeom_calculate_gbox_geodetic(lwgeom1, &gbox1);
                        lwgeom_free(lwgeom1);
                }
@@ -127,6 +127,7 @@ CircTreePIP(const CIRC_NODE* tree1, const GSERIALIZED* g1, const LWGEOM* lwgeom2
                        pt2d_inside.y = in_point.y;
                        /* Calculate a definitive outside point */
                        gbox_pt_outside(&gbox1, &pt2d_outside);
+                       POSTGIS_DEBUGF(3, "p2d_inside=POINT(%g %g) p2d_outside=POINT(%g %g)", pt2d_inside.x, pt2d_inside.y, pt2d_outside.x, pt2d_outside.y);
                        /* Test the candidate point for strict containment */
                        POSTGIS_DEBUG(3, "calling circ_tree_contains_point for PiP test");
                        return circ_tree_contains_point(tree1, &pt2d_inside, &pt2d_outside, NULL);