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