Geometry Boost List: [ggl] Move a point_ll by a distance

Subject: [ggl] Move a point_ll by a distance
From: Barend Gehrels (Barend.Gehrels)
Date: 2009-10-21 05:57:57

Hi Richard,

Thanks for your description.

Sorry, forgot that our SVN needs accounts. I can give you an account if 
you wish. However, we're in preparation for the review, many things are 
changing, so you might wait for a few weeks. We're doing our best to 
have it in Boost :-)

The source code fits in the example 02_point_ll_example, I'll give the 
relevant pieces here:

// Formula to get the course (direction) between two points.
// This might be a GGL-function in the future.
template <typename P1, typename P2>
inline double *get_course*(P1 const& p1, P2 const& p2)
{
    double const& lat1 = ggl::get_as_radian<1>(p1);
    double const& lon1 = ggl::get_as_radian<0>(p1);
    double const& lat2 = ggl::get_as_radian<1>(p2);
    double const& lon2 = ggl::get_as_radian<0>(p2);
    // http://williams.best.vwh.net/avform.htm#Crs
    return atan2(sin(lon1-lon2)*cos(lat2),
       cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon1-lon2));
}

// Formula to calculate the point at a distance/angle from another point
// This might be a GGL-function in the future.
template <typename P1, typename P2>
inline void *point_at_distance*(P1 const& p1,
        double distance, double tc, double radius,
        P2& p2)
{
    double earth_perimeter = radius * ggl::math::two_pi;
    double d = (distance / earth_perimeter) * ggl::math::two_pi;
    double const& lat1 = ggl::get_as_radian<1>(p1);
    double const& lon1 = ggl::get_as_radian<0>(p1);

    // http://williams.best.vwh.net/avform.htm#LL
    double lat = asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc));
    double dlon = atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat));
    double lon = lon1 - dlon;

    ggl::set_from_radian<1>(p2, lat);
    ggl::set_from_radian<0>(p2, lon);
}

And this can be called like this (the direction is here calculated, but 
you will specified it directly):

    // Other way round: have Amsterdam and go 430 km to the south (i.e. 
first calculate direction)
    double tc = get_course(amsterdam, paris);
    std::cout << "Course: " << (tc * ggl::math::r2d) << std::endl;

    double const average_earth_radius = 6372795.0;
    point_ll_deg paris_calculated;
    point_at_distance(amsterdam, 430 * 1000.0, tc, average_earth_radius, 
paris_calculated);
    std::cout << "Paris calculated (degree): " << 
ggl::wkt(paris_calculated) << std::endl;

Regards, Barend

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/ggl/attachments/20091021/f1a7470c/attachment.html