Using boost::geometry::line_interpolate
with boost::geometry::srs::spheroid
, I'm calculating great circle navigation points along the shortest distance between 2 geographic points. The code below calculates the navigation points for the shortest distance around the great circle. In some rare cases, I need to generate the longer distance that wraps around the globe in the wrong direction. For example, when interpolating between a lon/lat of (20, 20) to (30, 20), there only 10 degrees of difference in the shorter direction and 350 degrees in the other. In some cases I would like the ability to want to interpolate in the longer direction (e.g. 350 deg).
This 2d map shows shows the 10 degree longitude difference in red, and 350 degrees green. I drew the green line by hand to the line is only an approximation. How can I get the points for this green line?
This code is based on the example from boost.org, line_interpolate_4_with_strategy
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
int main()
{
typedef boost::geometry::model::d2::point_xy<double, boost::geometry::cs::geographic<boost::geometry::degree> > Point_Type;
using Segment_Type = boost::geometry::model::segment<Point_Type>;
using Multipoint_Type = boost::geometry::model::multi_point<Point_Type>;
boost::geometry::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
boost::geometry::strategy::line_interpolate::geographic<boost::geometry::strategy::vincenty> str(spheroid);
Segment_Type const start_end_points { {20, 20}, {30, 20} }; // lon/lat, interpolate between these two points
double distance { 50000 }; // plot a point ever 50km
Multipoint_Type mp;
boost::geometry::line_interpolate(start_end_points, distance, mp, str);
std::cout << "on segment : " << wkt(mp) << "\n";
return 0;
}
Note that line_interpolate
interpolates points on a linestring where a segment between two points follows a geodesic.
Therefore, one workaround could be to create an antipodal point to the centroid of the original segment and create a linestring that follows the requested path. Then call line_interpolate
with this linestring. The following code could do the trick.
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
int main()
{
namespace bg = boost::geometry;
using Point_Type = bg::model::d2::point_xy<double, bg::cs::geographic<bg::degree>>;
using Segment_Type = boost::geometry::model::segment<Point_Type>;
using Linstring_Type = bg::model::linestring<Point_Type>;
using Multipoint_Type = bg::model::multi_point<Point_Type>;
bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
bg::strategy::line_interpolate::geographic<bg::strategy::vincenty> str(spheroid);
Segment_Type const start_end_points { {20, 20}, {30, 20} };
Point_Type centroid;
bg::centroid(start_end_points, centroid);
Point_Type antipodal_centroid;
bg::set<0>(antipodal_centroid, bg::get<0>(centroid) + 180);
bg::set<1>(antipodal_centroid, bg::get<1>(centroid) * -1);
Linstring_Type line;
line.push_back(start_end_points.first);
line.push_back(antipodal_centroid);
line.push_back(start_end_points.second);
double distance { 50000 }; // plot a point ever 50km
Multipoint_Type mp;
bg::line_interpolate(line, distance, mp, str);
std::cout << "on segment : " << wkt(mp) << "\n";
return 0;
}
Note that since the spheroid you are constructing is non-spherical then there is no great circle (apart from equator and meridians) and the geodesic segment is not closed but looks like this. Therefore you will notice that the last interpolated point will be different from the segment's endpoint.