• Webinar Procesamientos GIS Distribuidos
  • ES
  • EN
  • Avda. de Grecia 8. 41012 Sevilla, España
  • |
  • www.geographica.gs

Cayetano Benavent

Calculando la ruta más corta sobre un elipsoide con Python

Líneas Geodésicas y GIS

Escrito el 11/06/2015 por Cayetano Benavent

Una geodésica (o línea geodésica) es el camino más corto posible entre dos puntos ubicados sobre una superficie dada. En el mundo de la cartografía y la geodesia, hallar el camino más corto entre dos ubicaciones sobre un elipsoide de revolución supone calcular la línea geodésica que los une (ver más en Wikipedia).

Considerar la Tierra como una esfera simplifica los cálculos, siendo el camino más corto el arco de círculo máximo entre ambos puntos ( Wikipedia).

Desde hace algún tiempo, una de las principales librerías geo del ecosistema open source, Proj.4,  ha venido introduciendo profundos cambios en su algoritmia de cálculos geodésicos. Concretamente, a partir de la versión 4.9 se ha portado a C dentro de Proj.4 parte de la librería C++ GeographicLib escrita por Charles Karney (más sobre los algoritmos de C. Karney).

Hasta la versión 4.8 esta librería venía utilizando los algoritmos de Paul D. Thomas para ello. El intenso uso que hago diariamente de esta librería me ha conducido a realizar algunos testeos sobre cálculos de distancias geodésicas.

En este enlace podéis acceder a un benchmark que hicimos con varias librerías para calcular la distancia geodésica entre diversas ubicaciones, conociendo las coordenadas de dichas ubicaciones (lo cual conlleva resolver el problema geodésico inverso): https://github.com/cayetanobv/GeodeticMusings

Tras ello, y como apoyo de otras labores que estaba acometiendo, desarrollé una pequeña aplicación que, resolviendo el problema geodésico inverso con base en los algoritmos mencionados, calcula las líneas geodésicas entre los puntos de inicio y fin en formato GIS directamente (Shapefile y/o GeoJSON): GeodesicLinesToGIS.

Es importante resaltar que el programa solventa de forma eficiente el problema de líneas geodésicas que cruzan el antimeridiano a la hora de generar las geometrías en formato GIS. Este programa toma como base las excelentes librerías Pyproj, Fiona y Shapely.

En los ejemplos que mostramos a continuación se puede ver la diferencia entre la línea geodésica calculada (camino más corto; línea verde), la línea de rumbo (o loxodroma; línea roja) y una línea recta entre ambos puntos (línea discontinua negra).

Las diferencias máximas, como es lógico esperar, ocurren entre la proyección Mercator (ver Fig 1; la loxodroma es una línea recta) y la proyección Gnomónica (ver Fig 2; la geodésica es una línea recta). El problema de líneas geodésicas que cruzan la línea de 180º se observa solventado en las Fig 3 y 4.

Fig 1. Proyección Mercator. Líneas Geodésicas
Fig 1. Proyección Mercator.
Fig 2. Proyección Gnomónica (centrada en 50W y 60N). Líneas Geodésicas
Fig 2. Proyección Gnomónica (centrada en 50W y 60N).
Fig 3. Líneas geodésicas que cruzan 180º, Proyección Mercator.
Fig 3. Líneas geodésicas que cruzan 180º, Proyección Mercator.
Fig 4. Líneas geodésicas que cruzan 180º, Proyección Mercator (Centrada en 150E).
Fig 4. Líneas geodésicas que cruzan 180º, Proyección Mercator (Centrada en 150E).

Esta librería puede instalarse directamente desde el repositorio de aplicaciones de Python:

https://pypi.python.org/pypi/GeodesicLinesToGIS

Todo el código fuente está en GiHub:

https://github.com/GeographicaGS/GeodesicLinesToGIS

¡Hasta pronto!

GISOpen sourcePython @es
GEO_blog_cayetano
Escrito el 11/06/2015 por
 Cayetano Benavent
