¿Cómo crear un mapa interactivo con Folium?

En esta entrada voy a describir el proceso usado para crear https://kikocorreoso.github.io/datos_aemet/ con la ayuda de la librería Folium.

¿Qué es Folium?

Folium es una librería Python que permite crear mapas interactivos usando Leaflet.js. Lo que hace, de forma elegante, es crear código javascript que usa la maravillosa librería de mapas interactivos leaflet.

¿Crear el mapa?

En la rama gh-pages del repositorio git datos_aemet en github (https://github.com/kikocorreoso/datos_aemet/tree/gh-pages) hay una serie de ficheros. Los ficheros index.html, map.html y readme.html los generaremos a partir de los ficheros:

  • custom.css: Algo de css para que la página cuadre. Por debajo usa, además, bootstrap.
  • src/create_base_map.py: Esta es la madre del cordero y lo que vamos a comentar.
  • src/template.html: Aquí tenemos la estructura principal del HTML usado y que sirve de plantilla a index.html y readme.html.

El fichero src/create_base_map.py hace una serie de cosas:

  • Por un lado lee template.html y modifica una parte del mismo para crear readme.html. Lo importante ocurre de las líneas 11 a 53. Básicamente lee un texto identificativo que he dejado en template.html y lo reemplaza por código explicativo sobre la página.
  • Por otro lado, lee el fichero hdf5 de datos, aemet.h5 y de los datos diarios extrae los metadatos de la estación y las fechas de inicio y fin de los registros. Esta información se formatea e incluye en un marcador usando folium.Marker. Cada marcador tendrá un color en función del periodo de datos disponible. Si una estación tiene, por ejemplo, más de 50 años de registros se incluye en un grupo usando folium.FeatureGroup. Estos grupos de estaciones discretizadas por periodo de medidas se pueden manejar mediante un control que aparece en la parte superior derecha del mapa y que incluimos usando folium.LayerControl. Tanto los grupos de marcadores discretizados por periodos de medidas como el control de capas lo incluimos en el mapa creado usando folium.Map. Cuando tenemos todo colocadito guardamos el mapa en formato html. El guardado del mapa nos crea una página donde el mapa ocupa el 100% de la misma.
  • Por último, la página con el mapa la incluimos en index.html mediante un IFrame usando el mismo template.html que hemos usando para el readme.html.

Creación de un mapa en detalle

Primero los imports

import folium
import branca

branca sirve para ayudarnos a meter HTML en los popup de los marcadores. Sin ello solo he conseguido ver texto plano.

Creamos el mapa:

mi_mapa = folium.Map(location=(39.7, 2.2), zoom_start=8)

Indicamos donde estará centrado el mapa usando location y el nivel de zoom inicial.

Podéis guardar el mapa

mi_mapa.save("mapa.html")

y abrirlo en vuestro navegador favorito y veréis un bonito mapa centrado cerca de unas hermosas islas.

Pero este mapa está muy tímido sin mostrar muchas cosas. Vamos a crear varios marcadores:

# creamos el mapa de nuevo para partir de 0
mi_mapa = folium.Map(location=(39.7, 2.2), zoom_start=8)
# creamos 4 marcadores
marcador1 = folium.Marker(location=(40, 2.1))
marcador2 = folium.Marker(location=(40, 3.5))
marcador3 = folium.Marker(location=(39, 2.1))
marcador4 = folium.Marker(location=(39, 3.5))

Y los incluimos en el mapa y guardamos el mapa:

marcador1.add_to(mi_mapa)
marcador2.add_to(mi_mapa)
marcador3.add_to(mi_mapa)
marcador4.add_to(mi_mapa)
mi_mapa.save("mapa.html")

Si abrís el mapa veréis cuatro marcadores alrededor de la isla de Mallorca. Si pulsáis sobre los marcadores no harán nada.

Vamos a incluir información en un popup y vamos a cambiar el color de los iconos de los marcadores usando folium.Icon:

# creamos el mapa de nuevo para partir de 0
mi_mapa = folium.Map(location=(39.7, 2.2), zoom_start=8)
# La información de los popups la añadiremos usando branca
# La información solo será la posición del marcador
# os dejo a vosotros la innovación
html = "<p>Latitud: 40.0</p><p>Longitud: 2.1</p>"
iframe1 = branca.element.IFrame(html=html, width=500, height=300)
html = "<p>Latitud: 40.0</p><p>Longitud: 3.5</p>"
iframe2 = branca.element.IFrame(html=html, width=500, height=300)
html = "<p>Latitud: 39.0</p><p>Longitud: 2.1</p>"
iframe3 = branca.element.IFrame(html=html, width=500, height=300)
html = "<p>Latitud: 39.0</p><p>Longitud: 3.5</p>"
iframe4 = branca.element.IFrame(html=html, width=500, height=300)
# creamos 4 marcadores y añadimos la información del popup usando folium.Popup
# además, añadimos un icono que será de un color para los marcadores al este
# y de otro color para los marcadores del oeste.
marcador1 = folium.Marker(
    location=(40, 2.1),
    popup=folium.Popup(iframe1, max_width=500),
    icon=folium.Icon(color="black")
)
marcador2 = folium.Marker(
    location=(40, 3.5),
    popup=folium.Popup(iframe2, max_width=500),
    icon=folium.Icon(color="gray")
)
marcador3 = folium.Marker(
    location=(39, 2.1),
    popup=folium.Popup(iframe3, max_width=500),
    icon=folium.Icon(color="black")
)
marcador4 = folium.Marker(
    location=(39, 3.5),
    popup=folium.Popup(iframe4, max_width=500),
    icon=folium.Icon(color="gray")
)
# Añadimos los marcadores al mapa
marcador1.add_to(mi_mapa)
marcador2.add_to(mi_mapa)
marcador3.add_to(mi_mapa)
marcador4.add_to(mi_mapa)
# Y guardamos el mapa
mi_mapa.save("mapa.html")

Por último, vamos a modificar un poco todo esto para añadir los marcadores del este (grises) a una capa y los del oeste (negros) a otra capa y añadir, además, el control de capas. Añado, además, los imports del principio para tener un script completo que podéis modificar a vuestro gusto.

import folium
import branca

# creamos el mapa de nuevo para partir de 0
mi_mapa = folium.Map(location=(39.7, 2.2), zoom_start=8)
# La información de los popups la añadiremos usando branca
# La información solo será la posición del marcador
# os dejo a vosotros la innovación
html = "<p>Latitud: 40.0</p><p>Longitud: 2.1</p>"
iframe1 = branca.element.IFrame(html=html, width=500, height=300)
html = "<p>Latitud: 40.0</p><p>Longitud: 3.5</p>"
iframe2 = branca.element.IFrame(html=html, width=500, height=300)
html = "<p>Latitud: 39.0</p><p>Longitud: 2.1</p>"
iframe3 = branca.element.IFrame(html=html, width=500, height=300)
html = "<p>Latitud: 39.0</p><p>Longitud: 3.5</p>"
iframe4 = branca.element.IFrame(html=html, width=500, height=300)
# creamos 4 marcadores y añadimos la información del popup usando folium.Popup
# además, añadimos un icono que será de un color para los marcadores al este
# y de otro color para los marcadores del oeste.
marcador1 = folium.Marker(
    location=(40, 2.1),
    popup=folium.Popup(iframe1, max_width=500),
    icon=folium.Icon(color="black")
)
marcador2 = folium.Marker(
    location=(40, 3.5),
    popup=folium.Popup(iframe2, max_width=500),
    icon=folium.Icon(color="gray")
)
marcador3 = folium.Marker(
    location=(39, 2.1),
    popup=folium.Popup(iframe3, max_width=500),
    icon=folium.Icon(color="black")
)
marcador4 = folium.Marker(
    location=(39, 3.5),
    popup=folium.Popup(iframe4, max_width=500),
    icon=folium.Icon(color="gray")
)
# Creamos dos grupos para los marcadores
grp_este = folium.FeatureGroup(name='Este')
grp_oeste = folium.FeatureGroup(name='Oeste')
# Añadimos los marcadores AL GRUPO AL QUE CORRESPONDAN (NO AL MAPA)
marcador1.add_to(grp_oeste)
marcador2.add_to(grp_este)
marcador3.add_to(grp_oeste)
marcador4.add_to(grp_este)
# Y ahora añadimos los grupos al mapa
grp_este.add_to(mi_mapa)
grp_oeste.add_to(mi_mapa)
# Y añadimos, además, el control de capas
folium.LayerControl().add_to(mi_mapa)
# Y guardamos el mapa
mi_mapa.save("mapa.html")

Et voilà, tenemos un precioso mapa interactivo con mucha funcionalidad en unas pocas líneas de Python.

Fin del 'así se hizo' de https://kikocorreoso.github.io/datos_aemet/.

Usando google earth con ayuda de python y pykml (II)

Esta es la segunda parte de la entrada 'Usando google earth con ayuda de python y pykml (I)'. En el ejemplo de hoy vamos a ver algo más relacionado con la ciencia que lo visto anteriormente, que no era más que un aburrido ejemplo para introduciros en google earth, kml/kmz y pykml.

[Para esta entrada se ha usado numpy 1.6.1, matplotlib 1.1.1rc, netCDF4 1.0.2, pykml 0.1.0, lxml 2.3.2 (es un prerrequisito de pykml) en el notebook de IPython 0.13 corriendo todo sobre python 2.7.3]

Primero importamos una serie de cosas que usaremos.

import os
from zipfile import ZipFile
import datetime as dt
import numpy as np
from matplotlib import pyplot as plt
import netCDF4 as nc
from lxml import etree

En este caso vamos a representar el huracán Iván, que se generó en 2004 en el Atlántico occidental y llegó a convertirse en un huracán muy potente. Vamos a representar la evolución de la presión al nivel del mar y las zonas de viento superiores a 50 km/h a medida que vamos viendo la trayectoria del huracán. El valor de viento se obtiene según los datos que vamos a usar, estos datos no son muy rigurosos en este sentido por lo que esto solo sirve como ejemplo.

Descargaremos datos de reanálisis del centro europeo de predicción a medio plazo (ECMWF, por sus siglas en inglés, European Center for Medium-range Weather Forecasts) para las fechas de ocurrencia del huracán. Los datos usados en el ejemplo los podéis descargar de aquí (al fichero le he puesto extensión .gif para poder subirlo a wordpress pero es un fichero netcdf, para que funcione el ejemplo deberéis quitar la extensión gif o cambiar el nombre del fichero en el código).

Ahora importamos todo lo necesario de pykml para este ejemplo:

from pykml.factory import nsmap
from pykml.factory import KML_ElementMaker as KML
from pykml.factory import GX_ElementMaker as GX

Ahora vamos a detallar la trayectoria del huracán con su posición, fecha y a en qué categoría se encontraba en cada fecha. La categoría del huracán puede ser de 1 a 5 siendo el 5 el más extremo, según la escala Saffir-Simpson. En la categoría del huracán también vamos a detallar también los momentos en que no es huracán como tal sino que es tormenta tropical o depresión tropical.

# Desde la versión 5.0, google ha añadido una serie de extensiones
# a kml que lo hacen más potente e interactivo
# https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions
# Definimos una variable para el espacio de nombres de las
# extensiones de google que vamos a usar.
gxns = '{' + nsmap['gx'] + '}'
# Ponemos la información referente al  huracán:
# La posición en longitud cada 6h desde que es tormenta tropical
lon_hur = [-30.3, -32.1, -33.6, -35, -36.5, -38.2, -39.9,
           -41.4, -43.4, -45.1, -46.8, -48.5, -50.5, -52.5,
           -54.4, -56.1, -57.8, -59.4, -61.1, -62.6, -64.1,
           -65.5, -67, -68.3, -69.5, -70.8, -71.9, -72.8,
           -73.8, -74.7, -75.8, -76.5, -77.6, -78.4, -79,
           -79.6, -80.4, -81.2, -82.1, -82.8, -83.5, -84.1,
           -84.7, -85.1, -85.6, -86, -86.5, -87, -87.4,
           -87.9, -88.2, -88.2, -87.9, -87.7, -87.4, -86.5,
           -85.7, -84, -82.3, -80.5, -78.5, -76.7, -75.5,
           -74, -74, -74.5, -75.8, -77.5, -78.5, -78.7, -79.1,
           -79.7, -80.6, -81.7, -82.8, -84.1, -86.1, -87.3,
           -88.6, -89.5, -91, -92.2, -92.7, -93.2, -94.2]
# La posición en latitud cada 6h desde que es tormenta tropical
lat_hur = [9.7, 9.5, 9.3, 9.1, 8.9, 8.9, 9, 9.3, 9.5, 9.8,
          10.2, 10.6, 10.8, 11, 11.3, 11.2, 11.3, 11.6,
          11.8, 12, 12.3, 12.6, 13, 13.3, 13.7, 14.2, 14.7,
          15.2, 15.7, 16.2, 16.8, 17.3, 17.4, 17.7, 18,
          18.2, 18.4, 18.8, 19.1, 19.5, 19.9, 20.4, 20.9,
          21.6, 22.4, 23, 23.7, 24.7, 25.6, 26.7, 27.9,
          28.9, 30, 31.4, 32.5, 33.8, 34.7, 35.4, 36.2,
          37, 37.7, 38.4, 38, 37.5, 36, 34.5, 32.8, 31,
          29, 27.5, 26.4, 26.1, 25.9, 25.8, 25.2, 24.8,
          25.1, 26, 26.5, 27.1, 27.9, 28.9, 29.2, 29.6, 30.1]
# El estado del huracan cada 6h
tipo = ['TS','TS','TS','TS','TS','TS','TS',
        'TS','H1','H2','H3','H4','H3','H3',
        'H2','H2','H2','H3','H3','H4','H4',
        'H4','H4','H4','H5','H5','H4','H4',
        'H4','H4','H4','H4','H4','H4','H5',
        'H5','H4','H4','H4','H5','H5','H5',
        'H5','H5','H5','H4','H4','H4','H4',
        'H4','H4','H3','H3','H1','TS','TD',
        'TD','TD','TD','TD','TD','TD','TD',
        'TS','TS','TS','TS','TS','TS','TD',
        'TD','TD','TD','TD','TD','TD','TD',
        'TD','TD','TS','TS','TS','TS','TD','TD']
# La fecha cada 6h desde el 21/09/1998 a las 18Z
fecha_inicial = dt.datetime(2004, 9, 3, 6, 0, 0)
deltat = dt.timedelta(hours = 6)
fechas = [fecha_inicial + deltat * i for i in range(len(tipo))]
# El icono que vamos a asignar a cada 6h, en función del estado de la tormenta.
iconos = {'TD':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/td.gif',
          'TS':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/ts.gif',
          'H1':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h1.gif',
          'H2':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h2.gif',
          'H3':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h3.gif',
          'H4':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h4.gif',
          'H5':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h5.gif'}

Ahora vamos a definir una función que es la que usaremos para dibujar los datos de presión y viento que luego irán en el fichero kml/kmz

# Función que crea las imágenes
def pinta_mapas(nombre, lon, lat, mslp, u, v, loni, lati):
    fig = plt.figure()
    fig.set_size_inches(8, 8. * lat.shape[0] / lon.shape[0])
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    fig.add_axes(ax)
    ax.contour(lon, lat, mslp / 100., np.arange(900,1100.,2.),
               colors='y',linewidths=2, aspect='normal')
    wspd = np.sqrt(u * u + v * v) * 3.6
    ax.contourf(lon, lat, wspd, np.arange(50, 250, 25),
                colors = plt.hot(), aspect='normal')
    plt.savefig(nombre, dpi = 80, transparent=True)

Y ahora vamos a leer los datos que os hayáis descargado:

# Leemos los datos de presión, velocidad, latitud y longitud.
data = nc.Dataset('Usando google earth con ayuda de python y pykml (II)/data.nc')
lon = data.variables['longitude'][:]
lat = data.variables['latitude'][:]
# presión = mslp por mean sea level pressure,
# presión media al nivel del mar
mslp = data.variables['msl'][1:-2,:,:]
lon[lon > 180] = lon - 360
u = data.variables['u10'][1:-2,:,:]
v = data.variables['v10'][1:-2,:,:]

Y ahora es cuando vamos a crear el fichero kml haciendo un bucle donde se van a ir añadiendo cosas tanto a la carpeta (folder) como al Tour. Voy explicándo lo que hace cada cosa mediante comentarios en el código:

# Como vamos a crear un fichero kmz que incluirá todas las imágenes
# abrimos el fichero kmz donde iremos guardándolo todo
fich_kmz = ZipFile('HuracanIvan.kmz', 'w')
# Desde la versión 5.0, google ha añadido una serie de extensiones
# a kml que lo hacen más potente e interactivo
# https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions
# Definimos una variable para el espacio de nombres de las
# extensiones de google que vamos a usar.
gxns = '{' + nsmap['gx'] + '}'
# start with a base KML tour and playlist
fich_kml = KML.kml(
    KML.Document(
      GX.Tour(
        KML.name(u"¡Reprodúceme!"),
        GX.Playlist(),
      ),
      KML.Folder(
        KML.name('Ruta del huracan Ivan'),
        id='lugares',
      ),
    )
)
# Hacemos un bucle para recorrer todos los 'momentos'
# del huracán definidos en el bloque de código anterior
for i in range(len(tipo)):
    # Primero volamos hasta cada posición del huracán (cada 6h)
    # y nos quedamos observándolo desde el espacio
    # A unos 2000 km de altura.
    # Para ello usamos FlyTo donde especificamos
    # La duración en segundos, si queremos que vaya más rápido o más
    # lento, el modo de vuelo, que puede ser más suave (smooth) o
    # más rápido (bounce) y la posición.
    # La posición la determinan longitude, latitude y altitude,
    # tilt es el ángulo desde la vertical (en este caso lo vemos
    # verticalmente (0º)) y range es la distancia desde la posición
    # determinada por longitude, latitude y altitude, teniendo en
    # cuenta el ángulo tilt (en este caso hemos puesto 2.000 km)
    # Las unidades para latitude y longitude son grados, para altitude
    # y range son metros y para tilt son grados desde la vertical.
    # En el gráfico de este enlace se verá mejor:
    # https://developers.google.com/kml/documentation/images/lookAt.gif
    fich_kml.Document[gxns+"Tour"].Playlist.append(
      GX.FlyTo(
        GX.duration(1),
        GX.flyToMode("bounce"),
        KML.LookAt(
          KML.longitude(lon_hur[i]),
          KML.latitude(lat_hur[i]),
          KML.altitude(0),
          KML.heading(0),
          KML.tilt(0),
          KML.range(10000000.0),
          KML.altitudeMode("relativeToGround"),
        )
      ),
    )
    fich_kml.Document[gxns+"Tour"].Playlist.append(GX.Wait(GX.duration(0.1)))
    # Añadimos información de cada lugar en la carpeta
    # como 'placemarks' o marcas de posición. En el nombre y en la descripción
    # del placemark se puede meter HTML y CSS, por si queréis ser un poco más
    # creativos y dejar la información de una forma visible y bonita.
    # extrude == 1 indica que el nombre del lugar se hará visible. Si queremos
    # hacerlo visible solo cuando lo seleccionemos pues le ponemos el valor 0
    desc = '<![CDATA[{0}
lon: {1}
lat: {2}
Estado del huracan: {3}]]>'
    desc = desc.format(fechas[i].isoformat()[0:13], lon_hur[i], lat_hur[i], tipo[i])
    fich_kml.Document.Folder.append(
      KML.Placemark(
        KML.name(fechas[i].isoformat()[0:13]),
        KML.description(
            "{0}<br/>{1}".format('IVAN', desc,)
        ),
        KML.Style(
          KML.IconStyle(
            KML.Icon(KML.href(iconos[tipo[i]]),)
          ),
        ),
        KML.Point(
          KML.extrude(1),
          KML.altitudeMode("relativeToGround"),
          KML.coordinates("{0},{1},0".format(lon_hur[i],
                                             lat_hur[i],)
          ),
        ),
        id=fechas[i].isoformat()[0:13]
      )
    )
    # Aquí le indicamos a la reproducción que nos
    # muestre el 'globo' o 'viñeta' con la información del sitio,
    # nombre y descripción. Visibility a 1 para
    # hacerlo visible.
    fich_kml.Document[gxns+"Tour"].Playlist.append(
        GX.AnimatedUpdate(
          GX.duration(1.0),
          KML.Update(
            KML.targetHref(),
            KML.Change(
              KML.Placemark(
                KML.visibility(1),
                GX.balloonVisibility(1),
                targetId=fechas[i].isoformat()[0:13]
              )
            )
          )
        )
    )
    # creamos el mapa
    nombre = '{0}.png'.format(fechas[i].isoformat()[0:13])
    pinta_mapas(nombre, lon, lat, mslp[i,:,:], u[i,:,:], v[i,:,:], lon_hur[i], lat_hur[i])
    # añadimos el mapa al fichero kmz
    fich_kmz.write(nombre)
    # y eliminamos el fichero png
    os.remove(nombre)
    # añadimos la ruta de las imágenes en el
    # fichero kml (que estará contenido dentro
    # del fichero kmz.
    fich_kml.Document.Folder.append(
      KML.GroundOverlay(
        KML.name(nombre),
        KML.visibility(0),
        KML.Icon(KML.href(nombre)),
        KML.LatLonBox(KML.north(np.max(lat)),
                      KML.south(np.min(lat)),
                      KML.east(np.max(lon)),
                      KML.west(np.min(lon)),
        ),
        id = nombre
      )
    )
    fich_kml.Document[gxns+"Tour"].Playlist.append(GX.Wait(GX.duration(0.1)))
    # Aquí le indicamos a la reproducción que nos
    # muestre el 'mapa'. Visibility a 1 para
    # hacerlo visible.
    fich_kml.Document[gxns+"Tour"].Playlist.append(
        GX.AnimatedUpdate(
          GX.duration(1.0),
          KML.Update(
            KML.targetHref(),
            KML.Change(
              KML.GroundOverlay(
                KML.visibility(1),
                #GX.balloonVisibility(1),
                targetId=nombre
              )
            )
          )
        )
    )
    # Esta parte de código hará que desaparezca el 'globo'
    # con la información así no tendremos el 'globo'
    fich_kml.Document[gxns+"Tour"].Playlist.append(
        GX.AnimatedUpdate(
          GX.duration(0.1),
          KML.Update(
            KML.targetHref(),
            KML.Change(
              KML.Placemark(
                GX.balloonVisibility(0),
                targetId=fechas[i].isoformat()[0:13]
              )
            )
          )
        )
    )
    # Esta parte de código hará que desaparezca el 'mapa'
    fich_kml.Document[gxns+"Tour"].Playlist.append(
        GX.AnimatedUpdate(
          GX.duration(1.0),
          KML.Update(
            KML.targetHref(),
            KML.Change(
              KML.GroundOverlay(
                KML.visibility(0),
                targetId=nombre
              )
            )
          )
        )
    )

