Archivo del Autor: RealIvanSanchez

The Null Island Algorithm

We geomaticians like to gather around a mythical place called Null Island. This island has everything: airports, train stations, hotels, postcodes, all kinds of shops, a huge lot of geocoded addresses, and whatever geographical feature ends up with null coordinates due to whatever buggy geoprocessing pipeline and ends up in the (0,0) coordinates.

But earlier this year, some geonerds such as @mizmay and @schuyler realized that there is no one Null Island, but one Null Island per datum / coordinate system (depending on who you ask). And @smathermather had the spare time to find out how the “Null Archipielago” looks like:

(Null archipielago image by @smathermather, containing Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under CC BY SA)

Fast forward a few months. I received an e-mail from one of my university peers, asking for help with a puzzle:

A friend of mine received a puzzle with some coordinates. He has to find a place on earth represented by 861126.41, 941711.64.

It’s supposed to be a populated place. Any ideas?

Well, off the top of my head, those looked slightly like UTM coordinates – two digits after the decimal point, suggesting centimeter precision… but the easting is way off the valid range for UTM coordinates.

And I realized this is the Null Archipielago problem, all over again; but instead of plotting (0,0) on a map, let’s plot all points having (861126.41, 941711.64) coordinates in any reference system.

Cue PostGIS. We can create a point in every CRS like so:

select srid, ST_GeomFromText('POINT(861126.41 941711.64)',srid) as geom
 from spatial_ref_sys;

Note the complete absence of PL/SQL in there.

But it will be much easier to work with the data if all the points are in our beloved EPSG:4326 latitude-longitude coordinate system. And while we’re at it, let’s materialize that data into a table:

select srid, ST_Transform(ST_GeomFromText('POINT(861126.41 941711.64)',srid),4326) as geom
 from spatial_ref_sys;

But there is a problem with this – the PostGIS query will crash due to some CRSs having an empty Proj4 string. This took me a while to trace and fix:

select srid, ST_Transform(ST_GeomFromText('POINT(861126.41 941711.64)',srid),4326) as geom
 from spatial_ref_sys where proj4text!='';

And now we can take this data out into a file… but once again, there’s a catch: some of the coordinates are out of bounds and represented by (∞, ∞) coordinate pair. Even though file formats can handle ∞/-∞ values (good thing we know how IEEE floating point format works, right folks?), some mapping software can not accommodate for these values. And I’m looking at you, CartoDB upload page.

In this particular case, there are only points for (∞, ∞) so the data can be cleaned up in just one pass:

delete from archipielago where ST_X(geom)>180;

Then just add a tiny bit of CartoDB magic, and publish a map:


I still don’t know if the original puzzle has anything to do with any obscure used-in-the-real-world CRS, but at least it’s worth a try.

Instalando MapProxy en windows, paso a paso

La semana pasada tuve el placer de formar parte de los formadores de los voluntarios de EUROSHA, un grupo de 25 jóvenes destinados a levantar cartografía en diversos países de África, como parte de las actividades del HOT. Uno de los problemas a los que se enfrentan estos voluntarios es una conexión a internet no muy fiable.

Es perfectamente posible editar datos de OSM offline (guardando los datos a fichero, editando, y resolviendo conflictos de versiones a posteriori), pero lo que no se puede hacer es consultar cartografía de fondo para comparar. Había que hacer algo al respecto. Y la solución fue instalar MapProxy, que permite tomar imágenes ráster de varias fuentes y servirlas como WMS, en local. En un portátil con linux (y python, python-pil y python-pip), instalarlo y probar la configuración por defecto fue una cuestión de minutos.

Ahora bien, los ordenadores que el HOT va a desplegar en África van con windows, principalmente por no disponer del tiempo suficiente para hacer una instalación completa con las herramientas adecuadas para la situación. Improvisemos pues, e instalemos MapProxy tal y como sugiere el manual

We advise you to install MapProxy into a virtual Python environment.

Bueno, pues no hagáis esto. Al instalar python desde cero, lo más probable es que os encontréis con problemas a la hora de instalar las librerías necesarias, en concreto PIL (Python Imaging Library). La manera sencilla de instalar Python para hacer funcionar MapProxy encima es OSGeo4W. Así que descargamos el instalador, elegimos una instalación avanzada, y nos aseguramos de que al menos los paquetes para python y python-pil se van a instalar:

El siguiente paso es descargarse y ejecutarlo dentro de una shell de OSGeo4W como administrador:

En esa misma consola, ejecutamos un easy_install mapproxy, y justo después un easy_install pyproj:

En este punto, los ejecutables de MapProxy ya están instalados. Lo podemos comprobar ejecutando mapproxy-util:

Ahora bien, MapProxy es inútil sin un fichero de configuración que le diga qué servicios tiene que cachear. Así que hacemos copia-pega de una configuración de MapProxy para OpenStreetMap, guardamos el fichero resultante como (por ejemplo) C:\OSGeo4W\mapproxy.yaml, y lanzamos mapproxy-util:

¡Et voilà! Nuestro MapProxy está funcionando y respondiendo a peticiones desde localhost:8080, cacheando tiles de OSM para convertirlas en un servicio WMS:

El resto de opciones se pueden consultar en el manual de MapProxy, pero hay unas cuantas cosas a tener en cuenta:

  • MapProxy siempre debe ejecutarse dentro del entorno de OSGeo4W.
  • … lo que quiere decir que si queremos que se ejecute automáticamente, se puede hacer un .bat haciendo copia-pega de C:\OSGeo4W\osgeo4w.bat, y modificando el comando que se lanza en la última línea de ese script.
  • La utilidad para inicializar o refrescar la caché, mapproxy-seed.exe, ha de ejecutarse también dentro del entorno de OSGeo4W.
  • Los datos cacheados se almacenan en el directorio que se especifique en el fichero de configuración, y es relativo a la ruta donde se lanza mapproxy.