GIS Analyst
COMPARTIR
Fecha: GIS, Open source, Python | Tagged: gis, Líneas Geodésicas, SIG

Análisis de datos con Matplotlib

Mapeando el mundo con Matplotlib

Escrito el 02/03/2015 por Cayetano Benavent

Estas últimas semanas he estado “exprimiendo” un poco más una de mis herramientas favoritas para análisis de datos: Matplotlib. Esta potente librería open source del ecosistema pythoniano, junto a otras como Numpy, Scipy, IPython, Pandas, Sympy, etc., es actualmente core package de las Scientific Computing Tools for Python, más conocidas como The Scipy Stack (http://www.scipy.org/).

A pesar de que llevo años trabajando con Matplotlib, indagar poco a poco en sus entrañas me va permitiendo conocer cada vez más su verdadero potencial.

Como geógrafo, no podía ser de otra manera que el toolkit de Matplotlib que más he usado es Basemap, con el que he pasado muchas horas ploteando resultados antes de integrarlos en un Sistema de Información Geográfica. Cabe destacar que la elevada robustez de Basemap se debe a que se sustenta sobre dos pesos pesados: PROJ4 y GEOS.

Basemap, aunque se encaja directamente encima de Matplotlib como un toolkit, es realmente una librería independiente desarrollada por Jeffrey S. Whitaker, un meteorólogo del Earth System Research Laboratory de la National Oceanic and Atmospheric Administration (NOAA).

Me gustaría recordar la importancia crucial que tiene Jeffrey S. Whitaker en el mundo del software geocientífico, ya que ha regalado a la comunidad otras valiosas “perlas” como el binding Python a la Librería C de NetCDF4 de UNIDATA o los bindings Python a las librerías más importantes de decodificación de ficheros GRIB (NOAA y ECMWF).

Un ejemplo de mis últimos escarceos con Basemap han quedado materializados en una pequeña librería que he construido sobre ella, y que he llamado daynight2geojson. El objetivo de esta pequeña aplicación es obtener la distribución espacial mundial de la noche para una fecha exacta que especifiquemos (que se asume que es UTC). Si no se especifica ninguna fecha, el programa calcula la cobertura nocturna para la fecha actual (UTC, claro está).

Una vez extraídos los datos, se procesan y almacenan en un fichero GeoJSON, que puede por supuesto ser cargado en cualquier cliente capaz de leer este formato (su simpleza hace que hoy día casi cualquier cliente geo lo consuma sin ningún problema).

Todo el procesado de las geometrías se realiza con ayuda de la librería Python para tratar con el formato GeoJSON y la librería para manipulación y análisis de objetos geométricos Shapely (que funciona sobre GEOS). El sistema de referencia de coordenadas (CRS) de la cobertura de salida es EPSG:4326.

Para el que quiera profundizar, todo el fundamento matemático sobre el que descansa el cálculo de la geometría mundial de la noche se encuentra en el módulo solar de Basemap: el cálculo de la GHA (Greenwich hour angle), la declinación solar, etc.

A continuación se muestran algunos ejemplos de cálculos realizados con la librería daynight2geojson para diferentes fechas:

  • Cálculos para el día 15 de enero de 2015 a las 12:00 horas UTC.

Acceso al GeoJSON calculado.

 

  • Cálculos para el día 15 de enero de 2015 a las 18:00 horas UTC (el mismo día a distinta hora).

Acceso al GeoJSON calculado.

Todo el código fuente está accesible en Github: https://github.com/GeographicaGS/daynight2geojson

El que tenga algo que aportar, corregir o discutir, lo puede hacer a través de la apertura de un nuevo issue o forkeándolo y lanzando un pull request.

¡Hasta la próxima!

DatascienceGISOpen sourcePython @es
GEO_blog_cayetano
Escrito el 02/03/2015 por
 Cayetano Benavent
GIS Analyst
COMPARTIR
Fecha: Datascience, GIS, Open source, Python | Tagged: Datascience, geographica

© 2022 Geographica. Todos los derechos reservados.

  • English
  • Español