From: Paul Ramsey Date: Wed, 11 Jul 2012 04:05:57 +0000 (+0000) Subject: Fix issue with projecting from the poles, retain the source longitude for more sensib... X-Git-Tag: 2.1.0beta2~813 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e47e3940e97684ba07f62eeb08003d14db85e0d4;p=postgis Fix issue with projecting from the poles, retain the source longitude for more sensible result. 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 6e35fcf91..f472c4e54 100644 --- 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 17a009a08..3bab17cc0 100644 --- 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 1b03e83cf..8c1fcd95a 100644 --- 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);