En este apartado se van a realizar algunas pruebas de medición de tiempo de ejecución de predicados y operadores espaciales con el tipo geography.

Es una obviedad que cualquier operación sobre el esferoide que requiera el cálculo de una distancia u obtención de un punto a partir de un azimut empleará mucho más tiempo que el simple cálculo mediante geometría plana que realiza en el tipo geometry.

En este artículo se realiza una simple pero suficiente comparación para que el lector tenga una idea de que orden de magnitud de aumento del tiempo empleado en el cálculo computacional se está hablando.

Para el ejercicio vamos a cargar en PostGIS, las capas rioschegeog (rioschegeog.sql) que es la misma capa que riosche pero al igual que se ha realizado con la capa meteoche y meteochgeog, se ha convertido al tipo geography en 4326.

De la misma forma se carga la capa nucleoschegeog (nucleoschegeog.sql) que contiene los núcleos de población (capa poligonal).

Pues encontrar los ficheros sql de las capas en los datos del libro de PostGIS.

Se van a realizar dos cálculos de obtención de los 5 vecinos más próximos entre dos capas.

• Cálculo A. Mediante consultas laterales y radio de búsqueda (comandos ST_DWithin y ST_Distance). Hemos creado la siguiente consulta (sustituid capa1 y capa2 por meteoche-riosche, nucleosche-riosche, meteochegeog-rioschegeog y nucleoschjegeog-rioschegeog para obtener las consultas de la tabla de resultados. 
• select m.gid as gida, gidb, distance,indice
       from capa1 m, lateral
        ( select r.gid as gidb, st_distance (r.geom, m.geom) as distance,
                                ROW_NUMBER() over () as indice
            from capa2 r
               where st_dwithin (r.geom, m.geom, 100)
            order by st_distance (r.geom, m.geom) limit 5            
        ) as tabla
        order by gida;
• Cálculo B. Mediante el operador <-> de tipo KNN que no necesita radio de búsqueda. Igual que antes sustituid los nombres de las capas por las correspondientes.
  select m.gid as gida, gidb, distance, indice
       from capa1 m, lateral
        ( select r.gid as gidb, st_distance (r.geom, m.geom) as distance,
                                ROW_NUMBER() over () as indice
            from capa2 r
             order by m.geom <-> r.geom limit 5                
        ) as tabla
     order by gida;

Para los cálculos con el tipo geometry se calcula para cada punto de meteoche los 5 tramos de ríos más cercanos y también para cada polígono de nucleosche los 5 tramos de ríos más cercanos.

La misma tarea se realiza con el tipo geography pero con las capas meteochegeog-rioschegeog y nucleoschegeog-rioschegeog.

Se puede ver como en el mejor de los casos el tipo geography es unas 6 veces más lento y en el peor hasta 6200 veces más lento. Cuando se utiliza la capa de núcleos en lugar de las estaciones los resultados son más extremos, ya que cada geometría tiene muchos puntos (no un único punto como cada estación) y por tanto produce un cálculo más intensivo de las distancias.

Con el tipo geography, el método B (KNN) es más lento que el método A (radio de búsqueda) en todos los casos menos uno (radio 10000 m. y capa meteochegeog). Lo cual hace pensar sobre si su uso es adecuado con el tipo geography. Con el tipo geometry ocurre algo parecido pero menos marcado y en menos casos.

Por último, los operadores espaciales ST_Buffer, ST_Intersection, etc. aplicados al tipo geography, etc. además de no dar resultados exactos sobre el esferoide, pueden ejecutarse varias veces más lentos como mínimo debido al doble proceso de transformación de coordenadas que llevan implícito.

English Translation:

In this section, we will perform some tests of the running time of predicates and space operators with the geography type.

It is obvious that any operation on the spheroid that requires the calculation of a distance or obtaining a point from an azimuth will use much more time than the simple calculation using flat geometries (geometry type).

In this article, a simple but sufficient comparison is made so that the reader has an idea about the increasing time used with the calculations (geography type).

For the exercise, we will load in PostGIS, the layers rioschegeog, meteochegeog, nucleoschegeog which are the same layer as riosche, meteoche and nucleosche but they have been transformed to the geography type (srid 4326).

You can download de SQL files from PostGIS Book.

Two calculations will be made to obtain the 5 closest neighbors between two layers.

• Calculation A. Through lateral queries and search radius (commands ST_DWithin and ST_Distance). We have created the following query (just change the layer1 and layer2 names for meteoche-riosche, nucleosche-riosche, meteochegeog-rioschegeog and nucleoschjegeog-rioschegeog to obtain the queries.
  select m.gid as gida, gidb, distance,indice
       from capa1 m, lateral
        ( select r.gid as gidb, st_distance (r.geom, m.geom) as distance,
                                ROW_NUMBER() over () as indice
            from capa2 r
               where st_dwithin (r.geom, m.geom, 100)
            order by st_distance (r.geom, m.geom) limit 5            
        ) as tabla
        order by gida;
• Calculation B. Using the <-> KNN operator that does not need a search radius. As before, replace the names of the layers with the corresponding ones.
  select m.gid as gida, gidb, distance, indice
       from capa1 m, lateral
        ( select r.gid as gidb, st_distance (r.geom, m.geom) as distance,
                                ROW_NUMBER() over () as indice
            from capa2 r
             order by m.geom <-> r.geom limit 5                
        ) as tabla
     order by gida;

For calculations with the geometry type, the 5 closest river are calculated for each meteoche point and for each nucleosche polygon.

The same task is done with the geography type but with the meteochegeog–rioschegeog and nucleoschegeog–rioschegeog layers.

You can see how in the best case the geography type is about 6 times slower and at worst up to 6200 times slower. When the nucleos layer is used instead of the stations, the results are more extreme, since each geometry has many points (not a single point like each station) and therefore produces a more intensive calculation of the distances

With the geography type, the method B (KNN) is slower than method A (search radius) in all cases but one (radius 10000 m and meteochegeog layer). With the geometry type, something similar occurs but less exaggerated.

Finally, the spatial operators ST_Buffer, ST_Intersection, etc. applied to geography type, etc., in addition to not giving exact results on the spheroid, they can be executed several times slower at least due to the double process of transformation of coordinates.

Con la incorporación de los métodos ST_Project y ST_Azimuth más el ya estudiado ST_Distance es posible resolver el problema directo e inverso de la Geodesia de forma sencilla con PostGIS.

With the incorporation of the ST_Project and ST_Azimuth methods plus the already studied ST_Distance it is possible to solve the direct and inverse Geodesy problem in a simple way with PostGIS.

A partir de la versión 2.2.-2.4, PostGIS utiliza la biblioteca GeographicLib [1]  para los cálculos sobre el esferoide en algunos de sus métodos como: ST_Area, ST_Project, ST_Azimuth, ST_Distance, etc. Esta biblioteca aporta más precisión y fiabilidad a PostGIS en este tipo de cálculos sobre el esferoide.

As of version 2.2.-2.4, PostGIS uses the GeographicLib library [1] for calculations on the spheroid in some of its methods such as ST_Area, ST_Project, ST_Azimuth, ST_Distance, etc. This library provides more precision and reliability to PostGIS in this type of calculations on the spheroid.

Resuelve el problema directo de la Geodesia. A partir de un punto A de coordenadas Longitud = 40.5º y Latitud = -60.4º en el esferoide WGS84, calcula las coordenadas de un punto B situado a 500000 m de distancia del A y con un azimut (AB) de 24.6º.

Solve the direct problem of Geodesy. From a point A of coordinates Longitude = 40.5º and Latitude = -60.4º on the spheroid WGS84, calculate the coordinates of a point B located at 500000m distance from A and with an azimuth (AB) of 24.6º.

El problema inverso de la geodesia vendría dado por el cálculo del azimut (AB y BA) y la distancia sobre el esferoide entre dos puntos A y B.

The inverse problem of geodesy would be given by the calculation of the azimuth (AB and BA) and the distance on the spheroid between two points A and B.

[1] La biblioteca GeographicLib ( https://geographiclib.sourceforge.io/ ) implementa cálculos fiables sobre el esferoide. El autor Charles Karney aporta versiones en C, Python, Java, Matlab, etc. así como ejemplos y una buena documentación. Además el sitio web incorpora herramientas online como una muy buena calculadora geodésica: https://geographiclib.sourceforge.io/cgi-bin/GeoConvert

 [1] The GeographicLib library (https://geographiclib.sourceforge.io/) implements reliable calculations on the spheroid. The author Charles Karney provides versions in C, Python, Java, Matlab, etc. as well as examples and good documentation. In addition, the website incorporates online tools such as a very good geodetic calculator: https://geographiclib.sourceforge.io/cgi-bin/GeoConvert

Este artículo ha sido extraido de libro del cual soy autor: PostGIS: Análisis Espacial Avanzado