1/*  Part of SWI-Prolog
    2
    3    Author:        Willem Robert van Hage
    4    E-mail:        W.R.van.Hage@vu.nl
    5    WWW:           http://www.few.vu.nl/~wrvhage
    6    Copyright (c)  2009-2012, Vrije Universiteit Amsterdam
    7    All rights reserved.
    8
    9    Redistribution and use in source and binary forms, with or without
   10    modification, are permitted provided that the following conditions
   11    are met:
   12
   13    1. Redistributions of source code must retain the above copyright
   14       notice, this list of conditions and the following disclaimer.
   15
   16    2. Redistributions in binary form must reproduce the above copyright
   17       notice, this list of conditions and the following disclaimer in
   18       the documentation and/or other materials provided with the
   19       distribution.
   20
   21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   22    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
   24    FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
   25    COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
   26    INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
   27    BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
   28    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   29    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
   30    LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
   31    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
   32    POSSIBILITY OF SUCH DAMAGE.
   33*/
   34
   35:- module(space,
   36          [
   37           set_space/1,               % +Option
   38           set_space/2,               % +IndexName, +Option
   39           space_setting/1,           % ?Option
   40
   41           space_index_all/1,         % +IndexNames
   42           space_index_all/0,         % uses default index
   43
   44           space_assert/3,            % +URI, +Shape, +IndexName
   45           space_assert/2,            % uses default index
   46           space_retract/3,           % +URI, +Shape, +IndexName
   47           space_retract/2,           % uses default index
   48           space_index/1,             % +IndexName
   49           space_index/0,             % uses default index
   50           space_clear/1,             % +IndexName
   51           space_clear/0,             % uses default index
   52
   53           space_bulkload/2,          % +CandidatePred, +IndexName
   54           space_bulkload/1,          % +CandidatePred (uses default index and 'user' module)
   55           space_bulkload/0,          % also defaults to uri_shape
   56
   57           space_contains/3,          % +Query, -URI, +IndexName
   58           space_contains/2,          % uses default index
   59           space_intersects/3,        % +Query, -URI, +IndexName
   60           space_intersects/2,        % uses default index
   61           space_nearest/3,           % +Query, -URI, +IndexName
   62           space_nearest/2,           % uses default index
   63           space_within_range/4,      % +Query, -URI, +WithinRange, +IndexName
   64           space_within_range/3,      % uses default index
   65           space_nearest_bounded/4,   % +Query, -URI, +WithinRange, +IndexName
   66           space_nearest_bounded/3,   % uses default index
   67
   68           shape/1,                   % +Shape
   69           uri_shape/2,               % ?URI, ?Shape
   70           uri_shape/3,               % ?URI, ?Shape, +Source
   71           uri_shape_graph/3,         % +Graph, ?URI, ?Shape
   72
   73           space_distance/3,              % +Feature1, +Feature2, -Distance
   74           space_distance/4,              % +Feature1, +Feature2, -Distance, +IndexName
   75           space_distance_pythagorean/3,  % +Feature1, +Feature2, -Distance
   76           space_distance_greatcircle/4,  % +Feature1, +Feature2, -Distance, +Unit
   77           space_distance_greatcircle/3,  % +Feature1, +Feature2, -Distance (nm)
   78
   79           space_bearing/3            % +Point, +Point, -Heading (degrees)
   80
   81          ]).   82
   83:- use_module(wkt).   84:- use_module(kml).   85:- use_module(georss). % also GML support
   86:- use_module(wgs84).   87:- use_module(freebase).   88:- use_module(dbpedia).   89:- use_module(library(error)).   90:- use_module(library(shlib)).   91
   92:- rdf_meta(space_nearest(t,?)).   93:- rdf_meta(space_nearest(t,?,?)).   94:- rdf_meta(space_intersects(t,?)).   95:- rdf_meta(space_intersects(t,?,?)).   96:- rdf_meta(space_contains(t,?)).   97:- rdf_meta(space_contains(t,?,?)).   98:- rdf_meta(space_nearest_bounded(t,?,?)).   99:- rdf_meta(space_nearest_bounded(t,?,?,?)).  100:- rdf_meta(space_within_range(t,?,?)).  101:- rdf_meta(space_within_range(t,?,?,?)).  102:- rdf_meta(space_assert(r,?)).  103:- rdf_meta(space_assert(r,?,?)).  104:- rdf_meta(space_retract(r,?)).  105:- rdf_meta(space_retract(r,?,?)).  106:- rdf_meta(uri_shape(r,?)).  107:- rdf_meta(uri_shape(r,?,?)).  108
  109:- dynamic space_queue/4.  110:- dynamic shape/1. % allows you to adapt space_index_all.
  111:- dynamic uri_shape/2. % allows you to adapt space_index_all.
  112:- dynamic uri_shape/3. % allows you to adapt space_index_all.
  113
  114:- use_foreign_library(foreign(space)).

