18. Geografía

Es muy común tener datos en los que las coordenadas son “geográficas” o “latitud/longitud”.

A diferencia de las coordenadas en Mercator, UTM o Stateplane, las coordenadas geográficas no son coordenadas cartesianas. Las coordenadas geográficas no representan una distancia lineal desde un origen trazado en un plano. Más bien, estas coordenadas esféricas describen coordenadas angulares sobre un globo. En coordenadas esféricas un punto se especifica por el ángulo de rotación desde un meridiano de referencia (longitud), y el ángulo desde el ecuador (latitud).

Puedes tratar las coordenadas geográficas como coordenadas cartesianas aproximadas y continuar realizando cálculos espaciales. Sin embargo, las mediciones de distancia, longitud y área serán sin sentido. Dado que las coordenadas esféricas miden distancia angular, las unidades están en “grados”. Además, los resultados aproximados de índices y pruebas verdadero/falso como intersects y contains pueden resultar completamente erróneos. La distancia entre puntos se vuelve más grande a medida que se aproximan a zonas problemáticas como los polos o la línea internacional de cambio de fecha.

Por ejemplo, aquí están las coordenadas de Los Ángeles y París.

• Los Ángeles: POINT(-118.4079 33.9434)

• París: POINT(2.3490 48.8533)

El siguiente cálculo mide la distancia entre Los Ángeles y París usando la función cartesiana estándar de PostGIS ST_Distance(geometry, geometry). Obsérvese que el SRID 4326 declara un sistema de referencia espacial geográfico.

SELECT ST_Distance(
  'SRID=4326;POINT(-118.4079 33.9434)'::geometry, -- Los Angeles (LAX)
  'SRID=4326;POINT(2.5559 49.0083)'::geometry     -- Paris (CDG)
  );

121.898285970107

¡Ajá! 122. Pero, ¿qué significa eso?

Las unidades para el sistema de referencia espacial 4326 son grados. Así que nuestra respuesta es 122 grados. Pero (otra vez), ¿qué significa eso?

En una esfera, el tamaño de un “cuadrado de un grado” es muy variable, volviéndose más pequeño a medida que uno se aleja del ecuador. Piensa en los meridianos (líneas verticales) del globo acercándose entre sí al ir hacia los polos. Así que, una distancia de 122 grados no significa nada. Es un número sin sentido.

Para calcular una distancia significativa, debemos tratar las coordenadas geográficas no como coordenadas cartesianas aproximadas sino como verdaderas coordenadas esféricas. Debemos medir las distancias entre puntos como trayectorias reales sobre una esfera —una porción de un círculo máximo.

PostGIS proporciona esta funcionalidad a través del tipo geography.

Nota

Diferentes bases de datos espaciales tienen distintos enfoques para “manejar coordenadas geográficas”

• Oracle intenta disimular las diferencias haciendo cálculos geográficos de manera transparente cuando el SRID es geográfico.

• SQL Server usa dos tipos espaciales: “STGeometry” para datos cartesianos y “STGeography” para geográficos.

• Informix Spatial es una extensión puramente cartesiana de Informix, mientras que Informix Geodetic es una extensión puramente geográfica.

• De forma similar a SQL Server, PostGIS usa dos tipos: “geometry” y “geography”.

Usando el tipo geography en lugar de geometry, probemos de nuevo medir la distancia entre Los Ángeles y París.

SELECT ST_Distance(
  'SRID=4326;POINT(-118.4079 33.9434)'::geography, -- Los Angeles (LAX)
  'SRID=4326;POINT(2.5559 49.0083)'::geography     -- Paris (CDG)
  );

9124665.27317673

¡Un número grande! Todos los valores devueltos por cálculos con geography están en metros, así que nuestra respuesta es 9125 km.

Las versiones antiguas de PostGIS soportaban cálculos muy básicos sobre la esfera mediante la función ST_Distance_Spheroid(point, point, measurement). Sin embargo, ST_Distance_Spheroid es bastante limitada. La función solo funciona con puntos y no proporciona soporte para indexación a través de los polos o la línea internacional de cambio de fecha.

La necesidad de soportar geometrías que no son puntos se hace muy clara al plantear una pregunta como: “¿Qué tan cerca pasará un vuelo de Los Ángeles a París de Islandia?”

Trabajar con coordenadas geográficas en un plano cartesiano (la línea púrpura) da una respuesta muy errónea. Usar rutas de círculo máximo (las líneas rojas) da la respuesta correcta. Si convertimos nuestro vuelo LAX-CDG en una line string y calculamos la distancia a un punto en Islandia usando geography, obtendremos la respuesta correcta (en metros).

