=Paper=
{{Paper
|id=Vol-3157/paper6
|storemode=property
|title=Navigating the Earth with Pure SPARQL
|pdfUrl=https://ceur-ws.org/Vol-3157/paper6.pdf
|volume=Vol-3157
|authors=Damien Graux
|dblpUrl=https://dblp.org/rec/conf/esws/Graux22
}}
==Navigating the Earth with Pure SPARQL==
https://ceur-ws.org/Vol-3157/paper6.pdf
Navigating the Earth with Pure SPARQL
Damien Graux
Inria, Université Côte d’Azur, CNRS, I3S, France
Abstract
Let’s assume you are on a boat and have to navigate only basic items such as a map or a compass, and
a . . . . . . SPARQL engine! How to set a course? How to compute distances? In this article, we present a set
of SPARQL code-blocks to be used in such situation and more generally in use-cases where practitioners
need to compute mathematical calculus on geo-data represented by (lat,lon) pairs.
Keywords
Navigation computation, Latitude and Longitude, SPARQL 1.1, Mathematical formula, Geodata
1. Introduction
During the past two decades, Semantic Web technologies for the Web have been developed and
it is now possible to produce, share, analyze and interlink large knowledge graphs (sometimes
containing billions of facts) structured using the RDF W3C standard [1]. Additionally, the W3C
has standardized SPARQL [2], the de facto query language dedicated to RDF which has been
more recently improved to add new features, see e.g. [3] for its current version. In parallel,
hundreds of resources have been published online1 and are often accompanied by associated
SPARQL endpoints on which users are able to send queries in order to retrieve pieces of data.
Among the various types of data able to be found within these resources (e.g. dates, texts
in natural language), it is common to have geo-spatial data2 ; in particular, they are usually
represented using their coordinates through the pair (latitude,longitude).
Consequently, use cases sometimes require the computation involving lat/lon points ac-
cording to specific patterns, to e.g. filter by distance. However, in the current version of the
standard3 , only the four basic mathematical operators are available (+, −, *, /) and some
basic predefined functions, such as CEIL or FLOOR. To address this lack in the standard, some
popular evaluators allow extensions to the SPARQL language to cover popular mathematical
functions (e.g. trigonometric operations). Nonetheless, this results in queries specifically built to
be executed by a specific system and which therefore cannot be shared between users without
GeoLD 2022: 5th International Workshop on Geospatial Linked Data co-located with ESWC, May 30 2022, Hersonissos,
Greece
$ damien.graux@inria.fr (D. Graux)
https://dgraux.github.io/ (D. Graux)
0000-0003-3392-3162 (D. Graux)
© 2022 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR
Workshop
CEUR Workshop Proceedings (CEUR-WS.org)
Proceedings
http://ceur-ws.org
ISSN 1613-0073
1
See the current number of listed resources gathered on the Linked Open Data cloud: 1 301 as of May 2021.
https://lod-cloud.net/
2
For example, WikiData contains information about specific places such as cities.
3
https://www.w3.org/TR/2013/REC-sparql11-query-20130321/#expressions
Figure 1: A loxodrome is represented in red, passing through two points; the shortest great circle arc
defined by the same two points is displayed in blue (shortest path).
adapting them beforehand.
Regarding geo-spatial RDF data in particular, the Open Geospatial Consortium has been
active to propose standardized methods to structure geo-data and to retrieve information stored
as such. In 2011, Battle & Kolas presented GeoSPARQL [4]: a geographic query language for
RDF data. As of now, multiple triplestores offer (partial) support of GeoSPARQL [5].
Nevertheless, even if GeoSPARQL provides a wide set of functions4 , most of them focus
on topological geometries. In this study, we want to mainly focus on navigation-related com-
putations on Earth, and this implies to design new specific functions. Moreover, in order to
guarantee interoperability between triplestores, (and to serve as a showcase of SPARQL’s poten-
tial), we choose to make only use of SPARQL 1.1 standard operators. Therefore, in this article,
we provide a set of ready to be used SPARQL instructions to compute various operations on
(lat,lon)-represented points. In short, a SPARQL practitioner can directly copy and paste the
block of code in order to enrich her query with her own additional specific computations.
The rest of the article is organised as follows: first we provide the reader with brief background
notions about Earth-surface geometry (often used by sailors) in Section 2. Second, in Section 3,
we describe how standard SPARQL can be used to compute advanced mathematical functions.
Then, in the main Section 4, we provide the reader with several implementations of useful
geometrical use-cases. Finally, after reminding the related work in Section 5, we conclude in
Section 6.
2. Background on Earth-surface geometry
Great-circle path.
A great circle, also known as an orthodrome, of a sphere is the intersection of the sphere and
a plane that passes through the center point of the sphere. A great circle is the largest circle
4
http://www.opengis.net/def/function/geosparql/
that can be drawn on any given sphere. Any diameter of any great circle coincides with a
diameter of the sphere, and therefore all great circles have the same center and circumference
as each other. For most pairs of distinct points on the surface of a sphere, there is a unique great
circle through the two points5 . Great-circle navigation or orthodromic navigation consists of
navigating a vessel (ship or aircraft) along a great circle. Such routes yield the shortest distance
between two points on the globe.
Rhumb line.
In navigation, a rhumb line, or loxodrome is an arc crossing all meridians of longitude at the
same angle, that is, a path with constant bearing as measured relative to true north.
Comparison.
Figure 1 visually shows the difference between the path following the great-circle (in cyan) and
the one on a rhumb line (in red). Practically, on a great circle, the bearing6 to the destination
point does not remain constant.
3. Geometrical computations with only standard SPARQL
Depending on the use cases, one might have to deal with geospatial data and more particularly
with (latitude,longitude) coordinates. If e.g. the discovery of points is the only thing needed,
then a classic SPARQL query can be built in order to retrieve the relevant pieces of information.
However, practitioners -when dealing with geo-data- often have to compare records based on
the distance between points for example. For these use cases, it can be very hard to implement
the filters directly within the SPARQL queries, as the SPARQL 1.1 standard only allows the four
basic mathematical operators (+, −, *, /). The chosen approach is thereby to (1) either rely
on built-in specific mathematical functions provided by the SPARQL engine itself or (2) to treat
results afterwards (i.e. locally once having obtained results from the engine). The first case
breaks the interoperability of the SPARQL query itself as it becomes engine-dependent. The
second case moves computations out from the SPARQL engine and implies the query-designer
to build a query having a larger scope so to further refine/filter the results later on, losing
therefore in performance and network traffic.
In particular, geodata related computations, since (lat,lon) pairs use a spherical coordinate
system, rely on trigonometric functions most of the time. However, such functions (e.g. cos,
tan. . . ) are not available by default in SPARQL. In order to compute mathematical expressions
in SPARQL, one solution is to approximate the results using Taylor series as presented by Graux
et al. [6]. For example, close to 0, sin can be developed using the following series (which makes
5
Exception: a pair of antipodal points, for which there are infinitely many great circles.
6
In navigation, bearing is the horizontal angle between the direction of an object and another object, or between it
and that of true north.
only use of the 4 “basic” operators):
+∞
∑︁ 𝑥2𝑘+1
sin 𝑥 = (−1)𝑘
(2𝑘 + 1)!
𝑘=0
Easily, this trigonometric function can be coded in SPARQL and using the first terms could be
enough to for example compare/sort entities. The first term of the sin series can be expressed
as follows using SPARQL primitives7 :
VALUES ?2PI {"6.28318530718"^^xsd:double}
BIND ( (?X-?2PI*FLOOR(?X/?2PI) ) AS ?X_in_0_2PI )
BIND ( (1*( 1* ?X_in_0_2PI )/1.0 ) AS ?sin_first_term )
In order to gain in precision, the value of the variable ?X is projected within [0, 2𝜋] taking
advantage of the 2𝜋-periodicity of the sin.
More generally, this approach can be declined to various mathematical functions which admit
Taylor series. This approach has the advantage of being fully compliant with the SPARQL
standard. Nevertheless, it is practically prone to error if one has to write down manually the
SPARQL bindings. Therefore, to ease the process, in Section 4, we provide the readers with the
SPARQL code blocks for several common geospatial computations.
4. Computations on (Latitude,Longitude) Pairs
In the current section, we present a set of formulas and SPARQL-code blocks8 to compute oper-
ations on (lat,lon) pairs. It is worth noting that these formulas are for calculations considering a
spherical earth, i.e. ignoring the ellipsoidal effects.
General SPARQL Bindings.
In the rest of the Section, we present several block of bindings to operate geospatial computations.
In order to factorise the syntax, there are several bindings that could be put upfront the following
blocks, such as constants like 𝜋 or the Earth radius:
# Useful values.
BIND ( xsd:double("3.14159265359") AS ?PI ) # 𝜋 with 11 digits.
BIND ( xsd:double("6.28318530718") AS ?2PI ) # 2𝜋 with 11 digits.
BIND ( xsd:double("6371") AS ?E_radius ) # Earth's radius, in km.
Similarly, we assume that the considered variables in the SPARQL query for respectively the
latitude and the longitude are ?lat and ?lon when represented in degrees and ?lar and ?lor
when in radians. More generally, we can also use the following bindings to convert a variable
in degrees to its radian-equivalent or to designate the deltas between 2 pairs of coordinates, i.e.
?dellar and ?dellor:
7
Relying on MINDS to generate the bindings from the mathematical formula [6].
8
For space limitation, we do not provide all and limit the number of terms in series.
BIND ( (xsd:double(?lat) * ?PI/180) AS ?lar ) # degrees → radians.
# Having two pairs of coordinates, below are the deltas in radians.
BIND ( ((xsd:double(?lat2)-xsd:double(?lat1)) * ?PI/180) AS ?dellar )
BIND ( ((xsd:double(?lon2)-xsd:double(?lon1)) * ?PI/180) AS ?dellor )
4.1. Considering the great-circle
Distance between two points.
In order to calculate the great-circle distance between two points9 , the result is given using the
haversine10 formula:
𝑎 = sin2 (𝑑𝑒𝑙𝑙𝑎𝑟/2) + cos(𝑙𝑜𝑟1 ). cos(𝑙𝑜𝑟2 ). sin2 (𝑑𝑒𝑙𝑙𝑜𝑟/2)
√ √
𝑐 = 2.𝑎𝑡𝑎𝑛2( 𝑎, 1 − 𝑎)
𝑑 = E_radius . 𝑐
It is worth noting that this formula involves 7 trigonometric functions and 2 square roots. In
particular, it requires an 𝑎𝑡𝑎𝑛2 which is defined as follows:
⎧ 𝑦
⎨ 2. arctan √︀ 2
⎪
⎪ if 𝑥 > 0,
𝑥 + 𝑦2 + 𝑥
𝑎𝑡𝑎𝑛2(𝑦, 𝑥) = 𝜋 if 𝑥 < 0 and 𝑦 = 0,
⎪
⎩ undefined
⎪
if 𝑥 = 0 and 𝑦 = 0.
# A
BIND ((0+1*(1* ((?dellar/2)-?2PI*FLOOR((?dellar/2)/?2PI) ))/1.0 +-1*(1*
˓→ ((?dellar/2)-?2PI*FLOOR((?dellar/2)/?2PI) )*
˓→ ((?dellar/2)-?2PI*FLOOR((?dellar/2)/?2PI) )*
˓→ ((?dellar/2)-?2PI*FLOOR((?dellar/2)/?2PI) ))/6.0 )AS ?Asub1)
BIND ((0+1*(1)/1.0 +-1*(1* (?lar1)-?2PI*FLOOR(?lar1)/?2PI) )*
˓→ (?lar1)-?2PI*FLOOR(?lar1)/?2PI) ))/2.0 )AS ?Asub2)
BIND ((0+1*(1)/1.0 +-1*(1* (?lar2)-?2PI*FLOOR(?lar2)/?2PI) )*
˓→ (?lar2)-?2PI*FLOOR(?lar2)/?2PI) ))/2.0 )AS ?Asub3)
BIND ((0+1*(1* ((?dellor/2)-?2PI*FLOOR((?dellor/2)/?2PI) ))/1.0 +-1*(1*
˓→ ((?dellor/2)-?2PI*FLOOR((?dellor/2)/?2PI) )*
˓→ ((?dellor/2)-?2PI*FLOOR((?dellor/2)/?2PI) )*
˓→ ((?dellor/2)-?2PI*FLOOR((?dellor/2)/?2PI) ))/6.0 )AS ?Asub4)
BIND ( ( FLOOR(( (1*(?Asub1)*(?Asub1)) +?Asub2*?Asub3* (1*(?Asub4)*(?Asub4))
˓→ )*10000)/10000 ) AS ?A )
√
# 𝐴
BIND ((0+(1*(((?A)-1)/((?A)+1)))/1.0
˓→ +(1*(((?A)-1)/((?A)+1))*(((?A)-1)/((?A)+1))*(((?A)-1)/((?A)+1)))/3.0
˓→ +(1*(((?A)-1)/((?A)+1))*(((?A)-1)/((?A)+1))*(((?A)-1)/((?A)+1))*
˓→ (((?A)-1)/((?A)+1))*(((?A)-1)/((?A)+1)))/5.0 )AS ?sub1)
9
The shortest distance over the Earth’s surface, considering an “as-the-crow-flies” distance, while ignoring hills and
topography (︂ )︂
10 1 𝜃 1 − cos(𝜃)
The half versed sine: ℎ𝑎𝑣(𝜃) = .𝑣𝑒𝑟𝑠𝑖𝑛𝑒(𝜃) = sin2 = .
2 2 2
BIND ((0+(1)/1.0 +(1*?sub1)/1.0 +(1*?sub1*?sub1)/2.0 )AS ?sub2)
√
# 1−𝐴
BIND ((0+(1*(((1-?A)-1)/((1-?A)+1)))/1.0
˓→ +(1*(((1-?A)-1)/((1-?A)+1))*(((1-?A)-1)/((1-?A)+1))*
˓→ (((1-?A)-1)/((1-?A)+1)))/3.0+(1*(((1-?A)-1)/((1-?A)+1))*
˓→ (((1-?A)-1)/((1-?A)+1))*(((1-?A)-1)/((1-?A)+1))*(((1-?A)-1)/
˓→ ((1-?A)+1))*(((1-?A)-1)/((1-?A)+1)))/5.0 )AS ?sub3)
BIND ((0+(1)/1.0 +(1*?sub3)/1.0 +(1*?sub3*?sub3)/2.0 )AS ?sub4)
√︁ √ √
# ( 𝐴)2 + ( 1 − 𝐴)2
BIND (((?sub2)*(?sub2)+(?sub4)*(?sub4))AS ?sub5)
BIND ((0+(1*(( ?sub5-1)/( ?sub5+1)))/1.0 +(1*(( ?sub5-1)/( ?sub5+1))*(( ?sub5-1)/(
˓→ ?sub5+1))*(( ?sub5-1)/( ?sub5+1)))/3.0 +(1*(( ?sub5-1)/( ?sub5+1))*(( ?sub5-1)/(
˓→ ?sub5+1))*(( ?sub5-1)/( ?sub5+1))*(( ?sub5-1)/( ?sub5+1))*(( ?sub5-1)/(
˓→ ?sub5+1)))/5.0 )AS ?sub6)
BIND ((0+(1)/1.0 +(1*?sub6)/1.0 +(1*?sub6*?sub6)/2.0 )AS ?sub7)
√
⎛ ⎞
𝐴
# 𝑎𝑟𝑐𝑡𝑎𝑛 ⎝ √︁ √
√ √
⎠
( 𝐴)2 + ( 1 − 𝐴)2 + 1 − 𝐴
BIND (((?sub2)/( ?sub7+(?sub4)))AS ?sub8)
BIND ((0+1*(1* ?sub8)/1.0 +-1*(1* ?sub8* ?sub8* ?sub8)/3.0 +1*(1* ?sub8* ?sub8*
˓→ ?sub8* ?sub8* ?sub8)/5.0 )AS ?sub9)
# Grand finale
BIND ( ( FLOOR((2 * ?sub9)*10000)/10000 ) AS ?C )
BIND ( (2 * ?E_radius * ?C) AS ?distance )
As presented displayed above, this haversine formula can be converted in SPARQL bindings,
using series approximations. For clarity reason, we only used the first 3 terms for the series.
Practically, such set of bindings is tedious to write down manually and the use of a tool (such as
MINDS [6]) is much safer to avoid typos and guarantee a SPARQL compliant output.
Initial bearing.
Most of the time, while navigating following a defined course, the current bearing will vary to
follow a great circle path. The following formula computes the initial bearing11 in radians:
𝜃 = 𝑎𝑡𝑎𝑛2(sin(𝑑𝑒𝑙𝑙𝑜𝑟) cos(𝑙𝑎𝑟2 ) , cos(𝑙𝑎𝑟1 ) sin(𝑙𝑎𝑟2 ) − sin(𝑙𝑎𝑟1 ) cos(𝑙𝑎𝑟2 ) cos(𝑑𝑒𝑙𝑙𝑜𝑟))
where (𝑙𝑎𝑟1 , 𝑙𝑜𝑟1 ) is the starting point and (𝑙𝑎𝑟2 , 𝑙𝑜𝑟2 ) the end point.
11
Also known as forward azimuth.
Mid-point.
Similarly, the half-way point (𝑙𝑎𝑡𝑚 , 𝑙𝑜𝑛𝑚 ) of a course from (𝑙𝑎𝑡1 , 𝑙𝑜𝑛1 ) to (𝑙𝑎𝑡2 , 𝑙𝑜𝑛2 ) may be
calculated using:
𝐵𝑥 = cos(𝑙𝑎𝑟2 ). cos(𝑑𝑒𝑙𝑙𝑜𝑟)
𝐵𝑦 = cos(𝑙𝑎𝑟(︂2 ). sin(𝑑𝑒𝑙𝑙𝑜𝑟) )︂
√︁
𝑙𝑎𝑟𝑚 = 𝑎𝑡𝑎𝑛2 sin(𝑙𝑎𝑟1 ) + sin(𝑙𝑎𝑟2 ) , (cos(𝑙𝑎𝑟1 ) + 𝐵𝑥 )2 + 𝐵𝑦 2
𝑙𝑜𝑟𝑚 = 𝑙𝑜𝑟1 + 𝑎𝑡𝑎𝑛2 (𝐵𝑦 , cos(𝑙𝑎𝑟1 ) + 𝐵𝑥 )
Intermediate point.
More generally, an intermediate point at any fraction 𝑓 along the great circle path between two
points (𝑙𝑎𝑡1 , 𝑙𝑜𝑛1 ) and (𝑙𝑎𝑡2 , 𝑙𝑜𝑛2 ) can be calculated:
sin((1 − 𝑓 ).𝛿)
𝑎 =
sin 𝛿
sin(𝑓.𝛿)
𝑏 =
sin 𝛿
𝑥 = 𝑎 cos(𝑙𝑎𝑟1 ) cos(𝑙𝑜𝑟1 ) + 𝑏 cos(𝑙𝑎𝑟2 ) cos(𝑙𝑜𝑟2 )
𝑦 = 𝑎 cos(𝑙𝑎𝑟1 ) sin(𝑙𝑜𝑟1 ) + 𝑏 cos(𝑙𝑎𝑟2 ) sin(𝑙𝑜𝑟2 )
𝑧 = 𝑎 sin(𝑙𝑎𝑟
(︁ 1 ) + 𝑏 sin(𝑙𝑎𝑟)︁2 )
√︀
𝑙𝑎𝑟𝑖 = 𝑎𝑡𝑎𝑛2 𝑧 , 𝑥2 + 𝑦 2
𝑙𝑜𝑟𝑖 = 𝑎𝑡𝑎𝑛2(𝑦 , 𝑥)
𝑑
with 𝛿 the angular distance between the two points. We note that for 𝑓 = 0 and
𝐸_𝑟𝑎𝑑𝑖𝑢𝑠
𝑓 = 1 respectively, we have point 1 and point 2.
Closest point to the poles.
Using Clairaut’s formula12 , we can obtain the maximum latitude 𝑙𝑎𝑟𝑀 of a great circle path,
given a bearing 𝐵 and latitude 𝑙𝑎𝑟 on the great circle:
𝑙𝑎𝑟𝑀 = acos(| sin 𝐵. cos(𝑙𝑎𝑟)|)
Destination point.
Considering an initial point (𝑙𝑎𝑟1 , 𝑙𝑜𝑟1 ), an initial bearing 𝐵 and a distance 𝑑, it is possible to
know the arrival point following a great-circle route:
𝑙𝑎𝑟𝑓 = asin(sin(𝑙𝑎𝑟1 ). cos 𝛿 + cos(𝑙𝑎𝑟1 ). sin 𝛿. cos 𝐵)
𝑙𝑜𝑟𝑓 = 𝑙𝑎𝑟1 + 𝑎𝑡𝑎𝑛2(sin 𝐵. sin 𝛿. cos(𝑙𝑎𝑟1 ) , cos 𝛿 − sin(𝑙𝑎𝑟1 ). sin(𝑙𝑎𝑟𝑓 ))
𝑑
with 𝛿 the angular distance .
𝐸_𝑟𝑎𝑑𝑖𝑢𝑠
12
https://en.wikipedia.org/wiki/Clairaut%27s_relation_(differential_geometry)
4.2. Considering rhumb lines
Most of the features presented for great-circle ways are applicable to rhumb lines too. As
explained in Section 2, rhumb lines differ from great-circle paths in the sense that they consider
constant bearing, crossing all meridians at the same angle. Traditionally, sailors used to navigate
along rhumb lines since it is easier to follow a constant compass bearing than to be continually
adjusting the bearing, as is needed to follow a great circle. In addition, rhumb lines are straight
lines on Mercator projection maps. However, rhumb lines are generally longer than great-circle
routes.
Practically, most of rhumb-related computations are using the inverse Gudermannian func-
tion13 , which gives the height on a Mercator projection map of a given latitude 𝜙:
(︁ (︁ 𝜋 𝜙 )︁)︁
𝑔𝑑−1 : 𝜙 ↦→ 𝑔𝑑−1 (𝜙) = 𝑙𝑛 tan +
4 2
𝜋 𝜙
# + renamed ?X
4 2
BIND ( ( FLOOR((?PI/4 + ?lar/2)*10000)/10000 ) AS ?X )
# Series of sin(?X)
BIND ((0+1*(1* ((?X)-?2PI*FLOOR((?X)/?2PI) ))/1.0 +-1*(1* ((?X)-?2PI*FLOOR((?X)/?2PI) )*
˓→ ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI) ))/6.0 +1*(1*
˓→ ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI)
˓→ )* ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI) ))/120.0 )AS ?sub1)
# Series of cos(?X)
BIND ((0+1*(1)/1.0 +-1*(1* ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI) ))/2.0
˓→ +1*(1* ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI) )*
˓→ ((?X)-?2PI*FLOOR((?X)/?2PI) )* ((?X)-?2PI*FLOOR((?X)/?2PI) ))/24.0 )AS ?sub2)
(︂ )︂
sin(?𝑋)
# Series for 𝑙𝑛(tan(?𝑋)) = 𝑙𝑛
cos(?𝑋)
BIND (( 2*(0+(1*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1)))/1.0
˓→ +(1*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/
˓→ ((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1)))/3.0
˓→ +(1*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/
˓→ ((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/
˓→ ((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1)))/5.0 ))AS ?sub3)
BIND ( ( FLOOR((?sub3)*10000)/10000 ) AS ?Gudermannian )
Distance.
Since a rhumb line is a straight line on a Mercator projection, the distance between two points
along a rhumb line is the length of that line (using Pythagoras’ theorem); but the distortion
of the projection needs to be compensated for. On a constant latitude course (travelling East-
West), this compensation
(︂ is simply)︂cos(𝑙𝑎𝑟); in the general case, it is 𝑑𝑒𝑙𝑙𝑎𝑟/𝑑𝑒𝑙𝑝𝑟𝑜𝑗 where
tan(𝜋/4 + 𝑙𝑎𝑟2 /2)
𝑑𝑒𝑙𝑝𝑟𝑜𝑗 = 𝑙𝑛 i.e. the projected latitude difference. Which leads to the
tan(𝜋/4 + 𝑙𝑎𝑟1 /2)
13
https://en.wikipedia.org/wiki/Gudermannian_function
following formula:
𝑑𝑒𝑙𝑙𝑎𝑟
𝑞 = (or cos(𝑙𝑎𝑟) for E-W lines)
𝑑𝑒𝑙𝑝𝑟𝑜𝑗
𝑑 = 𝑑𝑒𝑙𝑙𝑎𝑟2 + 𝑞 2 .𝑑𝑒𝑙𝑙𝑜𝑟2 . E_radius (Pythagoras)
√︀
Bearing.
As a rhumb line is a straight line on a Mercator projection, with an angle on the projection
equal to the compass bearing.
𝐵 = 𝑎𝑡𝑎𝑛2(𝑑𝑒𝑙𝑙𝑜𝑟, 𝑑𝑒𝑙𝑝𝑟𝑜𝑗)
4.3. Distance to the horizon from the “crow’s nest”
At a height ℎ above the ground, the distance to the horizon 𝑑, is given by: 𝑑 = (2 * 𝑅 * ℎ/𝑏)
√︀
with 𝑏 = 0.8279 is a factor that accounts for atmospheric refraction and depends on the
atmospheric temperature lapse rate, which is taken to be standard14 , leading us to the following
SPARQL bindings:
BIND ( "0.8279" AS ?b )
BIND ( (2*xsd:double(?E_radius)*xsd:double(?h)/xsd:double(?b)) AS ?int )
BIND ((0+(1*(((?int)-1)/((?int)+1)))/1.0
+(1*(((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1)))/3.0
+(1*(((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1))*
˓→ (((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1)))/5.0
+(1*(((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1))*
˓→ (((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1))*
˓→ (((?int)-1)/((?int)+1))*(((?int)-1)/((?int)+1)))/7.0
)AS ?sub1)
BIND ((0+(1)/1.0+(1*?sub1)/1.0+(1*?sub1*?sub1)/2.0 + (1*?sub1*?sub1*?sub1)/6.0 )AS ?sub2)
BIND ( ( FLOOR((?sub2)*10000)/10000 ) AS ?distance )
4.4. Faster approximations for great-circle distances
In practice, if performance is an issue and accuracy less important, for small distances Pythagoras’
theorem15 can be used on an equirectangular projection16 :
(︂ )︂
𝑙𝑎𝑟1 + 𝑙𝑎𝑟2
𝑥 = 𝑑𝑒𝑙𝑙𝑜𝑟 . cos
2
𝑦 = 𝑑𝑒𝑙𝑙𝑎𝑟
𝑑 = E_radius . 𝑥2 + 𝑦 2
√︀
This approximation uses one trigonometric and one square root function (as against the 7
trigonometric plus 2 square roots for haversine fescribed above). Technically, the accuracy is of
14
See Table 12 from the American Practical Navigator [7].
15
https://en.wikipedia.org/wiki/Pythagorean_theorem
16
https://en.wikipedia.org/wiki/Equirectangular_projection
the equirectangular method varies: along meridians there are no errors, otherwise it depends
on distance, bearing, and latitude; however the errors are small enough for many purposes.
Such a method would lead to the following SPARQL code block of bindings:
BIND ((0+1*(1)/1.0 +-1*(1* (((?lar1+?lar2)/2)-?2PI*FLOOR(((?lar1+?lar2)/2)/?2PI) )*
˓→ (((?lar1+?lar2)/2)-?2PI*FLOOR(((?lar1+?lar2)/2)/?2PI) ))/2.0 +1*(1*
˓→ (((?lar1+?lar2)/2)-?2PI*FLOOR(((?lar1+?lar2)/2)/?2PI) )*
˓→ (((?lar1+?lar2)/2)-?2PI*FLOOR(((?lar1+?lar2)/2)/?2PI) )*
˓→ (((?lar1+?lar2)/2)-?2PI*FLOOR(((?lar1+?lar2)/2)/?2PI) )*
˓→ (((?lar1+?lar2)/2)-?2PI*FLOOR(((?lar1+?lar2)/2)/?2PI) ))/24.0 )AS ?sub1)
BIND ( ( FLOOR((xsd:double(?dellor)*?sub1)*10000)/10000 ) AS ?X )
BIND ( ( FLOOR(( (1*?X*?X) + (1*?dellar*?dellar) )*10000)/10000 ) AS ?root )
BIND ((0+(1*(((?root)-1)/((?root)+1)))/1.0
˓→ +(1*(((?root)-1)/((?root)+1))*(((?root)-1)/((?root)+1))*(((?root)-1)/((?root)+1)))/3.0
˓→ +(1*(((?root)-1)/((?root)+1))*(((?root)-1)/((?root)+1))*(((?root)-1)/((?root)+1))
˓→ *(((?root)-1)/((?root)+1))*(((?root)-1)/((?root)+1)))/5.0 )AS ?sqsub1)
BIND ((0+(1)/1.0 +(1*?sqsub1)/1.0 +(1*?sqsub1*?sqsub1)/2.0 )AS ?sqsub2)
BIND ( ( FLOOR((xsd:double(?E_radius)*?sqsub2)*10000)/10000 ) AS ?d )
Alternatively, the polar coordinate “flat-earth” formula can be used. Considering the co-latitudes
𝑐𝑜𝑙𝑎𝑟1 = 𝜋2 − 𝑙𝑎𝑟1 and 𝑐𝑜𝑙𝑎𝑟2 = 𝜋2 − 𝑙𝑎𝑟2 , then:
√︁
𝑑 = E_radius . 𝑐𝑜𝑙𝑎𝑟1 2 + 𝑐𝑜𝑙𝑎𝑟2 2 − 2.𝑐𝑜𝑙𝑎𝑟1 .𝑐𝑜𝑙𝑎𝑟2 . cos(𝑑𝑒𝑙𝑙𝑜𝑟)
5. Related Work
Due to SPARQL’s lack of essential basic math functions in its standard17 , different approaches
have emerged to serve this need.
Some SPARQL evaluators actually do not give the possibility of computing mathematical
functions inside queries at all. This is for instance the case with 4store [8], RDF3X [9] or
SPARQLGX [10] which are nonetheless popular evaluators from the literature known for their
performances evaluating conjunctive queries. However, arguably, the research focus of these
systems was on optimization of joins and indexes and less on feature completeness.
Currently, all practical relevant SPARQL evaluators offer the opportunity of computing
mathematical functions inside the BIND elements and projections. While the SPARQL standard
defines the built-in functions as part of the syntax18 , the widely adopted approach by evaluator
developers is to take advantage of the Function Call rule, which allows arbitrary IRIs to be
used as function names. Hence, function extensions typically require no changes to the SPARQL
syntax. However, the lack of standardization implies two drawbacks:
• Firstly, the namespaces, local names and signatures of functions may vary between
SPARQL engines, which makes it tedious –if not prohibitive– to exchange backends.
• Secondly, the means of computation of a function and therefore the results may differ
between evaluators.
17
Currently, the SPARQL 1.2 Community Group which aims to advance SPARQL functionalities, is describing several
mathematical operators that could be added in the next iteration of the standard. https://github.com/w3c/sparql-12/
18
https://www.w3.org/TR/sparql11-query/#grammar
All popular SPARQL evaluators –often used to serve public endpoints– such as Virtuoso [11],
Jena-Fuseki [12], GraphDB19 and Stardog20 feature mathematical functions, yet, using different
IRIs. For instance, Virtuoso uses the bif: namespace, whereas Stardog reuses the XPath
function namespace21 . Using such an approach of naming differently similar function/operator22
implies a loss of interoperability, especially, it make the design of federated SPARQL queries far
more complex.
Regarding the spatial aspects, several solutions have been proposed by the community
to represent and query geo-data in RDF. For instance, stRDF and stSPARQL [13, 14] were
proposed to deal with spatiotemporal data. Similarly, the Open Geospatial Consortium runs
the GeoSPARQL initiatives [4] to provide a geographic query language for RDF data together
with methods to structure and represent geo-information. Nowadays, considering its adoption
within the community, GeoSPARQL is the de facto solution to deal with geo-RDF-data. Indeed,
multiple evaluators implement GeoSPARQL giving then access to spatial functions for use in
SPARQL queries such as finding a distance or computing a convex hull having sets of points
or geometries. Therefore, GeoSPARQL provides explicit namespaces (and signatures) for the
functions, and the calculations are based on other standards previously defined by the OGC,
which should guarantee the use and the results of these functions from one triplestore to an
other one. Unfortunately, these supports from the triplestores are not yet complete, as shown
by Jovanovik et al. in [5]: “none of these widely used RDF storage solutions fully comply with
the GeoSPARQL standard”. Furthermore, triplestores may come with specific spatial functions
like Jena23 . As a consequence, to guarantee full interoperability of our code blocks and since
some of the computations we wanted were not provided by default, we focused on standard
SPARQL 1.1 operators.
6. Conclusion
In this article, taking the example of navigating/sailing the Earth, we provide the reader with a
list of SPARQL fragments which can be used to compute complex (often tedious and prone-to-
error to write) computations on latitude and longitude pairs while combining trigonometric
functions. To go further, finer-grained and additional code blocks (with typically more terms in
the Taylor series to gain in accuracy) are openly available for copy-pasting purposes on the
following repository:
https://github.com/dgraux/Navigating-with-SPARQL
As our expressions make only use of standard SPARQL 1.1, the suggested code-blocks are able
to be run on all endpoints. In addition, from a query optimisation point of view, our solution
helps to keep as much computation on the endpoint side as possible and therefore gives the
opportunity to the endpoints for bespoke optimisation or caching strategies. Overall, while
19
https://ontotext.com/products/graphdb/
20
https://www.stardog.com/
21
https://www.w3.org/2005/xpath-functions/math#
22
Implementations for built-in STDEV in Virtuoso, Fuseki, Stardog, Sesame: https://gist.github.com/albertmeronyo/
c6ab285d0b73b05392e2f9b8a5bbea82
23
https://jena.apache.org/documentation/geosparql/index.html
basing this study on navigating features on (lat,lon) pairs, we hope it showcases the potential of
SPARQL and that it will foster SPARQL practitioners to rely on endpoints as much as possible.
References
[1] F. Manola, E. Miller, B. McBride, et al., RDF primer, W3C recommendation 10 (2004) 6.
[2] E. Prud’Hommeaux, A. Seaborne, et al., SPARQL query language for RDF, W3C recom-
mendation 15 (2008). Www.w3.org/TR/rdf-sparql-query/.
[3] W3C SPARQL Working Group, et al., SPARQL 1.1 overview, 2013.
Http://www.w3.org/TR/sparql11-overview/.
[4] R. Battle, D. Kolas, GeoSPARQL: enabling a geospatial semantic web, Semantic Web
Journal 3 (2011) 355–370.
[5] M. Jovanovik, T. Homburg, M. Spasić, A GeoSPARQL compliance benchmark, ISPRS
International Journal of Geo-Information 10 (2021) 487.
[6] D. Graux, G. Sejdiu, C. Stadler, G. Napolitano, J. Lehmann, MINDS: a translator to em-
bed mathematical expressions inside SPARQL queries, in: International Conference on
Semantic Systems, Springer, Cham, 2020, pp. 104–117.
[7] N. Bowditch, The american practical navigator, 1995 edition, Deffense Mapping Agency
Hydrographic/Topographic Center, Maryland (1995) 363–371.
[8] S. Harris, N. Lamb, N. Shadbolt, 4store: The design and implementation of a clustered RDF
store, in: 5th International Workshop on Scalable Semantic Web Knowledge Base Systems
(SSWS2009), 2009, pp. 94–109.
[9] T. Neumann, G. Weikum, The RDF-3X engine for scalable management of RDF data, The
VLDB Journal—The International Journal on Very Large Data Bases 19 (2010) 91–113.
[10] D. Graux, L. Jachiet, P. Geneves, N. Layaïda, SPARQLGX: Efficient distributed evaluation
of SPARQL with apache spark, in: International Semantic Web Conference, Springer, 2016,
pp. 80–87.
[11] O. Erling, I. Mikhailov, RDF support in the virtuoso DBMS, in: Networked Knowledge-
Networked Media, Springer, 2009, pp. 7–24.
[12] A. Jena, Apache jena fuseki, The Apache Software Foundation (2014).
[13] M. Koubarakis, K. Kyzirakos, Modeling and querying metadata in the semantic sensor
web: The model stRDF and the query language stSPARQL, in: Extended Semantic Web
Conference, Springer, 2010, pp. 425–439.
[14] K. Kyzirakos, M. Karpathiotakis, M. Koubarakis, Strabon: A semantic geospatial DBMS,
in: International Semantic Web Conference, Springer, 2012, pp. 295–311.