Core spatial database

*/

 configuration
  124:- dynamic space_setting_aux/1.  125space_setting_aux(rtree_default_index(space_index)).
  126
  127% foreign language predicate:
  128% set_space(IndexName,Option)
 set_space(+Option) is det
 set_space(+IndexName, +Option) is det
This predicate can be used to change the options of a spatial index (or de default index for set_space/1). Some options, like rtree_storage(S) where S is disk or memory only have effect after clearing or bulkloading. Others, take effect immediately on a running index. More documentation will be provided in the near future.
  139set_space(Option) :-
  140    space_setting(rtree_default_index(I)),
  141    rtree_set_space(I, Option).
  142
  143set_space(I, Option) :-
  144    rtree_set_space(I, Option).
  145
  146/*
  147set_space(Option) :-
  148        functor(Option,Name,1),
  149        functor(Term,Name,1),
  150        with_mutex(space_mutex,
  151                   (   retractall(space_setting_aux(Term)),
  152                       assert(space_setting_aux(Option))
  153                   )).
  154*/
  155
  156% FIXME: make bidirectional for settings stored in C++
  157space_setting(Option) :-
  158    with_mutex(space_mutex,space_setting_aux(Option)).
 space_assert(+URI, +Shape, +IndexName) is det
 space_assert(+URI, +Shape) is det
Insert URI with associated Shape in the queue to be inserted into the index with name IndexName or the default index. Indexing happens lazily at the next call of a query or manually by calling space_index/1.
  169space_assert(URI,Shape) :-
  170    space_setting(rtree_default_index(I)),
  171    space_assert(URI,Shape,I).
  172space_assert(URI,Shape,IndexName) :-
  173    dimensionality(Shape,Dimensionality),
  174    must_be(between(1,3), Dimensionality),
  175    (   space_queue(IndexName,retract,_,_)
  176    ->  space_index(IndexName)
  177    ; true
  178    ),
  179    % add URI/Shape only if it is not already queued
  180    % and does not exist in the C++ spatial index
  181    (  rtree_uri_shape(URI, Shape, IndexName)
  182    -> true
  183    ;  space_queue(IndexName,assert,URI,Shape)
  184    -> true
  185    ;  assert(space_queue(IndexName,assert,URI,Shape))
  186    ).
 space_retract(+URI, +Shape, +IndexName) is det
 space_retract(+URI, +Shape) is det
Insert URI with associated Shape in the queue to be removed into the index with name IndexName or the default index. Indexing happens lazily at the next call of a query or manually by calling space_index/1.
  196space_retract(URI,Shape) :-
  197    space_setting(rtree_default_index(I)),
  198    space_retract(URI,Shape,I).
  199space_retract(URI,Shape,IndexName) :-
  200    shape(Shape),
  201    (   space_queue(IndexName,assert,_,_)
  202    ->  space_index(IndexName)
  203    ; true
  204    ),
  205    assert(space_queue(IndexName,retract,URI,Shape)).
 space_index(+IndexName) is det
 space_index is det
