Geometry Boost List: Re: [geometry] Integrating OGR library with Boost.Geometry

Subject: Re: [geometry] Integrating OGR library with Boost.Geometry
From: Eric MSP Veith (eveith_at_[hidden])
Date: 2015-01-10 12:40:01

Hi Adam, 

thank you for your extensive reply! It already helped me a lot.

I guess that there are actually two topics that are covered by my problem and 
your answer. First, the simple question of how to overwrite bg::distance() 
with OGRGeometry::Distance, and second, the underlying question of why the 
results of bg::distance() and OGRGeometry::Distance actually differ.

Since my data source is a PostGIS database, and PostGIS relies on GEOS just as 
OGR does, I'll go for the low-hangig fruit first before diving deeper into the 
actual problem, if that's OK with you.

I've tried overloading bg::distance(), and tried to overload it for all 
instances of OGRGeometry's derived classes. OGRPoint is derived from 
OGRGeometry, OGRLineString, too, and all the over ones (OGRPolygon, 
OGRLinearRing, ...) the same.

C++11 templates are not (yet) my league, so I've run into problems using 
enable_if<> and is_base_of<>. What I tried is the following:

---%<---

namespace boost {
    namespace geometry {

        template <>
        struct default_distance_result<OGRGeometry, OGRGeometry>
        {
            typedef qreal type;
        };

        template <
                typename Geometry1,
                typename Geometry2,
                typename std::enable_if<
                    std::is_base_of<OGRGeometry, Geometry1>::value
                        and std::is_base_of<
                            OGRGeometry,
                            Geometry2>::value, int>::type = 0>
        inline typename default_distance_result<Geometry1, Geometry2>::value
        distance(Geometry1 const &geometry1, Geometry2 const &geometry2)
        {
            return geometry1.Distance(&geometry2);
        }
    }
}

--->%---

Replacing the specialization's parameters with, e.g. <OGRPoint, OGRLineString> 
does what I want, but it also means I need to create a specialized struct for 
every possible combination of OGRGeometry-derived classes, i.e., 
<OGRLineString, OGRPoint>, <OGRLineString, OGRPolygon>, <OGRPolygon, 
OGRLineString>, and so on.

Do you know a more general way of doing it for all of OGR's geometry classes?

Now to the deeper-lying root cause:

On Thursday 08 January 2015, 19:08:56, Adam Wulkiewicz wrote:
> If you have problems with this function this probably mean that you're
> using geographic CS. Is that right?

Yes, that is right.

> Boost.Geometry supports Point/Point, Point/Linestring, and Point/Box
> distance in cartesian and spherical_equatorial. Couldn't you just switch
> to e.g. spherical equatorial?

I've tried using different coordinate systems, because adapting the 
OGRGeometry classes to Boost's range concept was not that hard. To make things 
easier to explain, I've set up a github.com repository that contains the 
adapter code at <https://github.com/eveith/ogr-boost-adapter>. Perhaps there's 
just an error hidden somewhere in my adapter code.  ;-)

However, I'm presented with different results for OGRGeometry::Distance() and 
bg::distance().

For example, consider the following geometries:

ST_GeomFromText('LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0 4.0)', 
4326),
ST_GeomFromText('LINESTRING(15.0 15.0, 16.0 16.0)', 4326)
ST_GEomFromText('POINT(15.0 15.0)', 4326)

For this, the database gives:

---%<---

db=# SELECT ST_Distance(
    ST_GeomFromText('LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0 4.0)', 
4326),
    ST_GeomFromText('POINT(15.0 15.0)', 4326)); 

   st_distance
-----------------
 15.556349186104

db=# SELECT ST_Distance(
    ST_GeomFromText('LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0 4.0)', 
4326),
    ST_GeomFromText('LINESTRING(15.0 15.0, 16.0 16.0)', 4326)); 
   st_distance
-----------------
 15.556349186104

--->%---

The 15.556... is of course the same result OGR gives me, since both PostGIS 
and OGR use GEOS.

