1:- module(
2 geo,
3 [
4 geo_area/2, 5 geo_is_shape/1, 6 geo_max/2, 7 geo_min/2, 8 geo_property/1, 9 geo_shape_dimension/2, 10 geo_shape_type/2, 11 geo_transform/2, 12 geo_type/1 13 ]
14).
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_).
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).
63geo_is_shape(Shape) :-
64 geo_shape_type(Shape, _).
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]).
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]).
142geo_property(Property) :-
143 geo_property__(Property),
144 geo_property_(Property).
145
146geo_property__(geos_version(_)).
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])).
169geo_shape_type(shape(_,_,_,Term), Type) :-
170 Term =.. [Type,_],
171 geo_type(Type).
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, 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').
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')
Geospatial predicates
*/