Processes all asserts or retracts in the space queue for index IndexName or the default index if no index is specified.
  213space_index :-
  214    space_setting(rtree_default_index(I)),
  215    space_index(I).
  216space_index(IndexName) :-
  217    (   space_queue(IndexName,assert,_,_)
  218    ->  (   empty_nb_set(Assertions),
  219            findall(object(URI,Shape),
  220                    (   space_queue(IndexName,assert,URI,Shape),
  221                        add_nb_set(space_assert(URI,Shape),Assertions)
  222                    ),
  223                    List),
  224            rtree_insert_list(IndexName,List),
  225            retractall(space_queue(IndexName,assert,_,_)),
  226            size_nb_set(Assertions,N),
  227            format('% Added ~w URI-Shape pairs to ~w\n',[N,IndexName])
  228        )
  229    ;   (   space_queue(IndexName,retract,_,_)
  230        ->  (   empty_nb_set(Retractions),
  231                findall(object(URI,Shape),
  232                        (   space_queue(IndexName,retract,URI,Shape),
  233                            add_nb_set(space_retract(URI,Shape),Retractions)
  234                        ),
  235                        List),
  236                rtree_delete_list(IndexName,List),
  237                retractall(space_queue(IndexName,retract,_,_)),
  238                size_nb_set(Retractions,N),
  239                format('% Removed ~w URI-Shape pairs from ~w\n',[N,IndexName])
  240            )
  241        ;   true
  242        )
  243    ).
 space_clear(+IndexName) is det
 space_clear is det
Clears index IndexName or the default index if no index is specified, removing all of its contents.
  251space_clear :-
  252    space_setting(rtree_default_index(I)),
  253    space_clear(I).
  254space_clear(IndexName) :-
  255    retractall(space_queue(IndexName,_,_,_)),
  256    rtree_clear(IndexName).
 space_bulkload(:Closure, +IndexName) is det
 space_bulkload(:Closure) is det
 space_bulkload is det
Fast loading of many Shapes into the index IndexName. Closure is called with two additional arguments: URI and Shape, that finds candidate URI-Shape pairs to index in the index IndexName.

space_bulkload/0 defaults to uri_shape/2 for :Closure.

See also
- the uri_shape/2 predicate for an example of a suitable functor.
  271:- meta_predicate space_bulkload(2), space_bulkload(2,+).  272
  273space_bulkload :-
  274    space_bulkload(uri_shape).
  275space_bulkload(Functor) :-
  276    space_setting(rtree_default_index(I)),
  277    space_bulkload(Functor,I).
  278space_bulkload(Functor, IndexName) :-
  279    once(call(Functor, _Uri, Shape)),
  280    dimensionality(Shape, Dimensionality),
  281    must_be(between(1,3), Dimensionality),
  282%   space_clear(IndexName),  % FIXME: is this ok to skip?
  283    rtree_bulkload(IndexName, Functor, Dimensionality).
 space_contains(+Query, ?Cont, +IndexName) is nondet
 space_contains(+Query, ?Cont) is nondet
Containment query. Unifies Cont with shapes contained in Query Shape (or shape of Query URI) according to index IndexName or the default index.
  292space_contains(Q, Cont) :-
  293    (   shape(Q)
  294    ->  space_shape_contains(Q, Cont)
  295    ;   uri_shape(Q, Shape),
  296        space_shape_contains(Shape, Cont)
  297    ).
  298space_contains(Q, Cont, IndexName) :-
  299    (   shape(Q)
  300    ->  space_shape_contains(Q, Cont, IndexName)
  301    ;   uri_shape(Q, Shape),
  302        space_shape_contains(Shape, Cont, IndexName)
  303    ).
  304
  305space_shape_contains(Shape, Cont) :-
  306    space_setting(rtree_default_index(I)),
  307    space_contains(Shape,Cont,I).
  308space_shape_contains(Shape, Cont, IndexName) :-
  309    shape(Shape),
  310    space_index(IndexName),
  311    (   ground(Cont)
  312    ->  (   bagof(Con, rtree_incremental_containment_query(Shape, Con, IndexName), Cons),
  313            member(Cont, Cons), !
  314        )
  315    ;   rtree_incremental_containment_query(Shape, Cont, IndexName)
  316    ).
 space_intersects(+Query, ?Inter, +IndexName) is nondet
 space_intersects(+Query, ?Inter) is nondet
