|
|
Geometry : |
Subject: Re: [geometry] Integrating OGR library with Boost.Geometry
From: Adam Wulkiewicz (adam.wulkiewicz_at_[hidden])
Date: 2015-01-10 21:47:01
Hi Eric,
Eric Msp Veith wrote:
> 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);
> }
> }
> }
>
> --->%---
It looks good. Well, your default_distance_result<> specialization
wouldn't be instantiated since it's specialized for base type. But it
probably isn't needed really, the default one should probably work. It
creates the type the distance taking into account coordinate types of
both geometries. The function overload looks ok, if you can use C++11
function default template parameters. In C++98 you could use Boost.Core
(enable_if), Boost.MPL (and_), Boost.TypeTraits (is_base_of) and
implement a return type as enable_if<...>::type. Anyway, the problem is
that it's ambiguous since there are 2 overloads of bg::distance(), both
taking 2 arbitrary types, the original and yours. So this can't be done
this way or at least I don't have any idea.
> 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?
I could think about some hacks but you'd be forced to mess with the
internals of the library so I'd discourage that.
You could implement some convenient macro and then use it like this:
DEFINE_DISTANCE(OGRPoint,OGRPoint)
DEFINE_DISTANCE(OGRPoint,OGRLinestring)
DEFINE_DISTANCE(OGRPoint,OGRBox)
/*...*/
> 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. ;-)
Well, I see now that you've defined your own coordinate system tag. In
this case you should support it by implementation of a set of strategies
and maybe more. In general, it's playing with the internals of the
library which is something that a user shouldn't do. But if you're
interested below is more information about the strategies and how
algorithms use them.
> 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
(Below BG means Boost.Geometry)
In BG, using cartesian CS, the result is 15.556349186104 so the same.
AFAIR if you're using PostGIS like in your example all calculations are
done using cartesian CS.
Instead of ST_GeomFromText you could try to use ST_GeogFromText
(http://postgis.net/docs/ST_GeogFromText.html). But in this case you
should notice that the functionality is limited.
Furthermore, I'm not sure if it's possible to calculate distances (or
perform other operations) in other coordinate systems than cartesian in
GEOS. I'd guess that they probably support SRIDs somehow (store?), but
besides that all calculations are done in cartesian.
Note that GEOS =/= PostGIS, AFAIR in PostGIS there are additional things
implemented.
Some time ago I checked PostGIS code and I found some functions for the
Point/Point distance calculation in geographic CS, maybe there are other
like Point/Segment but I didn't see them. It's possible that they
calculate the closest distance between the Points of geometries, not the
"true"/"correct" one. You could easily check it defining some very long
segment and a Point very close to it somewhere near the middle.
In BG, in spherical_equatorial<degree> the result is 0.26942677292131906
[radians] for a unit sphere, or 1716517.97 [meters] for spherical Earth
with radius 6371000 [meters]. To get the distance in meters you may
multiply the resulting radians by the radius or pass the strategy
explicitly into the distance() function specifying the radius.
FYI, in 1.57, distance(Point, Linestring) or even distance(Point, Point)
doesn't compile for geographic CS. However in develop branch (probably
1.58) the distance(POINT(4, 4), POINT(15, 15)) does compile and the
result is 1712939.35 meters for WGS84 ellipsoid. Currently by default
Andoyer-Lambert formula is used. If Vincenty strategy is passed
explicitly, the result is 1712955.94 meters. You can verify the results
e.g. using this form:
http://www.ga.gov.au/geodesy/datums/vincenty_direct.jsp
> 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.
And Boost.Geometry returns the same distances as the ones mentioned above.
However distance(Linestring, Linestring) doesn't compile for
spherical_equatorial in 1.57 but it compiles in develop so probably
expect the same in 1.58.
>
> 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. :-)
So I guess the answer is already above. First you should probably check
if you really need this workadound.
The right way would be to properly implement missing elements in
Boost.Geometry :)
What bothers me is why your results are slightly different than mine. Do
you get the same for the following program?
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>
namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;
int main()
{
typedef bg::model::point<double, 2, bg::cs::cartesian> point;
typedef bg::model::linestring<point> linestring;
point pt;
linestring ls, ls2;
bg::read_wkt("LINESTRING(0.0 0.0, 1.0 1.0, 2.0 2.0, 3.0 3.0, 4.0
4.0)", ls);
bg::read_wkt("POINT(15.0 15.0)", pt);
bg::read_wkt("LINESTRING(15.0 15.0, 16.0 16.0)", ls2);
std::cout << std::setprecision(24) << bg::distance(pt, ls) <<
std::endl;
std::cout << std::setprecision(24) << bg::distance(ls, ls2) <<
std::endl;
}
>> 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.
Yes, as long as one coordinate system is used and the distances
calculated for various geometries are compatible you should be safe. You
probably shouldn't mix 2 libraries calculating similar things in a
different way, e.g. Boost.Geometry using one formula for Point/Point and
GEOS/GDAL/XXX for Point/Linestring calculating the distance differently,
e.g. using different formula or using cartesian CS, etc. I'm not sure
how this library does it.
I'm guessing that if you used GEOS/GDAL/XXX for all distance
calculations the rtree should work properly, probably... I can't really
guarantee anything.
Furthermore assuming that this library uses cartesian CS (I'm not sure)
it really doesn't have sense to replace Boost.Geometry. So this probably
needs some testing. If you plan to perform them could you please share
the results?
> >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.
It's a different algorithm, it measures how two geometries are similar
to each other.
>> The rest if you have much time to spare :)
> Yes, I have. :-)
Ok :)
>
>> 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?
This, or OGRPoint, every type adapted to the Point concept. In
Boost.Geometry a Point bound with some geometry like Linestring or
Polygon defines the coordinate type, dimension and coordinate system.
This Point can be retrieved from a Geometry by
bg::point_type<Geometry>::type internally using bg::traits::point_type<>
specialized by you when you were adapting your legacy geometries to
Boost.Geometry Concepts.
>
>> 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?
Yes, the Andoyer-Lambert formula is a part of a Strategy. This
particular Strategy defines how a distance between 2 points is
calculated in geographic CS.
The algorithm is the distance(), but internally there are many versions
of it, dispatched for various pairs of input Geometries, by using the
tags of those Geometries. Each dispatched version requires some specific
kind of strategy.
For instance, algorithm dispatched for Point/Point or
MultiPoint/MultiPoint requires Point-Point Strategy. It knows how to use
only this kind of Strategy. For example, Point/Point version calls the
strategy once for a pair of Points. MultiPoint/MultiPoint version could
for instance check all possible pairs of Points and return the smallest
distance found (this is just an example, actually different algorithm is
used since this'd have complexity O(N^2)).
Version for Point/Linestring or Point/Polygon requires Point/Segment
strategy. The former checks all segments in a Linestring, the latter
checks the exterior ring and all interior rings, etc.
For the algorithm it doesn't matter what is the coordinate system, the
required operation is defined in the Strategy.
>
>> 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"?
The short answer is yes, this would be the most simple way of achieving
what you planned to do.
Of course this'd be needed assuming that using cartesian or
spherical_equatorial CS wouldn't be enough.
>> 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.
>
Sorry, I should be more clear before. I already described above what a
strategy and algorithm is. So what an already implemented algorithm
(e.g. distance(Point, Point) or distance(Point, Linestring)) needs is a
Strategy defined for some coordinate system. E.g. to support
distance(PointLike, PointLike) the Point/Point strategy is required
(this is already done in develop for geographic CS as mentioned above).
To support distance(PointLike, LinearOrAreal) we'd need a strategy
calculating the distance between a Point and a Segment on ellipsoid.
This would be required to support all combinations like
Point/Linestring, Point/Polygon, MultiPoint/MultiPolygon etc. to
calculate the "true" distance. Well, actually we'd need more things but
that's another story.
Regards,
Adam
Geometry list run by mateusz at loskot.net