SELECT ST_Distance(
  ST_GeographyFromText('LINESTRING(-118.4079 33.9434, 2.5559 49.0083)'), -- LAX-CDG
  ST_GeographyFromText('POINT(-22.6056 63.9850)')                        -- Iceland (KEF)
);

502454.906643729

Así, la aproximación más cercana a Islandia (medida desde su aeropuerto internacional) en la ruta LAX-CDG es relativamente pequeña: 502 km.

El enfoque cartesiano para manejar coordenadas geográficas se rompe completamente en el caso de entidades que cruzan la línea internacional de cambio de fecha. La ruta de círculo máximo más corta de Los Ángeles a Tokio cruza el Océano Pacífico. La ruta cartesiana más corta cruza el Atlántico y el Índico.

SELECT ST_Distance(
  ST_GeometryFromText('Point(-118.4079 33.9434)'),  -- LAX
  ST_GeometryFromText('Point(139.733 35.567)'))     -- NRT (Tokyo/Narita)
    AS geometry_distance,
ST_Distance(
  ST_GeographyFromText('Point(-118.4079 33.9434)'), -- LAX
  ST_GeographyFromText('Point(139.733 35.567)'))    -- NRT (Tokyo/Narita)
    AS geography_distance;

 geometry_distance | geography_distance
-------------------+--------------------
  258.146005837336 |   8833954.76996256

18.1. Usando Geography

Para cargar datos de geometría en una tabla de tipo geography, primero la geometría debe proyectarse a EPSG:4326 (longitud/latitud), y luego debe convertirse a geography. La función :command:ST_Transform(geometry,srid) convierte coordenadas a geográficas y la función :command:Geography(geometry) o el sufijo ::geography realiza el casting a geography.

CREATE TABLE nyc_subway_stations_geog AS
SELECT
  ST_Transform(geom,4326)::geography AS geog,
  name,
  routes
FROM nyc_subway_stations;

Construir un índice espacial en una tabla geography es exactamente igual que en una tabla geometry:

CREATE INDEX nyc_subway_stations_geog_gix
ON nyc_subway_stations_geog USING GIST (geog);

La diferencia está en el funcionamiento interno: el índice geography manejará correctamente consultas que crucen los polos o la línea internacional de cambio de fecha, mientras que el de geometry no lo hará.

Aquí hay una consulta para encontrar todas las estaciones de metro dentro de 500 metros del Empire State Building.

WITH empire_state_building AS (
  SELECT 'POINT(-73.98501 40.74812)'::geography AS geog
)
SELECT name,
  ST_Distance(esb.geog, ss.geog) AS distance,
  degrees(ST_Azimuth(esb.geog, ss.geog)) AS direction
FROM nyc_subway_stations_geog ss,
     empire_state_building esb
WHERE ST_DWithin(ss.geog, esb.geog, 500);

Existen solo unas pocas funciones nativas para el tipo geography:

• ST_AsText(geography) devuelve text

• ST_GeographyFromText(text) devuelve geography

• ST_AsBinary(geography) devuelve bytea

• ST_GeogFromWKB(bytea) devuelve geography

• ST_AsSVG(geography) devuelve text

• ST_AsGML(geography) devuelve text

• ST_AsKML(geography) devuelve text

• ST_AsGeoJson(geography) devuelve text

• ST_Distance(geography, geography) devuelve double

• ST_DWithin(geography, geography, float8) devuelve boolean

• ST_Area(geography) devuelve double

• ST_Length(geography) devuelve double

• ST_Covers(geography, geography) devuelve boolean

• ST_CoveredBy(geography, geography) devuelve boolean

• ST_Intersects(geography, geography) devuelve boolean

• ST_Buffer(geography, float8) devuelve geography 1

• ST_Intersection(geography, geography) devuelve geography 1

18.2. Crear una tabla Geography

El SQL para crear una nueva tabla con una columna geography es muy similar al de geometry. Sin embargo, geography permite especificar el tipo de objeto directamente en la creación de la tabla. Por ejemplo:

CREATE TABLE airports (
    code VARCHAR(3),
    geog GEOGRAPHY(Point)
  );

INSERT INTO airports
  VALUES ('LAX', 'POINT(-118.4079 33.9434)');
INSERT INTO airports
  VALUES ('CDG', 'POINT(2.5559 49.0083)');
INSERT INTO airports
  VALUES ('KEF', 'POINT(-22.6056 63.9850)');

