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. 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 according 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
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 [
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 computations 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 potential), 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.
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 4http://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.
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.
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 trafic.
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. [
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.
In the current section, we present a set of formulas and SPARQL-code blocks8 to compute operations on (lat,lon) pairs. It is worth noting that these formulas are for calculations considering a spherical earth, i.e. ignoring the ellipsoidal efects.
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:
7Relying on MINDS to generate the bindings from the mathematical formula [
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 )
In order to calculate the great-circle distance between two points9, the result is given using the haversine10 formula: = 2.2(√, √1 = E_radius .
− ) = sin2(/2) + cos(1). cos(2). sin2(/2) 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(, ) =
⎧ ⎪ ⎨ ⎪ 2. arctan ⎪ ⎪⎩ undefined
√︀2 + 2 + if > 0, if < 0 and = 0, 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) 9The shortest distance over the Earth’s surface, considering an “as-the-crow-flies” distance, while ignoring hills and topography 10The half versed sine: ℎ( ) =
.( ) = sin2 1 # √︁(√)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 [
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. 11Also known as forward azimuth. calculated using: Similarly, the half-way point (, ) of a course from (1, 1) to (2, 2) may be = cos(2). cos() = cos(2). sin() = 2 sin(1) + sin(2) ,
(cos(1) + )2 + 2 = 1 + 2 ( , cos(1) + ) √︁ ︂)
More generally, an intermediate point at any fraction along the great circle path between two points (1, 1) and (2, 2) can be calculated: = cos(1) cos(1) + cos(2) cos(2) = cos(1) sin(1) + cos(2) sin(2) ︂( = = sin((1 − ). )
sin sin(. ) sin
_
︁( = 2 , √︀2 + 2
︁) = sin(1) + sin(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.
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()|)
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
_
. 12https://en.wikipedia.org/wiki/Clairaut%27s_relation_(diferential_geometry)
Most of the features presented for great-circle ways are applicable to rhumb lines too. As explained in Section 2, rhumb lines difer 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 function13, which gives the height on a Mercator projection map of a given latitude : − 1 : ↦→ − 1( ) = tan ︁( ︁( 4 + )︁ 2 # 4 + 2
renamed ?X # Series of sin(?X) BIND ( ( FLOOR((?PI/4 + ?lar/2)*10000)/10000 ) AS ?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) # Series for (tan(?)) = BIND (( 2*(0+(1*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1)))/1.0 ︂( sin(?) )︂
cos(?) ((?sub1/?sub2)+1))*(((?sub1/?sub2)-1)/((?sub1/?sub2)+1)))/5.0 ))AS ?sub3) BIND ( ( FLOOR((?sub3)*10000)/10000 ) AS ?Gudermannian ) = ︂( tan(/ 4 + 2/2) )︂
tan(/ 4 + 1/2) 13https://en.wikipedia.org/wiki/Gudermannian_function 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 EastWest), this compensation is simply cos(); in the general case, it is / where i.e. the projected latitude diference. Which leads to the following formula: = = √︀2 + 2.2 . E_radius (Pythagoras) (or cos() for E-W lines)
As a rhumb line is a straight line on a Mercator projection, with an angle on the projection equal to the compass bearing.
= 2(, )
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 )
In practice, if performance is an issue and accuracy less important, for small distances Pythagoras’
theorem15 can be used on an equirectangular projection16:
= . cos
=
= E_radius . √︀2 + 2
︂( 1 + 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
14See Table 12 from the American Practical Navigator [
√︁ = E_radius .
12 + 22 − 2.1.2. cos()
Due to SPARQL’s lack of essential basic math functions in its standard17, diferent 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 [
Currently, all practical relevant SPARQL evaluators ofer 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 syntax 18, 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 difer between evaluators. 17Currently, 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/ 18https://www.w3.org/TR/sparql11-query/#grammar
All popular SPARQL evaluators –often used to serve public endpoints– such as Virtuoso [
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 [
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-toerror 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: 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 19https://ontotext.com/products/graphdb/ 20https://www.stardog.com/ 21https://www.w3.org/2005/xpath-functions/math# 22Implementations for built-in STDEV in Virtuoso, Fuseki, Stardog, Sesame: https://gist.github.com/albertmeronyo/ c6ab285d0b73b05392e2f9b8a5bbea82 23https://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.