Geometry Boost List: [geometry] Question about boost::geometry::index::rtree and the nearest predicate

Subject: [geometry] Question about boost::geometry::index::rtree and the nearest predicate
From: Branimir Betov (bbetov_at_[hidden])
Date: 2014-01-22 16:23:49

Hi,

 I posted in boost-users, but looks like this list is better. Sorry about
the double post.

I don't know if I don't understand how the rtree works or there is a
problem with the indexing.

I have an rtree with about 20,000 points and I tried creating two rtrees
with different strategies (linear and rstar) defined as follows:

typedef boost::geometry::model::point<double, 2,
boost::geometry::cs::spherical_equatorial<boost::geometry::degree>> point;
boost::geometry::index::rtree< point, boost::geometry::index::linear<16> >
coastPointsLinear;
boost::geometry::index::rtree< point, boost::geometry::index::rstar<16> >
coastPointsRStar;

When I query for the nearest 3 points, I get different results returned!

coastPointsRStar.query(bgi::nearest(point(latitude, longitude), 3),
std::back_inserter(nearestRStar));
coastPointsLinear.query(bgi::nearest(point(latitude, longitude), 3),
std::back_inserter(nearestLinear));

This seems like a bug. (I attached a text file with the list of the points
I put in the rtree; first two columns are longitude and latitude). The
results also seem to change if I change the maximum number of elements in
nodes from 16 to something else (i.e. 32, 64, etc.)

For what is worth, I also did an experiment with returning points in a box
and it seems like the intersects code returns the same points (works as I
expected) and they seem OK when I plot them in a GIS.

Please help me figure out what I might be doing wrong - here is the whole
test function.

namespace bgi = boost::geometry::index;
typedef boost::geometry::model::point<double, 2,
boost::geometry::cs::spherical_equatorial<boost::geometry::degree>> point;
boost::geometry::index::rtree< point, boost::geometry::index::linear<16> >
coastPointsLinear;
boost::geometry::index::rtree< point, boost::geometry::index::rstar<16> >
coastPointsRStar;
using namespace std;

void Test()
{
    ifstream in("c:\\_Playground\\CoastPoints.txt");
    CSVRow row;

    while( in >> row )
    {
        double lat = atof(row[1].c_str());
        double lon = atof(row[0].c_str());
        coastPointsRStar.insert( point( lat, lon ));
        coastPointsLinear.insert( point( lat, lon ));
    }
    in.close();

    double longitude = 141.029209074;
    double latitude = 38.2836106326;

    vector<point> nearestRStar;
    coastPointsRStar.query(bgi::nearest(point(latitude, longitude), 3),
std::back_inserter(nearestRStar));
    PrintStats(longitude, latitude, nearestRStar, "RStar");

    vector<point> nearestLinear;
    coastPointsLinear.query(bgi::nearest(point(latitude, longitude), 3),
std::back_inserter(nearestLinear));
    PrintStats(longitude, latitude, nearestLinear, "Linear");
}

Thank you