Intersection query. Unifies Inter with shapes intersecting with Query Shape (or Shape of Query URI) according to index IndexName or the default index. (intersection subsumes containment)
  325space_intersects(Q, Inter) :-
  326    (   shape(Q)
  327    ->  space_shape_intersects(Q, Inter)
  328    ;   uri_shape(Q, Shape),
  329        space_shape_intersects(Shape, Inter)
  330    ).
  331space_intersects(Q, Inter, IndexName) :-
  332    (   shape(Q)
  333    ->  space_shape_intersects(Q, Inter, IndexName)
  334    ;   uri_shape(Q, Shape),
  335        space_shape_intersects(Shape ,Inter, IndexName)
  336    ).
  337
  338space_shape_intersects(Shape, Inter) :-
  339    space_setting(rtree_default_index(I)),
  340    space_shape_intersects(Shape, Inter, I).
  341space_shape_intersects(Shape, Inter, IndexName) :-
  342    shape(Shape),
  343    space_index(IndexName),
  344    (   ground(Inter)
  345    ->  (   bagof(In, rtree_incremental_intersection_query(Shape, In, IndexName), Ins),
  346            member(Inter, Ins), !
  347        )
  348    ;   rtree_incremental_intersection_query(Shape, Inter, IndexName)
  349    ).
 space_nearest(+Query, -Near, +IndexName) is nondet
 space_nearest(+Query, -Near) is nondet
Incremental Nearest-Neighbor query. Unifies Near with shapes in order of increasing distance to Query Shape (or Shape of Query URI) according to index IndexName or the default index.
  359space_nearest(Q, Near) :-
  360    (   shape(Q)
  361    ->  space_shape_nearest(Q, Near)
  362    ;   uri_shape(Q, Shape),
  363        space_shape_nearest(Shape, Near)
  364    ).
  365space_nearest(Q, Near, IndexName) :-
  366    (   shape(Q)
  367    ->  space_shape_nearest(Q, Near, IndexName)
  368    ;   uri_shape(Q, Shape),
  369        space_shape_nearest(Shape, Near, IndexName)
  370    ).
  371
  372space_shape_nearest(Shape, Near) :-
  373    space_setting(rtree_default_index(I)),
  374    space_shape_nearest(Shape, Near, I).
  375space_shape_nearest(Shape, Near, IndexName) :-
  376    shape(Shape),
  377    space_index(IndexName),
  378    rtree_incremental_nearest_neighbor_query(Shape, Near, IndexName).
 space_nearest(+Query, ?Near, +WithinRange, +IndexName) is nondet
 space_nearest(+Query, ?Near, +WithinRange) is nondet
