Introducción

Este artículo se corresponde con las páginas 206-211 del libro “PostGIS 2. Análisis Espacial Avanzado”.

Los SIG y las bases de datos espaciales ofrecen generalmente toda su potencia de análisis espacial mediante geometría euclidiana[1], es decir, mediante cálculos de geometría plana donde la distancia más corta entre dos puntos es la recta que los une, entendiendo por ‘recta’ una línea de curvatura cero, ya se esté tratando geometrías bidimensionales (X,Y) o en tres dimensiones.

La realización de cálculos mediante geometría euclidiana tiene sentido cuando las capas cartográficas sobre las cuales aplicamos estos cálculos están georreferenciadas según algún sistema de coordenadas proyectadas u otro sistema de coordenadas planas.

Los CRS horizontales se pueden clasificar en dos grupos:

• Los definidos por coordenadas geodésicas o geográficas que se apoyan en un esfeoride[2] de referencia, como WGS84 (EPSG:4326) o ETRS89 (EPSG:4258). Sus coordenadas generalmente quedan definidas por una pareja de valores de longitud y latitud.
• Los definidos por coordenadas planas (cartesianas) generalmente convertidos desde un CRS definido por coordenadas geodésicas por medio de un determinado artefacto matemático (proyección cartográfica), sus coordenadas quedan definidas por una pareja de valores X, Y o N, E.

Por lo tanto, los cálculos de la mayoría de los programas que utilizan cálculos geoespaciales funcionarán correctamente únicamente en el caso de utilizar CRS definidos por coordenadas planas. Aun así, hay que tener en cuenta las características de la proyección cartográfica y su comportamiento respecto a la conservación bajo ciertas condiciones de las magnitudes angulares, lineales y superficiales.

Las distorsiones de una proyección se pueden medir y representar de forma gráfica utilizando la elipse de Tissot[3] que representa como se transforma un círculo infinitesimal en el esferoide tras ser transformado por una determinada proyección cartográfica.

En cualquier caso, no existe una proyección perfecta, y todas ellas presentan deformaciones de las geometrías proyectadas (ninguna proyección puede conservar los ángulos y las áreas al mismo tiempo en toda la zona representada) que son más grandes cuanto más extensión del terreno abarca el mapa a representar.

La siguiente tabla muestra algunas de las proyecciones más conocidas y sus propiedades (estas proyecciones son las recomendadas según la correspondiente especificación de datos[4] de la directiva Inspire Europea):

Nombre	Propiedades	Propósito general
Lambert Acimutal equivalente[5]	Conserva las áreas. Las distancias y los ángulos (formas) no se conservan.	Mapas cuya representación de áreas verdaderas sea necesaria.
Universal Transversa de Mercator (UTM)	Distorsión lineal máxima ≈ 0.1% Proyección conforme (conserva los ángulos en zonas locales). No distorsiona las superficies en gran medida (latitudes inferiores a 80º).	División en 60 husos (zonas distintas) para representar escalas superiores a 1:500.000 en zonas que se extienden norte-sur.
Lambert Cónica conforme[6]	Conserva los ángulos en las zonas de los paralelos estándar. Según la configuración de los paralelos puede cubrir grandes extensiones con distorsiones en las distancias ≈ 0.01%	Zonas de latitudes medias para representar escalas inferiores a 1:500.000.

Tabla 1 Propiedades de algunas proyecciones cartográficas

En la gran mayoría de las ocasiones se puede encontrar una proyección que se adapte a las precisiones necesarias y que también cubra toda la zona a cartografiar.

Si se necesita cartografiar grandes zonas geográficas como varios países o zonas continentales de forma continua a una escala grande, por ejemplo 1:25.000 no existe ninguna proyección que presente distorsiones inferiores al error cartográfico correspondiente a la escala. Si se opta por utilizar proyecciones UTM habría que utilizar varios husos lo cual rompe la continuidad geométrica y por tanto el análisis espacial entre geometrías en diferentes husos.