En el siguiente código lo que hacemos es guardar el fichero kml final. Si os acordáis, en algún momento hemos puesto algunas etiquetas HTML. Algunos símbolos especiales se han transformado por lo que los volvemos a dejar como estaban ya que si no google earth no entenderá ese kml correctamente.

# guardamos toda la información
# del kml en un fichero kml
outfile = file('HuracanIvan.kml','w')
salida = etree.tostring(fich_kml, pretty_print=True)
salida = salida.replace('<code>&lt;</code>', '<') salida = salida.replace('<code>&gt;</code>', '>')
outfile.write(salida)
outfile.close()
# guardamos el fichero kml dentro del
# fichero kmz
fich_kmz.write('HuracanIvan.kml')
# Eliminamos el fichero kml
os.remove('HuracanIvan.kml')
# Cerramos el fichero kmz
fich_kmz.close()

Por último, lanzamos google earth con el fichero que acabamos de crear. Las siguientes tres líneas funcionan en linux. Ahora mismo no tengo un windows cerca pero creo que deberéis cambiar googleearth por googleearth.exe y tener en cuenta los '' en lugar de '/' en la ruta al fichero.

ruta = os.getcwd()
rutafich = ruta + '/HuracanIvan.kmz'
os.system('googleearth {0}'.format(rutafich))