Incremental Nearest-Neighbor query with a bounded distance scope. Unifies Near with shapes in order of increasing distance to Query Shape (or Shape of Query URI) according to index IndexName or the default index. Fails when no more objects are within the range WithinRange.
  390space_nearest_bounded(Q, Near, WithinRange) :-
  391    (   shape(Q)
  392    ->  space_shape_nearest_bounded(Q, Near, WithinRange)
  393    ;   uri_shape(Q,Shape),
  394        space_shape_nearest_bounded(Shape, Near, WithinRange)
  395    ).
  396space_nearest_bounded(Q, Near, WithinRange, IndexName) :-
  397    (   shape(Q)
  398    ->  space_shape_nearest_bounded(Q, Near, WithinRange, IndexName)
  399    ;   uri_shape(Q,Shape),
  400        space_shape_nearest_bounded(Shape, Near, WithinRange, IndexName)
  401    ).
  402
  403space_shape_nearest_bounded(Shape, Near, WithinRange) :-
  404    space_setting(rtree_default_index(I)),
  405    space_nearest_bounded(Shape,Near,WithinRange,I).
  406space_shape_nearest_bounded(Shape, Near, WithinRange, _IndexName) :-
  407    shape(Shape),
  408    ground(Near),
  409    uri_shape(Near,NearShape),
  410    space_distance(Shape,NearShape,Distance),
  411    Distance < WithinRange, !.
  412space_shape_nearest_bounded(Shape, Near, WithinRange, IndexName) :-
  413    space_index(IndexName),
  414    rtree_incremental_nearest_neighbor_query(Shape, Near, IndexName),
  415    (   uri_shape(Near,NearShape,IndexName)
  416    ->  true
  417    ;   uri_shape(Near,NearShape)
  418    ),
  419    space_distance(Shape,NearShape,Distance),
  420    (   ground(WithinRange)
  421    ->  (   Distance > WithinRange
  422        ->  !, fail
  423        ;   true
  424        )
  425    ;   WithinRange = Distance
  426    ).
  427
  428% alias for OGC compatibility
  429space_within_range(Q,Near,WithinRange) :- space_nearest_bounded(Q,Near,WithinRange).
  430space_within_range(Q,Near,WithinRange,IndexName) :- space_nearest_bounded(Q,Near,WithinRange,IndexName).
  431
  432space_display(IndexName) :-
  433    rtree_display(IndexName).
  434
  435space_display_mbrs(IndexName) :-
  436    rtree_display_mbrs(IndexName).
 uri_shape(?URI, ?Shape) is nondet
Finds pairs of URIs and their corresponding Shapes based on WGS84 RDF properties (e.g. wgs84:lat), GeoRSS Simple properties (e.g. georss:polygon), and GeoRSS GML properties (e.g. georss:where).

uri_shape/2 is a dynamic predicate, which means it can be extended. If you use uri_shape/2 in this way, the URI argument has to be a canonical URI, not a QName.

  449uri_shape(URI,Shape) :-
  450    uri_shape(URI,Shape,_).
 uri_shape(?URI, ?Shape, +Source) is nondet
Finds pairs of URIs and their corresponding Shapes using RDF that was loaded from a given Source. Unify Shape with the shape present in the spatial index, if none is found there, try to find the shape using rhe rdf database and the builtin conversions. If the builtin conversions fail, try the rdfuri_shape/3 hook provided by the user.

Otherwise, fails silently.

  466uri_shape(URI,Shape,Source) :-
  467    % load from C++ database
  468    (  var(Source)
  469    -> space_setting(rtree_default_index(Index)),
  470       Source1 = Index
  471    ;  Source1 = Source
  472    ),
  473        % try C++ spatial index
  474    (   rtree_uri_shape(URI, Shape, Source1)
  475    *-> true
  476        % try built-in conversions and hook
  477    ;   uri_shape_(URI,Shape,Source)
  478    ).
  479
  480
  481
  482% Internal predicate that calls hook if available.
  483uri_shape_(URI,Shape,Source) :-
  484    georss_candidate(URI,Shape,Source).
  485uri_shape_(URI,Shape,Source) :-
  486    wgs84_candidate(URI,Shape,Source).
  487uri_shape_(URI,Shape,Source) :-
  488    freebase_candidate(URI,Shape,Source).
  489uri_shape_(URI,Shape,Source) :-
  490    dbpedia_candidate(URI,Shape,Source).
  491uri_shape_(URI,Shape,Source) :-
  492    % user provided hook
  493    rdfuri_shape(URI,Shape,Source).
 rdfuri_shape(URI, Shape)