Pero entonces, ¿cómo se puede realizar un análisis espacial con este tipo de mapas?

La repuesta consiste en realizar los cálculos directamente en la superficie del esferoide utilizando las coordenadas (generalmente longitud, latitud) de las geometrías. Para realizar estos cálculos no se puede utilizar geometría euclidiana sino geometría diferencial sobre el esferoide utilizado. La diferencia en la complejidad de los cálculos es enorme comparado con el uso de geometría euclídea, un simple cálculo de una distancia entre dos puntos, o una distancia de un punto a una recta, una intersección de dos rectas, etc. viene dada por algoritmos iterativos que en el mejor de los casos requieren varios miles de operaciones comparados con la simplicidad del cálculo euclidiano que se reduce a unas pocas operaciones simples.

Esa ha sido la principal razón por la cual las bibliotecas de manejo de geometrías y de análisis espacial han utilizado y utilizan aún en la actualidad geometría euclídea en sus cálculos y por lo tanto necesitan de cartografía proyectada.

Medida de distancias sobre el esferoide

En los últimos años las bases de datos espaciales han empezado a introducir cálculos geodésicos en algunos de sus métodos espaciales. PostGIS ya introdujo desde versiones tempranas los métodos ST_Distance_spheroid y ST_Length_spheroid para el cálculo geodésico de distancias. Dichas funciones toman como argumentos geometrías con coordenadas en longitud, latitud y devuelven distancias medidas sobre el esferoide especificado en el segundo argumento.

Los tres CRS utilizados 3857, 3575 y 32630 utilizan el mismo elipsoide de referencia WGS84. La única proyección de estas tres que permite medir distancias sin un error considerable es la proyección UTM como se aprecia en el ejemplo al comparar la distancia medida con la real calculada sobre el esferoide por la función ST_Length_Spheroid.

Tipo geography

En la versión 1.5 se creó un nuevo tipo llamado geography para trabajar con las geometrías utilizando cálculos geodésicos, dejando el tipo geometry para trabajar con las geometrías utilizando geometría euclídea.

Para definir una geometría como geography en lugar de geometry PostGIS dispone de los constructores ST_GeogFromText y ST_GeogFromWKB, o también se puede realizar un Cast desde un objeto de tipo geometry al tipo geography con el operador ‘::’.

Constructor ST_GeogFromText:

Utilización de la conversión a geography y comparación de áreas:

Descripción de los campos de la tabla resultado de la consulta anterior:

• Euclidea: Cálculo del área de una geometría en coordenadas geográficas. PostGIS calculará el área mediante geometría euclídea, lo cual carece de sentido físico y da un resultado erróneo no interpretable. En el mismo error incidiremos si tratamos de calcular una distancia, una intersección de dos líneas, o cualquier operación de análisis espacial por simple que sea utilizando geometrías con coordenadas latitud, longitud, es decir, geometrías no proyectadas.
• Euclidea_utm: Al no ser la proyección UTM una proyección que tenga como propiedad conservar las áreas, ésta aunque cercana al valor real presenta un error considerable, en este caso de 0.05%.
• Euclidea_lamberteq: Se ha proyectado las geometrías utilizando una proyección equivalente, la diferencia con la magnitud real discrepa en un 0.003%. Si en lugar de definir la geometría con cuatro puntos la hubiéramos definido con segmentos más cortos la proyección equivalente de Lambert hubiera dado prácticamente el mismo resultado que el área real.
• Esferoide: El cálculo del área se realiza directamente en el esferoide utilizando algoritmos iterativos con base teórica en la geometría diferencial sobre el esferoide.

El comando ST_Area de PostGIS (a partir de la versión 1.5) está sobrecargado de la siguiente manera:

double precision ST_Area(geometry g1);

double precision ST_Area(geography g1);