Using cs::cartesian yields 16.9706 on bg::distance for the very same 
geometries. When using cs::geographic, 0.294082 is returned. I'm at loss for 
an explanation why. That was the point where I decided to just use 
OGRGeometry::Distance.  Since I'm not under time pressure (and relatively new 
to anything GIS), I'd like to learn about the backgrounds and hope to come up 
with an actually useable wrapper for OGRGeometry/GEOS classes, using 
Boost.Geometry's facilities the right way --- just resorting to template 
specialization seems like a bad hack to me.   :-)

> Btw, during any operation Boost.Geometry is not doing any projections,
> and probably will never do, at least in 2D. In particular the distance
> in a specific coordinate system (defined by the Point-type of a
> geometry) is calculated (or approximated) using a
> method/algorithm/formula designed for this coordinate system. If the
> user wanted to perform some projection and then some operation for a
> projected geometry, that's ok. But the library will probably always try
> to calculate the "correct" result. This is why for example to calculate
> the distance between the Points the library is using or could use:
> - in cartesian - well-known cartesian distance
> - in spherical - haversine formula or a formula derived from spherical
> law of cosines*
> - in geographic - formulas of Vincenty**, Andoyer-Lambert**,
> Andoyer-Lambert-Thomas* etc.
> 
> * - currently not implemented
> ** - planned for release in 1.58

Perhaps I'm mistaken, but the obvious conclusion seems to be that as long as I 
use the same projection for all geometries I apply predicates to, e.g. force 
WGS84 for all geometries, and choose the right algorithm, I'm on the safe 
side? At least that's how I always understood it.

>From what I could find out, GEOS uses the Discrete Hausdorff Distance 
algorithm? But that doesn't make much sense, given that the Hausdorff Distance 
is the longest distance one has to travel to a Geometry G1 if a point P is 
chosen on another Geometry G0. But again, I lack the knowledge to be sure and 
am most probably mistaken.

> The rest if you have much time to spare :)

Yes, I have.  :-)

> Everything is ok here, see. e.g.
> https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/str
> ategies/distance.hpp line 85.
> Those are Point types from which the algorithm can for instance extract
> coordinate types.

Like boost::model::point?

> In Boost.Geometry operations are divided into 2 parts - an algorithm and
> a strategy. The strategy defines some critical parts of the calculation
> and the algorithms knows how to use the strategies to get the result.

I'm sorry, I cannot wrap my head around it. The Andoyer-Lambert formula, for 
example, would be part of ... the strategy? I.e., the algorithm is 
bg::distance() and the strategy is chosen according to the actual geometries 
and the coordinate system, arriving at, for example, the Pythagoras strategy 
for cs::cartesian or the haversine strategy for cs::geographic?

> So in order to support some combination of geometries there must be an
> algorithm and a strategy adapted to a concept which this algorithm can use.
> 
> E.g. distance(Point, Point) requires a strategy calculating the distance
> between 2 points, obviously. But in the case of distance(Point,
> Linestring) it requires a strategy calculating the distance between a
> Point and a Segment, not Point and Linestring. It's because the strategy
> defines only the most critical part of the calculation and the algorithm
> does the rest.
> 
> So the default distance(Point, Linestring) algorithm knows how to use a
> strategy calculating the distance between a Point and a Segment.
> In order to support your strategy calculating the distance between a
> Point and a Linestring you'd be forced to implement an algorithm that
> knows how to use it.

So, in order to use OGRLineString::Distance(OGRPoint const&), I'd have to 
implement my own bg::distance()? If so, that sounds like "back to square one, 
specialize that template"?

> So you could try to implement Point/Segment strategy instead for your
> coordinate system (geographic?), and e.g. arbitrary geometries, not only
> OGRPoint and OGRLineString.

Yes, geographic. But from that paragraph I fear I got the idea of 
strategy/algorithm wrong. :-(

Thank you very much for the help and the background info. I currently try to 
understand why there are different distance algorithms and when they are 
employed. If you can suggest sources to read up on this, I'll happily consume 
them.

Thanks!

      --- Eric