Geodesic intersection: proposed algorithm and error assessment of current software

Haz clic aquí para ver la versión en español.

I want to show you the summary of a research that we have recently published. We will only focus on the results, for the methodology details you can check the two published articles [1] [2].

As we know to perform any geospatial analysis it is necessary to use a projected coordinate reference system (CRS). The current available geospatial software implements certain operations on the ellipsoid, without using a projected SRS, such as the direct and inverse problems of geodesy, and in this way allows to calculate distances on the ellipsoid, or to calculate a point from an azimuth and a distance. But, these operations are not enough to carry out a full spatial analysis on the ellipsoid, and to address a simple intersection of two segments, or the minimum distance from a point to a line is still solved using a projected CRS, and thereforem any more complex operation such as a full geometric intersection, a buffer, etc. is calculated in the same way.

This way several issues appear:

  • What projected CRS do I use to perform the calculations? This will undoubtedly limit the geographic zone of the analysis.
  • And worst of all, what error am I making when using a projected CRS? And I go further, these unknown errors could be greater than the cartographic precision of my geo data?

In the articles [1] [2] two cases are studied: an intersection of two linear segments, and the minimum distance from a point to a line, but that is, considering that these lines are geodesic. Based on these two cases, it would be possible to build or improve a computational geometry library (e.g. JTS) to perform the full geo spatial analysis on the ellipsoid.

Problem to solve:

Suppose we want to obtain the intersection point of two segments AB and CD each one defined by their two initial and final points (figure 1). Thus, we consider the following table with the coordinates of four different examples:

blog_gi5

The distance AX and CX gives us an idea of the length of these segments, from less than 10 km in case 1, up to more than 10,000 km in case 4.

blog_gi2

Table 1 shows the intersection point (point X) calculated by some of the best-known geo spatial software on the market, some software mention that they work directly on the ellipsoid (supposedly) and others use a projected CRS for the calculation.

blog_gi1

It is scary at first to see how each software obtains different results (point X). So what software gives the most correct results? And, what’s more, do any of them give the exact result? And if not, what error does each software make?

Table 3 shows the point X calculated by the algorithm proposed by the authors of this article. Specifically with the implementation made in PostGIS (STX_GeodesicIntersection function), although any other implementation should give the same results.

As you can see, it is a fast convergence algorithm (see number of iterations). Furthermore, Test B (see article [1] for its definition) calculates the distance from the exact point X to the calculated point. An error of less than 100nm indicates that computational precision has been achieved (calculations with 64-bit floating types).

blog_gi4

And as a final result, we present the table with the errors made by the different software.

See as in case 2 with geodesics of only a few tens of km in length, applications such as PostGIS, ArcGIS, etc. make about 3m of error.

blog_gi3

The PostGIS and Java test implementation can be downloaded from [3].

Finally, mention that, in a normal computer, about 20,000 intersections per second have been obtained. Perhaps, now is the time to migrate certain geospatial libraries such as JTS (GEOS), GDAL, etc. so that they carry out the calculations in full on the ellipsoid and thus be able to forget about the projected CRS, and leave them purely for aesthetic purposes.

[1] https://www.mdpi.com/2076-3417/11/11/5129

Martínez-Llario JC, Baselga S, Coll E. Accurate Algorithms for Spatial Operations on the Spheroid in a Spatial Database Management System. Applied Sciences. 2021; 11(11):5129. https://doi.org/10.3390/app11115129

[2] https://link.springer.com/article/10.1007/s11200-017-1020-z

Baselga, S., Martínez-Llario, J.C. Intersection and point-to-line solutions for geodesics on the ellipsoid. Stud Geophys Geod 62, 353–363 (2018). https://doi.org/10.1007/s11200-017-1020-z

[3] https://github.com/jomarlla/geodesicSpatialOperators

Intersección de geodésicas: algoritmo propuesto y cálculo de errores de los aplicativos actuales

English versión here

Quiero mostraros el resumen de una investigación que hemos publicado recientemente. Solo nos centraremos en los resultados, para las explicaciones de la metodología podéis consultar los dos artículos publicados por nosotros [1] [2].