Es decir, dicho método admite geometrías de tipo geography. Muy pocos de momento son los métodos PostGIS que admiten argumentos de tipo geography. En la versión 1.5 o 2.0 de PostGIS éstos son:

Método	Descripción
ST_Area, ST_Distance, ST_Length	Cálculo del área, distancia y longitud sobre el esferoide (o la esfera si se indica). Las unidades son en metros.
ST_GeogFromText, ST_AsText ST_GeogFromWKB, ST_AsBinary	Constructores (devuelven una geometría de tipo geography) y lectores.
ST_AsGML, ST_AsGeoJSON ST_AsSVG, ST_AsKML	Diferentes conversores a otros formatos de geometría.
ST_CoveredBy, ST_Covers ST_Intersects, ST_DWithin	Predicados espaciales. ST_Intersects utiliza una tolerancia de 0.00001 metros, de forma que cualquier punto más cercano que esa cantidad se considera que interseca.
ST_Buffer, ST_Intersection	Operadores espaciales que admiten el tipo geography pero NO utilizan cálculos sobre el esferoide sino que aplican una transformación de las geometrías a un sistema proyectado, calculan la operación y entonces vuelven a convertir el resultado a coordenadas geográficas.
&&	El operador de comparación de cajas está definido entre tipos geography. Los índices espaciales también están soportados.

Tabla 2 Métodos espaciales que soportan el tipo geography

El tipo geography únicamente soporta el esferoide WGS84 (EPSG: 4326), aunque las geometrías tengan otro CRS éste será ignorado. En futuras versiones de PostGIS se pretende eliminar esta restricción.

Creación de columnas espaciales con el tipo geography

Las columnas de tipo geography no se registran en geometry_columns sino en una nueva vista denominada geography_columns. El tipo geography aprovecha la capacidad typmod de PostgreSQL. La manera de crear una columna geography en una tabla es similar a la creación de una columna de tipo geometry a partir de PostGIS 2.0 (uso de la característica typmod de PostgreSQL).

Una tabla espacial que contenga una o varias columnas de tipo geography se crea utilizando la capacidad typmod de PostGIS, de forma que no existe un método de tipo Addgeographycolumn y la vista geography_columns aparecerá actualizada en todo momento.

Crea una capa nueva ciudadesbis con la columna geom de tipo geography que contenga puntos 3D (PointZ) en el CRS 4326.

Creación de un índice espacial:

s1=# create index ciudadesbis_geom_idx on ciudadesbis using gist (geom);

Una columna de tipo geography se elimina de una tabla espacial y de la vista geography_columns simplemente borrando la columna.

alter table ciudadesbis drop column geom;

Hay que tener presente que siempre se puede convertir una geometría geography a geometry utilizando un Cast explícito ‘::geometry’, pero los métodos PostGIS tratarán dicha geometría con geometría euclídea. De todas formas, existen muchos métodos no relacionados con análisis espacial que se pueden utilizar tras convertir un objeto geography a geometry, algunos de ellos: ST_IsClosed, ST_NumGeometries, ST_NPoints, ST_NDims, ST_Srid, etc.

A fecha de publicación de este libro, pocas aplicaciones de SIG de escritorio o bibliotecas geoespaciales soportan el tipo geography de PostGIS (entre los SIG más conocidos solo QGIS lo soporta actualmente). Es de esperar sin embargo que rapidamente otros programas lo empiecen a soportar por lo menos con capacidades de visualización.

La importación de cartografía con el comando shp2pgsql directamente utilizando el tipo geography se puede realizar con la opción ‘-G’.

Aunque parece ser que los inconvenientes son superiores a las ventajas, la oportunidad de poder tratar la cartografía de forma ‘nativa’ sin las deformaciones que crean las proyecciones puede ser suficiente para inclinar la balanza a favor del tipo geography en ciertos proyectos cartográficos.

Esta entrada es un ejemplo (5 páginas) del libro “PostGIS 2. Análisis Espacial Avanzado”.