En la definición de la tabla, GEOGRAPHY(Point) especifica nuestro tipo de datos de aeropuertos como puntos. Los nuevos campos geography no se registran en la vista geometry_columns, sino en la vista geography_columns.

SELECT * FROM geography_columns;

          f_table_name    | f_geography_column | srid |   type
--------------------------+--------------------+------+----------
 nyc_subway_stations_geog | geog               |    0 | Geometry
 airports                 | geog               | 4326 | Point

Nota

Algunas columnas fueron omitidas en la salida de arriba.

18.3. Convertir a tipo Geometry

Aunque las funciones básicas para geography cubren muchos casos de uso, en ocasiones se necesitan funciones que solo soporta geometry. Afortunadamente, se pueden convertir objetos de geography a geometry y viceversa.

La convención de PostgreSQL para casting es añadir ::typename al final del valor. Por ejemplo: 2::text convierte el número dos en una cadena de texto “2”. Y 'POINT(0 0)'::geometry convierte la representación en texto de un punto a una geometría de tipo punto.

La función :command:ST_X(point) solo admite el tipo geometry. ¿Cómo leer la coordenada X de nuestras geographies?

SELECT code, ST_X(geog::geometry) AS longitude FROM airports;

 code | longitude
------+-----------
 LAX  | -118.4079
 CDG  |    2.5559
 KEF  |  -21.8628

Adjuntando ::geometry a nuestro valor geography lo convertimos a geometry con SRID 4326. Desde ahí podemos usar todas las funciones de geometry que queramos. Pero hay que recordar que, una vez convertido a geometry, las coordenadas se interpretan como cartesianas, no esféricas.

18.4. Por qué (no) usar Geography

Las coordenadas geográficas son universalmente aceptadas: todos entienden qué significan latitud y longitud, pero muy pocos entienden lo que significan coordenadas UTM. ¿Por qué no usar siempre geography?

• Primero, como se señaló antes, hay muchas menos funciones disponibles directamente para geography. Podrías pasar mucho tiempo buscando soluciones a sus limitaciones.

• Segundo, los cálculos en una esfera son mucho más costosos computacionalmente que los cálculos cartesianos. La fórmula cartesiana de distancia (Pitágoras) requiere una llamada a sqrt(). La fórmula esférica (Haversine) requiere dos sqrt(), una arctan(), cuatro sin() y dos cos(). Las funciones trigonométricas son costosas, y los cálculos esféricos involucran muchas de ellas.

¿Conclusión?

Si tus datos son geográficamente compactos (dentro de un estado, condado o ciudad), usa el tipo geometry con una proyección cartesiana adecuada. Para elegir una, visita http://epsg.io e introduce el nombre de tu región para obtener una selección de posibles sistemas de referencia.

Si necesitas medir distancias en un conjunto de datos geográficamente disperso (cubriendo gran parte del mundo), usa el tipo geography. La complejidad que te ahorras compensará los problemas de rendimiento. Además, convertir a geometry puede suplir la mayoría de limitaciones funcionales.

18.5. Lista de funciones

`ST_Distance(geometry, geometry)<http://postgis.net/docs/ST_Distance.html>`_: Para geometry devuelve la distancia cartesiana mínima bidimensional (basada en el sistema de referencia espacial) entre dos geometrías en unidades proyectadas. Para geography devuelve la distancia mínima esferoidal entre dos geographies en metros.

ST_GeographyFromText(text): Devuelve un valor geography a partir de la representación en WKT o extendida.

ST_Transform(geometry, srid): Devuelve una nueva geometría con sus coordenadas transformadas al SRID indicado.

ST_X(point): Devuelve la coordenada X del punto, o NULL si no está disponible. La entrada debe ser un punto.

ST_Azimuth(geography_A, geography_B): Devuelve la dirección de A a B en radianes.

ST_DWithin(geography_A, geography_B, R): Devuelve verdadero si A está dentro de R metros de B.

Nota al pie

1(1,2)

Las funciones de buffer e intersection en realidad son envoltorios sobre un casting a geometry, y no se ejecutan de forma nativa en coordenadas esféricas. Como resultado, pueden fallar al devolver resultados correctos para objetos con extensiones muy grandes que no pueden convertirse limpiamente a una representación planar.

Por ejemplo, la función :command:ST_Buffer(geography,distance) transforma el objeto geography a una proyección “óptima”, aplica el buffer y luego lo transforma de nuevo a geográficas. Si no existe una proyección “óptima” (el objeto es demasiado grande), la operación puede fallar o devolver un buffer malformado.