Como sabemos para realizar un análisis espacial completo es necesario utilizar un sistema de referencia de coordenadas (SRC) proyectado. El panorama del software actual implementa ciertas operaciones sobre el elipsoide, sin utilizar un SRS proyectado, como los problemas directo e inverso de la geodesia, y permite de esta manera calcular distancias sobre el elipsoide, o calcular un punto a partir de un azimut y una distancia. Pero, es que estas operaciones no son suficientes para realizar un análisis espacial completo sobre el elipsoide, y un problema tan sencillo como la intersección de dos segmentos lineales, o la distancia mínima de un punto a una recta se siguen resolviendo utilizando SRC proyectados, y por ende cualquier operación más compleja como una superposición de polígonos, intersección de capas lineales, un área de influencia, etc. se calcula de la misma manera.

Aparecen así varios problemas:

  • ¿Qué SRC proyectado utilizo para realizar los cálculos?, lo cual sin duda limitará la extensión geográfica del análisis.
  • Y lo peor, ¿qué errores estoy cometiendo al utilizar un SRC proyectado? Y voy más allá, ¿estos errores (que no los conozco) podrían ser superior a la precisión cartográfica de mis datos?

En los artículos [1] [2] se estudian dos casos: una intersección de dos segmentos lineales, y la mínima distancia de un punto a una línea, pero eso sí considerando que esas líneas son geodésicas. Apoyándose en estos dos casos sería posible construir o mejorar una biblioteca de geometría computacional (e. g. JTS) para que fuera capaz de calcular un análisis espacial completo sobre el elipsoide.

Problema a resolver:

Supongamos que queremos obtener el punto de intersección de dos segmentos AB y CD definidos cada uno por sus dos puntos inicial y final (figura 1). Así, consideramos el siguiente cuadro con las coordenadas de cuatro ejemplos diferentes:

blog_gi5

La distancia AX y CX nos da una idea de la longitud de dichos segmentos, de menos de 10 km en el ejemplo 1, hasta más de 10000 km en el ejemplo 4.

blog_gi2

La tabla 1, muestra el punto de intersección (punto X) calculado por algunos de los aplicativos y bases de datos espaciales más conocidas del mercado, algunos softwares mencionan que trabajan directamente en el elipsoide (supuestamente) y otros utilizan un SRC proyectado para el cálculo.

blog_gi1

Asusta en un primer momento ver como cada software obtiene diferentes resultados (punto X). Entonces, qué ¿software da los resultados más correctos? Y, es más, ¿alguno de ellos da el resultado exacto?, y si no es así ¿qué error comete cada software?

La tabla 3, muestra el punto X calculado por el algoritmo propuesto por los autores de este artículo. En concreto con la implementación realizada en PostGIS (función STX_GeodesicIntersection), aunque cualquier otra implementación debe dar los mismos resultados.

Como se puede ver es un algoritmo de rápida convergencia (ver número de iteraciones). Además, el Test B (ver artículo [1] para su definición) calcula la distancia del punto exacto X al punto calculado. El error inferior a 100nm indica que se ha llegado a alcanzar la precisión computacional (cálculos con tipos flotantes de 64 bits).

blog_gi4

Y como resultado final, presentamos la tabla con los errores cometidos por los diferentes aplicativos.

Ver como en el caso 2 con geodésicas de solo algunas decenas de km de longitud, aplicativos como PostGIS, ArcGIS, etc. comenten alrededor de 3m de error.

blog_gi3

La implementación de prueba realizada en PostGIS y Java se puede descargar desde [3].

Por último, mencionar que, en un ordenador normal, se han obtenido unas 20000 intersecciones por segundo. Quizás, sea el momento ya, de migrar ciertas librerías geoespaciales como JTS (GEOS), GDAL, etc. para que realicen los cálculos de forma íntegra sobre el elipsoide y así poder olvidarnos de los SRC proyectados, y dejarlos meramente para fines estéticos.

 

[1] https://www.mdpi.com/2076-3417/11/11/5129

Martínez-Llario JC, Baselga S, Coll E. Accurate Algorithms for Spatial Operations on the Spheroid in a Spatial Database Management System. Applied Sciences. 2021; 11(11):5129. https://doi.org/10.3390/app11115129

[2] https://link.springer.com/article/10.1007/s11200-017-1020-z

Baselga, S., Martínez-Llario, J.C. Intersection and point-to-line solutions for geodesics on the ellipsoid. Stud Geophys Geod 62, 353–363 (2018). https://doi.org/10.1007/s11200-017-1020-z

[3] https://github.com/jomarlla/geodesicSpatialOperators