1:- module(
    2  geo,
    3  [
    4    geo_area/2,            % +Shape, -Area
    5    geo_is_shape/1,        % @Term
    6    geo_max/2,             % +Shape, -Maximums
    7    geo_min/2,             % +Shape, -Minimums
    8    geo_property/1,        % ?Property
    9    geo_shape_dimension/2, % +Shape, -Dimension
   10    geo_shape_type/2,      % +Shape, -Type
   11    geo_transform/2,       % +FromShape, -ToShape
   12    geo_type/1             % ?Type
   13  ]
   14).

Geospatial predicates

*/

   20:- use_module(library(apply)).   21:- use_module(library(error)).   22:- use_module(library(lists)).   23:- use_module(library(shlib)).   24:- use_module(library(yall)).   25
   26:- use_module(library(call_ext)).   27:- use_module(library(default)).   28
   29:- use_foreign_library(foreign(geo)).   30
   31:- at_halt(geo_halt_).
 geo_area(+Shape:compound, -Area:float) is det
   37geo_area(shape(_,_,_,'Line'(_)), 0.0) :- !.
   38geo_area(shape(_,_,_,'MultiLine'(_)), 0.0) :- !.
   39geo_area(shape(_,_,_,'MultiPoint'(_)), 0.0) :- !.
   40geo_area(shape(_,_,_,'MultiPolygon'(Polygons)), Area) :- !,
   41  maplist(polygon_area_, Polygons, Areas),
   42  sum_list(Areas, Area).
   43geo_area(shape(_,_,_,'Point'(_)), 0.0) :- !.
   44geo_area(shape(_,_,_,'Polygon'(LineStrings)), Area) :-
   45  polygon_area_(LineStrings, Area).
   46
   47polygon_area_([LineString], Area) :- !,
   48  polygon_area_(LineString, 0.0, Area).
   49polygon_area_([LineString1,LineString2], Area) :-
   50  maplist(polygon_area_, [LineString1,LineString2], [Area1,Area2]),
   51  Area is abs(Area1 - Area2).
   52
   53polygon_area_([[X1,Y1|_],[X2,Y2|T2]|Coords], Sum1, Area) :- !,
   54  Sum2 is Sum1 + (X1 * Y2) - (Y1 * X2),
   55  polygon_area_([[X2,Y2|T2]|Coords], Sum2, Area).
   56polygon_area_(_, Sum, Area) :-
   57  Area is abs(Sum / 2.0).
 geo_is_shape(@Term) is semidet
   63geo_is_shape(Shape) :-
   64  geo_shape_type(Shape, _).
 geo_max(+Shape:compound, -Maximums:list(float)) is det
   70geo_max(shape(_,_,_,Term), Maxs) :-
   71  geo_shape_dimension_(Term, Dim),
   72  geo_max_(Dim, Term, Maxs).
   73
   74geo_max_(2, Term, [X,Y]) :- !,
   75  geo_max_2(Term, X, Y).
   76geo_max_(3, Term, [X,Y,Z]) :- !,
   77  geo_max_3(Term, X, Y, Z).
   78geo_max_(4, Term, [X,Y,Z,LRS]) :-
   79  geo_max_4(Term, X, Y, Z, LRS).
   80
   81geo_max_2('Point'([X,Y]), X, Y) :- !.
   82geo_max_2(Term, X, Y) :-
   83  Term =.. [Type,Terms],
   84  call_must_be(geo_type, Type),
   85  maplist(geo_max_2, Terms, Xs, Ys),
   86  maplist(max_list, [Xs,Ys], [X,Y]).
   87
   88geo_max_3('Point'([X,Y,Z]), X, Y, Z) :- !.
   89geo_max_3(Term, X, Y, Z) :-
   90  Term =.. [Type,Terms],
   91  call_must_be(geo_type, Type),
   92  maplist(geo_max_3, Terms, Xs, Ys, Zs),
   93  maplist(max_list, [Xs,Ys,Zs], [X,Y,Z]).
   94
   95geo_max_4('Point'([X,Y,Z,LRS]), X, Y, Z, LRS) :- !.
   96geo_max_4(Term, X, Y, Z, LRS) :-
   97  Term =.. [Type,Terms],
   98  call_must_be(geo_type, Type),
   99  maplist(geo_max_4, Terms, Xs, Ys, Zs, LRSs),
  100  maplist(max_list, [Xs,Ys,Zs,LRSs], [X,Y,Z,LRS]).
 geo_min(+Shape:compound, -Minimums:list(float)) is det
  106geo_min(shape(_,_,_,Term), Mins) :-
  107  geo_shape_dimension_(Term, Dim),
  108  geo_min_(Dim, Term, Mins).
  109
  110geo_min_(2, Term, [X,Y]) :-
  111  geo_min_2(Term, X, Y).
  112geo_min_(3, Term, [X,Y,Z]) :-
  113  geo_min_3(Term, X, Y, Z).
  114geo_min_(4, Term, [X,Y,Z,LRS]) :-
  115  geo_min_4(Term, X, Y, Z, LRS).
  116
  117geo_min_2('Point'([X,Y]), X, Y) :- !.
  118geo_min_2(Term, X, Y) :-
  119  Term =.. [Type,Terms],
  120  call_must_be(geo_type, Type),
  121  maplist(geo_min_2, Terms, Xs, Ys),
  122  maplist(min_list, [Xs,Ys], [X,Y]).
  123
  124geo_min_3('Point'([X,Y,Z]), X, Y, Z) :- !.
  125geo_min_3(Term, X, Y, Z) :-
  126  Term =.. [Type,Terms],
  127  call_must_be(geo_type, Type),
  128  maplist(geo_min_3, Terms, Xs, Ys, Zs),
  129  maplist(min_list, [Xs,Ys,Zs], [X,Y,Z]).
  130
  131geo_min_4('Point'([X,Y,Z,LRS]), X, Y, Z, LRS) :- !.
  132geo_min_4(Term, X, Y, Z, LRS) :-
  133  Term =.. [Type,Terms],
  134  call_must_be(geo_type, Type),
  135  maplist(geo_min_4, Terms, Xs, Ys, Zs, LRSs),
  136  maplist(min_list, [Xs,Ys,Zs,LRSs], [X,Y,Z,LRS]).
 geo_property(?Property:compound) is nondet
  142geo_property(Property) :-
  143  geo_property__(Property),
  144  geo_property_(Property).
  145
  146geo_property__(geos_version(_)).
 geo_shape_dimension(+Shape:compound, -Dimension:positive_integer) is det
  152geo_shape_dimension(shape(_,_,_,Term), Dim) :-
  153  geo_shape_dimension_(Term, Dim).
  154
  155geo_shape_dimension_('Point'([_]), 1) :- !.
  156geo_shape_dimension_('Point'([_,_]), 2) :- !.
  157geo_shape_dimension_('Point'([_,_,_]), 3) :- !.
  158geo_shape_dimension_('Point'([_,_,_,_]), 4) :- !.
  159geo_shape_dimension_(Term, Dim) :-
  160  Term =.. [Type,Shapes],
  161  call_must_be(geo_type, Type),
  162  maplist(geo_shape_dimension_, Shapes, [Dim|Dims]),
  163  (maplist(=(Dim), Dims) -> true ; type_error(geo_dimensions,[Dim|Dims])).
 geo_shape_type(+Shape:compound, -Type:atom) is det
  169geo_shape_type(shape(_,_,_,Term), Type) :-
  170  Term =.. [Type,_],
  171  geo_type(Type).
 geo_transform(+FromShape:shape, ?ToShape:shape) is det
  177geo_transform(shape(Z,LRS,FromCrs,FromShape), shape(Z,LRS,ToCrs,ToShape)) :-
  178  default_value(ToCrs, 'http://www.opengis.net/def/crs/EPSG/0/4326'),
  179  geo_transform(FromCrs, FromShape, ToCrs, ToShape).
  180
  181geo_transform(FromCrs, 'Point'(FromCoord), ToCrs, 'Point'(ToCoord)) :- !,
  182  geo_transform_coords(FromCrs, [FromCoord], ToCrs, [ToCoord]).
  183geo_transform(FromCrs, 'LineString'(FromCoords), ToCrs, 'LineString'(ToCoords)) :- !,
  184  geo_transform_coords(FromCrs, FromCoords, ToCrs, ToCoords).
  185geo_transform(FromCrs, 'Polygon'(FromLineStrings), ToCrs, 'Polygon'(ToLineStrings)) :- !,
  186  maplist(
  187    {FromCrs,ToCrs}/[FromLineString0,ToLineString0]>>
  188      geo_transform(FromCrs, FromLineString0, ToCrs, ToLineString0),
  189    FromLineStrings,
  190    ToLineStrings
  191  ).
  192geo_transform(FromCrs, FromShape, ToCrs, ToShape) :-
  193  gtrace,%DEB
  194  geo_transform(FromCrs, FromShape, ToCrs, ToShape).
  195
  196geo_transform_coords(FromCrs, FromCoords, ToCrs, ToCoords) :-
  197  maplist(proj_crs, [FromCrs,ToCrs], [FromCrs0,ToCrs0]),
  198  geo_transform_coords_(FromCrs0, FromCoords, ToCrs0, ToCoords).
  199
  200proj_crs('http://www.opengis.net/def/crs/EPSG/0/28992', 'EPSG:28992').
  201proj_crs('http://www.opengis.net/def/crs/EPSG/0/4326', 'EPSG:4326').
 geo_type(+Type:atom) is semidet
geo_type(-Type:atom) is multi
  208geo_type('CircularString').
  209geo_type('CompoundCurve').
  210geo_type('CurvePolygon').
  211geo_type('GeometryCollection').
  212geo_type('LineString').
  213geo_type('MultiCurve').
  214geo_type('MultiLineString').
  215geo_type('MultiPoint').
  216geo_type('MultiPolygon').
  217geo_type('MultiSurface').
  218geo_type('Point').
  219geo_type('Polygon').
  220geo_type('PolyhedralSurface').
  221geo_type('TIN').
  222geo_type('Triangle')