
Geometry : 
Subject: Re: [geometry] Integrating OGR library with Boost.Geometry
From: Eric MSP Veith (eveith_at_[hidden])
Date: 20150110 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 lowhangig 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 OGRGeometryderived 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 deeperlying 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/ogrboostadapter>. 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 Pointtype 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  wellknown cartesian distance
>  in spherical  haversine formula or a formula derived from spherical
> law of cosines*
>  in geographic  formulas of Vincenty**, AndoyerLambert**,
> AndoyerLambertThomas* 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 AndoyerLambert 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
Geometry list run by mateusz at loskot.net