Hook to convert uri into a shape. This is a hook that the user can provide. This hook is called if the conversions provided by the space library fail. The current conversions provided are:
  516:- multifile rdfuri_shape/3.  517
  518
  519uri_shape_graph(G,U,S) :- uri_shape(U,S,G).
 space_index_all(+IndexName) is det
 space_index_all is det
Loads all URI-Shape pairs found with uri_shape/2 into index IndexName or the default index name.
  528space_index_all :-
  529    space_setting(rtree_default_index(IndexName)),
  530    space_index_all(IndexName).
  531
  532space_index_all(IndexName) :-
  533    space_bulkload(uri_shape,IndexName).
  534
  535
  536/*
  537 * Utility predicates
  538 */
  539
  540box_polygon(box(point(Lx,Ly),point(Hx,Hy)),
  541            polygon([[point(Lx,Ly),point(Lx,Hy),point(Hx,Hy),point(Hx,Ly),point(Lx,Ly)]])).
 shape(+Shape) is det
Checks whether Shape is a valid supported shape.
  547shape(Shape) :- dimensionality(Shape,_).
  548
  549dimensionality(Shape,Dim) :- ground(Shape) -> functor(Shape,point,Dim); !, fail.
  550dimensionality(box(Point,_),Dim) :- dimensionality(Point,Dim).
  551dimensionality(polygon([[Point|_]|_]),Dim) :- dimensionality(Point,Dim).
  552dimensionality(circle(Point,_,_),Dim) :- dimensionality(Point,Dim).
  553dimensionality(linestring([Point|_]),Dim) :- dimensionality(Point,Dim).
  554dimensionality(multipoint([Point|_]),Dim) :- dimensionality(Point,Dim).
  555dimensionality(multipolygon([Poly|_]),Dim) :- dimensionality(Poly,Dim).
  556dimensionality(multilinestring([LS|_]),Dim) :- dimensionality(LS,Dim).
  557dimensionality(geometrycollection([Geom|_]),Dim) :- dimensionality(Geom,Dim).
 space_distance(+Point1, +Point2, -Distance) is det
Calculates the distance between Point1 and Point2 by default using pythagorean distance.
See also
- space_distance_greatcircle/4 for great circle distance.
  566space_distance(X, X, 0).
  567space_distance(A, B, D) :-
  568    (   atom(A)
  569    ->  uri_shape(A, As)
  570    ;   As = A
  571    ),
  572    (   atom(B)
  573    ->  uri_shape(B, Bs)
  574    ;   Bs = B
  575    ),
  576    space_shape_distance(As, Bs, D).
  577
  578space_distance(X, X, 0, _).
  579space_distance(A, B, D, IndexName) :-
  580    (   atom(A)
  581    ->  uri_shape(A, As, IndexName)
  582    ;   As = A
  583    ),
  584    (   atom(B)
  585    ->  uri_shape(B, Bs, IndexName)
  586    ;   Bs = B
  587    ),
  588    space_shape_distance(As, Bs, D, IndexName).
  589
  590space_shape_distance(point(A1,A2), point(B1,B2), D) :-
  591    space_distance_pythagorean(point(A1,A2), point(B1,B2), D), !.
  592space_shape_distance(A, B, D) :-
  593    space_setting(rtree_default_index(IndexName)),
  594    rtree_distance(IndexName, A, B, D1),
  595    pythagorean_lat_long_to_kms(D1,D).
  596space_shape_distance(A, B, D, IndexName) :-
  597    rtree_distance(IndexName, A, B, D1),
  598    pythagorean_lat_long_to_kms(D1,D).
  599
  600% for speed, first assume A and B are shapes, not URIs
  601% if this fails, proceed to interpret them as URIs.
  602space_distance_pythagorean(A, B, D) :-
  603    space_distance_pythagorean_fastest(A, B, D1),
  604    pythagorean_lat_long_to_kms(D1, D).
  605space_distance_pythagorean(A, B, D) :-
  606    (   atom(A)
  607    ->  uri_shape(A, As)
  608    ;   As = A
  609    ),
  610    (   atom(B)
  611    ->  uri_shape(B, Bs)
  612    ;   Bs = B
  613    ),
  614    space_distance_pythagorean(As, Bs, D).
  615
  616space_distance_pythagorean_fastest(point(A, B), point(X, Y), D) :-
  617    D2 is ((X - A) ** 2) + ((Y - B) ** 2),
  618    D is sqrt(D2).
  619
  620pythagorean_lat_long_to_kms(D1, D) :-
  621    D is D1 * 111.195083724. % to kms
 space_distance_greatcircle(+Point1, +Point2, -Dist) is det
 space_distance_greatcircle(+Point1, +Point2, -Dist, +Unit) is det