A la izquierda de la pantalla de google earth, en mis 'lugares', dentro de 'lugares temporales' podréis ver 'HuracanIvan.kmz'. Si lo desplegáis veréis 'Reprodúceme'. Si hacéis doble click sobre 'Reprodúceme' podréis ver algo parecido a lo siguiente:

[youtube=http://www.youtube.com/watch?v=TPKZo4ez1V0]

El resultado tampoco es maravilloso pero la idea no era el resultado, la idea era mostrar el uso de pykml e introducir un poco de kml/kmz a los que les pueda resultar útil.

Por último, como comenté en el anterior capítulo, os dejo un enlace con un uso creativo de google earth [1].

[1] - http://www.fosslc.org/drupal/content/pykml-python-kml-library

P.D.: Como puede que wordpress me haya 'estropeado' algo al hacer copy&paste os dejo el enlace al ipython notebook o lo podéis visualizar directamente en nbviewer.

Usando google earth con ayuda de python y pykml (I)

El tema de la visualización de datos es extremadamente importante puesto que te ayuda a comprender mejor al información que estás procesando. Si, además, esos datos se pueden ver de forma 'natural' en el medio donde suceden, donde se registran o donde se simulan pues mucho mejor. Por ejemplo, no tiene nada que ver el hecho de mirar por la tele un tornado o el verlo en vivo, de ver a una leona cazando en un documental o verlo en el cráter del Ngorongoro.

Los sistemas de información geográfica (SIG o GIS por sus sigles en inglés) permiten mostrar información georeferenciada. Google Earth permite ver información georeferenciada permitiendo verlo, además, de forma tridimensional por lo que, para determinados fenómenos, podemos ser conscientes de como se relacionan con su entorno. No es igual a verlo en vivo pero puede ser algo más cercano.

Hoy vamos a ver como usar pykml para crear ficheros kml (o kmz) para, posteriormente, ver el resultado en Google Earth. Los ficheros kml (pongo el enlace a la wikipedia en inglés puesto que en español la información es bastante limitada) no son más que ficheros basados en la sintáxis xml para mostrar información geográfica. Google compró la empresa que estaba detrás de ello y más tarde, en 2008, se hizo estándar de la Open Geospatial Consortium (OGC). Un fichero kmz no es más que un fichero kml comprimido y que puede contener otra información relacionada con el kml, por ejemplo, una imagen en formato jpeg georeferenciada. El fichero kmz es útil cuando queremos distribuir la información de forma autocontenida (aunque la imagen georeferenciada siempre podría estar en un disco duro local, en un servidor accesible a través de la red, etc).

Vamos a ver varios ejemplos. En el primero daremos una vuelta por el mundo visitando diferentes sitios interesantes y en el siguiente ejemplo veremos la evolución del huracán Ivan (2004) que llegó hasta huracán de categoría 5 en la escala Saffir-Simpson.

[Para esta entrada se ha usado pykml 0.1.0, lxml 2.3.2 (es un prerrequisito de pykml) en el notebook de IPython 0.13 corriendo todo sobre python 2.7.3]

import os
from lxml import etree

En un documento kml podemos tener 'carpetas' para organizar jerárquicamente nuestra información. Dentro de las 'carpetas' podemos meter 'marcas de posición', 'capas', otras 'carpetas', geometrías,... Aquí no voy a explicar mucho más puesto que nos desviamos de python. Si queréis algo concreto usad los comentarios o podéis ver toda la referencia de kml [aquí](https://developers.google.com/kml/documentation/kmlreference?hl=es).

También podemos hacer cosas muy potentes como mostrar líneas de tiempo y reproducir fenómenos de forma que podemos ver su evolución en el tiempo, mostrar modelos 3D,... Al final del segundo capítulo enlazaré una forma creativa de usar también google earth [1] pero no nos adelantemos.

Como he comentado, íbamos a ver dos ejemplos, Empecemos por el primero. Este primer ejemplo es un plagio remozado de un [ejemplo de la documentación](http://packages.python.org/pykml/examples/tour_examples.html) de pykml.

Primero importamos todo lo necesario de pykml para este ejemplo:

from pykml.factory import nsmap
from pykml.factory import KML_ElementMaker as KML
from pykml.factory import GX_ElementMaker as GX

Y ahora, con la ayuda de python, vamos a empezar a ir creando cosas en el documento kml y a explicarlas poco a poco:

# Desde la versión 5.0, google ha añadido una serie de extensiones
# a kml que lo hacen más potente e interactivo
# https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions
# Definimos una variable para el espacio de nombres de las
# extensiones de google que vamos a usar.
gxns = '{' + nsmap['gx'] + '}'
# Ponemos un listado de sitios que vamos a visitar
desc = """<![CDATA[
PyCon Spain 2013
<H3>Visitanos!!!</h3>
<img src="http://www.python.org/community/logos/python-logo-master-v3-TM.png"/>
]]>
"""
lugares = [
    {'name':"LHC",
     'desc':'En algun lugar ahi abajo esta el Gran Colisionador de Hadrones',
     'lon':6.053295,'lat':46.235339, 'tilt':45,
     'range':300, 'range2':10000000.0, 'vis':0,},
    {'name':"VLT",
     'desc':'ESO Very Large Telescope - VLT',
     'lon':-70.4042,'lat':-24.6275, 'tilt':45,
     'range':300, 'range2':10000000.0, 'vis':0,},
    {'name':"Rapa Nui",
     'desc':'Ahu Akivi',
     'lon':-109.395,'lat':-27.115, 'tilt':45,
     'range':50, 'range2':10000000.0, 'vis':0,},
    {'name':"Osaka Science Museum",
     'desc':'Museo de ciencia de Osaka',
     'lon':135.4915, 'lat':34.6912, 'tilt':45,
     'range':300, 'range2':10000000.0, 'vis':0,},
    {'name':"IMAX",
     'desc':'Cine IMAX de Glasgow',
     'lon':-4.2939, 'lat':55.8581, 'tilt':45,
     'range':300, 'range2':10000000.0, 'vis':0,},
    {'name':"PyConES",
     'desc':desc,
     'lon':-3.7, 'lat':40.4165, 'tilt':0,
     'range':15000, 'range2':10000.0, 'vis':1,},
]

En la anterior porción de código solo hemos detallado una serie de características de lugares que vamos a visitar. Ahora vamos a crear el esqueleto del documento que contendrá una carpeta (folder) donde se guardará información de los lugares que vamos a visitar (posición, descripción, nombre) y un elemento gx:Tour que hace uso de las extensiones de Google:

# start with a base KML tour and playlist
fich_kml = KML.kml(
    KML.Document(
      GX.Tour(
        KML.name(u"¡Reprodúceme!"),
        GX.Playlist(),
      ),
      KML.Folder(
        KML.name('Lugares visitados'),
        id='lugares',
      ),
    )
)

Y ahora vamos a hacer un bucle donde se van a ir añadiendo cosas tanto a la carpeta (folder) como al Tour. Voy explicándo lo que hace cada cosa mediante comentarios en el código:

# Hacemos un bucle para recorrer todos los lugares definidos en
# el bloque de código anterior
for lugar in lugares:
    # Primero volamos hasta cada lugar (transición entre lugares)
    # y nos quedamos observándolo desde el espacio.
    # Para ello usamos FlyTo donde especificamos
    # La duración en segundos, si queremos que vaya más rápido o más
    # lento, el modo de vuelo, que puede ser más suave (smooth) o
    # más rápido (bounce) y la posición.
    # La posición la determinan longitude, latitude y altitude,
    # tilt es el ángulo desde la vertical (en este caso lo vemos
    # verticalmente (0º)) y range es la distancia desde la posición
    # determinada por longitude, latitude y altitude, teniendo en
    # cuenta el ángulo tilt (en este caso hemos puesto 10.000 km)
    # Las unidades para latitude y longitude son grados, para altitude
    # y range son metros y para tilt son grados desde la vertical.
    # En el gráfico de este enlace se verá mejor:
    # https://developers.google.com/kml/documentation/images/lookAt.gif
    fich_kml.Document[gxns+"Tour"].Playlist.append(
      GX.FlyTo(
        GX.duration(1),
        GX.flyToMode("bounce"),
        KML.LookAt(
          KML.longitude(lugar['lon']),
          KML.latitude(lugar['lat']),
          KML.altitude(0),
          KML.heading(0),
          KML.tilt(45),
          KML.range(10000000.0),
          KML.altitudeMode("relativeToGround"),
        )
      ),
    )
    # Volamos desde el punto anterior (que estaba a 10.000 km
    # de altura) hasta el punto de vista desde donde queremos
    # visualizar cada lugar. Vamos a ver todos los lugares con un
    # ángulo (tilt) de 45º y a una altura de unos 300m de altura
    # Altura relativa al suelo, no al nivel del mar), excepto para
    # el caso de Rapa Nui, que bajaremos un poco más para que se vea mejor.
    fich_kml.Document[gxns+"Tour"].Playlist.append(
      GX.FlyTo(
        GX.duration(1),
        GX.flyToMode("bounce"),
        KML.LookAt(
          KML.longitude(lugar['lon']),
          KML.latitude(lugar['lat']),
          KML.altitude(0),
          KML.heading(0),
          KML.tilt(lugar['tilt']),
          KML.range(lugar['range']),
          KML.altitudeMode("relativeToGround"),
        )
      ),
    )
    # Damos una pequeña vuelta alrededor del lugar
    # en el que nos encontremos, para ello usamos
    # heading, que indica la dirección (azimut) de la
    # cámara
    for aspect in range(0,360,15):
        fich_kml.Document[gxns+"Tour"].Playlist.append(
          GX.FlyTo(
            GX.duration(0.1),
            GX.flyToMode("bounce"),
            KML.LookAt(
              KML.longitude(lugar['lon']),
              KML.latitude(lugar['lat']),
              KML.altitude(0),
              KML.heading(aspect),
              KML.tilt(lugar['tilt']),
              KML.range(lugar['range']),
              KML.altitudeMode("relativeToGround"),
            )
          )
        )
    fich_kml.Document[gxns+"Tour"].Playlist.append(GX.Wait(GX.duration(0.1)))
    # Añadimos información de cada lugar en la carpeta
    # como 'placemarks' o marcas de posición. En el nombre y en la descripción
    # del placemark se puede meter HTML y CSS, por si queréis ser un poco más
    # creativos y dejar la información de una forma visible y bonita.
    # extrude == 1 indica que el nombre del lugar se hará visible. Si queremos
    # hacerlo visible solo cuando lo seleccionemos pues le ponemos el valor 0
    fich_kml.Document.Folder.append(
      KML.Placemark(
        KML.name(lugar['name']),
        KML.description(
            "{name}<br/>{desc}".format(
                    name=lugar['name'],
                    desc=lugar['desc'],
            )
        ),
        KML.Point(
          KML.extrude(1),
          KML.altitudeMode("relativeToGround"),
          KML.coordinates("{lon},{lat},{alt}".format(
                  lon=lugar['lon'],
                  lat=lugar['lat'],
                  alt=0,
              )
          )
        ),
        id=lugar['name'].replace(' ','_')
      )
    )
    # Aquí le indicamos a la reproducción que en el
    # momento que haya dado la vuelta al lugar nos
    # muestre el 'globo' o 'viñeta' con la información del sitio,
    # nombre y descripción. Visibility a 1 para
    # hacerlo visible.
    fich_kml.Document[gxns+"Tour"].Playlist.append(
        GX.AnimatedUpdate(
          GX.duration(2.0),
          KML.Update(
            KML.targetHref(),
            KML.Change(
              KML.Placemark(
                KML.visibility(1),
                GX.balloonVisibility(1),
                targetId=lugar['name'].replace(' ','_')
              )
            )
          )
        )
    )
    fich_kml.Document[gxns+"Tour"].Playlist.append(GX.Wait(GX.duration(1.0)))
    # Esta parte de código hará que desaparezca el 'globo'
    # con la información así no tendremos el 'globo'
    fich_kml.Document[gxns+"Tour"].Playlist.append(
        GX.AnimatedUpdate(
          GX.duration(2.0),
          KML.Update(
            KML.targetHref(),
            KML.Change(
              KML.Placemark(
                GX.balloonVisibility(lugar['vis']),
                targetId=lugar['name'].replace(' ','_')
              )
            )
          )
        )
    )
    # Finalmente, esta parte nos separa de la tierra
    # para que tengamos una vista privilegiada del
    # planeta azul a 10.000 km de altura
    fich_kml.Document[gxns+"Tour"].Playlist.append(
      GX.FlyTo(
        GX.duration(5),
        GX.flyToMode("bounce"),
        KML.LookAt(
          KML.longitude(lugar['lon']),
          KML.latitude(lugar['lat']),
          KML.altitude(0),
          KML.heading(0),
          KML.tilt(0),
          KML.range(lugar['range2']),
          KML.altitudeMode("relativeToGround"),
        )
      ),
    )

En el siguiente código lo que hacemos es guardar el fichero kml final. Si os acordáis, en algún momento hemos puesto algunas etiquetas HTML. Algunos símbolos especiales se han transformado por lo que los volvemos a dejar como estaban ya que si no google earth no entenderá ese kml correctamente.

# output a KML file (named based on the Python script)
outfile = file('VolandoVoy.kml','w')
salida = etree.tostring(fich_kml, pretty_print=True)
salida = salida.replace('&lt;', '<') salida = salida.replace('&gt;', '>')
outfile.write(salida)
outfile.close()

Por último, lanzamos google earth con el fichero que acabamos de crear. Las siguientes tres líneas funcionan en linux. Ahora mismo no tengo un windows cerca pero creo que deberéis cambiar `googleearth` por `googleearth.exe` y tener en cuenta los '' en lugar de '/' en la ruta al fichero.

ruta = os.getcwd()
rutafich = ruta + '/VolandoVoy.kml'
os.system('googleearth {0}'.format(rutafich))

A la izquierda de la pantalla de google earth, en mis 'lugares', dentro de 'lugares temporales' podréis ver 'VolandoVoy.kml'. Si lo desplegáis veréis 'Reprodúceme'. Si hacéis doble click sobre 'Reprodúceme' podréis ver algo parecido a lo siguiente:

[youtube=http://www.youtube.com/watch?v=P-0mn_9ixjk]

No es un maravilla de vídeo pero solo es para que veáis un ejemplo en vivo del resultado.

AVISO: el otro día estuve probando todo el código anterior en un windows XP y me daba un error extraño. Cuando introducía la etiqueta `description` en el kml se escribía en mayúscula y google earth no entendía bien el kml y no me mostraba la información de la descripción... Estuve una hora larga con ello intentando adivinar qué es lo que estaba mal ya que el fichero kml se creaba sin problemas y google earth lo abría sin errores (pero no mostraba correctamente la descripción). Si a alguno le pasa, antes de guardar el fichero puede hacer salida = salida.replace('<Description>', '<description>').

Y esto es todo por hoy, en 24h tendréis la segunda parte disponible.

P.D.: Como puede que wordpress me haya 'estropeado' algo al hacer copy&paste os dejo el enlace al ipython notebook o lo podéis visualizar directamente en nbviewer.

Buscando esa playa en la isla a mediodía (usando Shapely)

Hay un cuento de Julio Cortazar titulado 'la isla a mediodía' en el cual el protagonista, asistente de vuelo, se queda maravillado con una isla en el mar Egeo que sobrevuela a menudo gracias a su trabajo.

Bien, imaginemos por un rato que somos el protagonista del cuento, Marini se llama, que en lugar de sobrevolar el Egeo estamos sobre el Mediterráneo, que la isla es Mallorca en lugar de ser una de las perlas griegas y que sobrevolamos la isla de noroeste a sudeste. Cada vez que pasamos nos quedamos fascinados con esa isla y en especial con una playa que vemos cuando llegamos a la isla.

¿Cómo podemos adivinar cual es esa cala?, ¿cómo podemos conocer el tamaño de la isla?, ¿cómo podemos saber el perímetro de la costa la cual nos gustaría recorrer en bicicleta en nuestras próximas vacaciones?

Una persona normal pensaría, me compro una guía lonely planet, busco en la wikipedia,... Pero ya sabéis que en pybonacci nos falta alguna tuerca y para ello vamos a hacer uso de Shapely, una librería que nos permite hacer una serie de operaciones típicas del mundo SIG (Sistemas de Información Geográfica, GIS en inglés) sobre espacios bidimensionales. Hay que dejar claro que se usan proyecciones y las operaciones que hacemos, si el espacio geográfico usado es grande, puede sufrir distorsiones importantes y los valores obtenidos podrían distar bastante de los teóricos (avisados estáis, leer más aquí).

[Para esta entrada usamos numpy 1.6, matplotlib 1.1 y  shapely 1.2 en python 2.7]

Shapely permite manejar polígonos, líneas, puntos,..., y permite hacer operaciones con ellos. Primero de todo vamos a descargar la línea de costa de la isla de Mallorca desde http://www.ngdc.noaa.gov/mgg/coast/. He usado datos del World Data Bank II a escala 1:2000000 y los he pasado a coordenadas UTM WGS84 y subido aquí para que los podáis descargar y reproducir los ejemplos de esta entrada). El fichero a descargar tiene extensión xls pero en realidad es un fichero csv (restricciones de wordpress.com, disculpen las molestias). Las unidades de los datos son metros.

Una vez que tenemos los datos vamos a comenzar los cálculos. Leemos los datos del fichero recién descargados:

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point, LineString, Polygon
data = np.loadtxt('Mallorca(UTM_WGS84).xls', delimiter=',')

Leídos los datos vamos a crear un objeto 'polígono' y a dibujarlo:

poligono = Polygon(data[:,0:2])
## Dibujamos el polígono accediendo a los datos (valores x e y del polígono
plt.ion()
plt.plot(list(poligono.exterior.xy)[0], list(poligono.exterior.xy)[1])

Calcular el área y la longitud del polígono es tan sencillo como hacer lo siguiente (los expresamos en km2 y en km, respectivamente):

print poligono.area/(1000.**2), poligono.length/1000.

(vuelvo a avisar, estos valores son una aproximación). Perfecto, ya tenemos nuestra isla soñada. Vamos a por la ruta de nuestro avión para saber por donde llega a la isla. Para ello, vamos a usar un objeto 'línea' y la vamos a dibujar sobre la isla:

linea = LineString([(440000,4410000),(560000,4350000)])
## Dibujamos la isla (polígono) y la ruta del avión (línea) accediendo a los valores xy
plt.plot(list(poligono.exterior.xy)[0], list(poligono.exterior.xy)[1])
plt.plot(list(linea.xy)[0], list(linea.xy)[1])

Genial, ahora, ¿cómo puedo saber la localización de la ruta del avión sobre la isla?, ¿cómo puedo saber cuantos kilómetros recorre el avión sobre la isla?. Nos ayudaremos de los objetos 'punto'.

## Para obtener la intersección entre la línea y el polígono
intersecciones = linea.intersection(poligono)
## Sabiendo donde intersecciona la ruta del avión podemos extraer
## el punto donde el avión llega a la isla y por donde la abandona
puntoNW = Point([intersecciones.xy[0][0], intersecciones.xy[1][0]])
puntoSE = Point([intersecciones.xy[0][1], intersecciones.xy[1][1]])
## Imprimimos en pantalla esos valores
print linea.intersection(poligono)
## Y, por último, calculamos la distancia entre los dos puntos
## (La distancia que el avión recorre sobre la isla en km)
print puntoNW.distance(puntoSE)/1000.

Bien, ya sé donde está esa playa que veo cuando estamos llegando a la isla (punto NW). Voy a dibujar la isla (polígono), la ruta del avión (línea) y los puntos donde el avión llega y sale de la isla (punto NW y punto SE, respectivamente):

plt.plot(list(poligono.exterior.xy)[0], list(poligono.exterior.xy)[1])
plt.plot(list(linea.xy)[0], list(linea.xy)[1])
plt.plot(intersecciones.xy[0][0], intersecciones.xy[1][0], 'rs')
plt.plot(intersecciones.xy[0][1], intersecciones.xy[1][1], 'ys')

Con este resultado:

Perfecto, ya conozco mi destino para las próximas vacaciones, sé donde está la playa de mis sueños, sé que la isla tiene un área considerable y sé que mi vuelta en bici alrededor de la isla la deberé hacer en varias etapas puesto que hay distancia suficiente para ello. Y todo ello gracias a python y a shapely!!!

Hasta la próxima.