Opened 6 years ago

Last modified 6 years ago

#12314 new Bugs

boost::geometry::within for a point in a geographic coordinate system fails

Reported by: oleary.simon@… Owned by: Barend Gehrels
Milestone: To Be Determined Component: geometry
Version: Boost 1.61.0 Severity: Showstopper
Keywords: geometry within geographic coordinate system cs Cc: oleary.simon@…

Description

I am using the geometry part of the library for a GIS application. The idea is to find points within defined regions (a rectangle in this case). This works sometimes but not always. Here is an example where the point should be within the rectangle but isn't...

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>
#include <boost/geometry/geometries/register/point.hpp>
#include <iostream>
#include <boost/geometry/io/wkt/wkt.hpp>

class wxPoint
{
public :
    double  getx() const
    {
        return m_x;
    }
    double gety() const
    {
        return m_y;
    }
    void setx( double in)
    {
        m_x = in;
    }
    void sety(double in)
    {
        m_y = in;
    }
private:
    double m_x;
    double m_y;
};

BOOST_GEOMETRY_REGISTER_POINT_2D_GET_SET(
    wxPoint,
    double,
    boost::geometry::cs::geographic<boost::geometry::degree>,
    wxPoint::getx,
    wxPoint::gety,
    wxPoint::setx,
    wxPoint::sety )

int main()
{
    boost::geometry::model::polygon< wxPoint > poly;

    boost::geometry::read_wkt( "POLYGON((0 89, 180 89, 180 0, 0 0, 0 89 ))", poly );  

    wxPoint point;
    point.setx( 150 );
    point.sety( 88 );

    bool within = boost::geometry::within( point, poly );


    return 0;
}

Change History (1)

comment:1 by awulkiew, 6 years ago

In non-cartesian coordinate systems the edges of Polygons are also shortest lines between points so e.g. the edge: (0 89, 180 89) intersects the north pole. In the case of the edge (180 0, 0 0) I'm not sure what the algorithm assumes but this edge is invalid. The reason for that is that there is infinite number of possible edges connecting these points in spherical and 2 possible edges in geographic (going through north or south pole). Either way I presume this is not what you had in mind.

You could fix the invalidity by adding additional point between the endpoints of the invalid edge, so: POLYGON((0 89, 180 89, 180 0, 90 0, 0 0, 0 89 )) but this is not what you wanted to do either.

AFAIU what you wanted to do is to create a "rectangular" region (in spherical/geographic having edges being parallels and meridians) and that you should use your own rectangle representation adapted to Boost.Geometry Box concept or bg::model::box for this:

boost::geometry::model::box< wxPoint > box;
boost::geometry::read_wkt("BOX(0 0, 180 89)", box);

or

boost::geometry::model::box< wxPoint > box(pt1, pt2);

So I don't think this is a bug.

Does the above explanation make sense and solution work for you?

Last edited 6 years ago by awulkiew (previous) (diff)
Note: See TracTickets for help on using tickets.