Calculates great circle distance between Point1 and Point2 in the specified Unit, which can take as a value km (kilometers) or nm (nautical miles). By default, nautical miles are used.
  631space_distance_greatcircle(A, B, D) :-
  632    (   atom(A)
  633    ->  uri_shape(A, As)
  634    ;   As = A
  635    ),
  636    (   atom(B)
  637    ->  uri_shape(B, Bs)
  638    ;   Bs = B
  639    ),
  640    space_shape_distance_greatcircle(As, Bs, D).
  641
  642space_distance_greatcircle(A, B, D, Unit) :-
  643    (   atom(A)
  644    ->  uri_shape(A, As)
  645    ;   As = A
  646    ),
  647    (   atom(B)
  648    ->  uri_shape(B, Bs)
  649    ;   Bs = B
  650    ),
  651    space_shape_distance_greatcircle(As, Bs, D, Unit).
  652
  653
  654space_shape_distance_greatcircle(point(A1,A2), point(B1,B2), D) :-
  655    space_shape_distance_greatcircle(point(A1,A2), point(B1,B2), D, nm).
  656
  657space_shape_distance_greatcircle(point(A1,A2), point(B1,B2), D, km) :-
  658    R is 6371, % kilometers
  659    space_distance_greatcircle_aux(point(A1,A2), point(B1,B2), D, R).
  660space_shape_distance_greatcircle(point(A1,A2), point(B1,B2), D, nm) :-
  661    R is 3440.06, % nautical miles
  662    space_distance_greatcircle_aux(point(A1,A2), point(B1,B2), D, R).
  663
  664% Haversine formula
  665space_distance_greatcircle_aux(point(Lat1deg, Long1deg), point(Lat2deg, Long2deg), D, R) :-
  666    deg2rad(Lat1deg,Lat1),
  667    deg2rad(Lat2deg,Lat2),
  668    deg2rad(Long1deg,Long1),
  669    deg2rad(Long2deg,Long2),
  670    DLat is Lat2 - Lat1,
  671    DLong is Long2 - Long1,
  672    A is (sin(DLat/2)**2) + cos(Lat1) * cos(Lat2) * (sin(DLong/2)**2),
  673    SqA is sqrt(A),
  674    OnemA is 1 - A,
  675    Sq1mA is sqrt(OnemA),
  676    C is 2 * atan(SqA,Sq1mA),
  677    D is R * C.
  678
  679
  680deg2rad(Deg,Rad) :-
  681    Rad is (Deg * pi) / 180.
  682rad2deg(Rad,Deg) :-
  683    Deg is (Rad * 180) / pi.
  684
  685space_bearing(point(Lat1deg, Long1deg), point(Lat2deg, Long2deg), Bearing) :-
  686    deg2rad(Lat1deg,Lat1),
  687    deg2rad(Lat2deg,Lat2),
  688    deg2rad(Long1deg,Long1),
  689    deg2rad(Long2deg,Long2),
  690    DLong is Long2 - Long1,
  691    Y is sin(DLong) * cos(Lat2),
  692    X is cos(Lat1) * sin(Lat2) - sin(Lat1) * cos(Lat2) * cos(DLong),
  693    Bearing0 is atan(Y, X),
  694    rad2deg(Bearing0, Bearing)