Basemap y Google Geocode para representar puntos sobre un mapa

Como siempre, en GitHub podréis encontrar el notebook y los ficheros adicionales que forman parte de éste artículo.

Las gráficas pueden proporcionar mucha información de un sólo vistazo, pero no siempre son el mejor método de representación. A veces es necesario dar un paso más; y ese es precisamente el caso que vamos a tratar en éste artículo.

Representar valores sobre un mapa geográfico nos permite ubicar la información sobre el terreno. Lo datos a representar pueden ser desde contornos de temperatura hasta vectores de velocidad del viento pasando por la identificación de puntos geográficos. Para facilitarnos todo ese trabajo —en Python— disponemos de una serie de librerías. PyNGL, CDAT o Basemap —del que Kiko ya ha hablado en ésta entrada— son librerías que nacieron para satisfacer las necesidades de ciertos colectivos de científicos, como meteorólogos u oceanógrafos.

En éste notebook vamos a utilizar:

Datos

Lo primero es tener claro qué es lo que vamos a representar. Eso nos permitirá definir el tipo de representación a utilizar, así como las características del mapa —su proyección y límites geográficos.

Circuitos de Formula 1

La Formula 1 —categoría reina del automovilísmo— comenzó sus andanzas en el año 1950, con una primera prueba en Silverstone. Desde entonces, y aunque su foco de actividad se encuentra principalmente en Europa, ha ido expandiéndose para celebrar Grandes Premios por todo el mundo.

En Wikipedia podemos encontrar una simple tabla con una lista de todos los circuitos que alguna vez han albergado un Gran Premio. La copiamos y generamos un fichero CSV o XLS en Excel o en un editor de texto como Notepad++.

import pandas as pd
data = pd.DataFrame.from_csv('F1-circuits.csv', header=0, sep=';', index_col=None, parse_dates=False, encoding='latin-1')
data.head(5)
Circuit Type Direction Location Current Length Grands Prix Season(s) Grands Prix held
0 Adelaide Street Circuit Street Clockwise Adelaide, Australia 3.780 km (2.349 mi) Australian Grand Prix 1985-1995 11
1 Ain-Diab Circuit Road Clockwise Casablanca, Morocco 7.618 km (4.734 mi) Moroccan Grand Prix 1958 1
2 Aintree Road Clockwise Liverpool, United Kingdom 4.828 km (3.000 mi) British Grand Prix 1955, 1957, 1959, 1961-1962 5
3 Albert Park Street Clockwise Melbourne, Australia 5.303 km (3.295 mi) Australian Grand Prix 1996-2014 19
4 AVUS Street Anti-clockwise Berlin, Germany 8.300 km (5.157 mi) German Grand Prix 1959 1

Google Geocode

Ya tenemos la tabla cargada con datos como el nombre del circuito, tipo y localización, así como el número de Grandes Premios albergados. Pero entre toda esa información no hay nada que le diga a Python dónde colocar el circuito en un mapa —aunque nosotros si sepamos ubicar Casablanca o Melbourne—. Aquí es donde entra en juego el API de codificación geográfica de Google.

La codificación geográfica es el proceso de transformar direcciones (como "1600 Amphitheatre Parkway, Mountain View, CA") en coordenadas geográficas (como 37.423021 de latitud y -122.083739 de longitud), que se pueden utilizar para colocar marcadores o situar el mapa. El API de codificación geográfica de Google proporciona una forma directa de acceder a un geocoder mediante solicitudes HTTP.

El uso del API de codificación geográfica de Google está sujeto a un límite de 2.500 solicitudes de codificación geográfica al día, más que suficientes para ubicar los cerca de 70 circuitos que han albergado alguna vez en su historia un Gran Premio de Fórmula 1.

Como bien indica la documentación, al geocoder se accede mediante solicitudes HTTP, y para ello nada mejor que Requests. Todas las consultas se realizan a una dirección HTTP http://maps.googleapis.com/maps/api/geocode/json a la que se añaden una serie de parámetros obligatorios:

  • address: es la dirección que quieres codificar de forma geográfica.
  • sensor: indica si la solicitud de codificación geográfica procede de un dispositivo con un sensor de ubicación. Este valor debe ser true o false.

Hay, además, una serie de parámetros opcionales, pero no nos harán falta.

La consulta —requests.get(url, params)— devuelve una respuesta en formato json donde se incluyen las coordenadas geográficas que buscamos. El formato json, Python lo interpreta como un conjunto de diccionarios y arrays, para lo que indicaremos los key y los índices hasta llegar al punto donde se encuentra la información que buscamos. En éste caso: ['results'][0]['geometry']['location'].

import requests
_GEOCODE_QUERY_URL = 'http://maps.googleapis.com/maps/api/geocode/json'
def geocode(address, sensor='false'):
    """
    Given a string 'address', return a dictionary of information about
    that location, including its latitude and longitude.
    """
    params = dict(address=address, sensor=sensor)
    response = requests.get(url=_GEOCODE_QUERY_URL, params=params)
    return response.json()
def address_to_latlng(address):
    """
    Given a string 'address', return a '(latitude, longitude)' pair.
    """
    location = geocode(address)['results'][0]['geometry']['location']
    return tuple(location.values())

Basemap

Basemap es un toolkit de matplotlib que nos facilita la tarea de representar información 2D sobre mapas. Ésta información pueden ser contornos, vectores o puntos entre otros como se puede ver en los ejemplos.

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt

Lo primero, vamos a definir es el tipo de proyección a emplear. Hay un montón de ellas descritas en la Wikipedia, y en Basemap disponemos de 24 entre las que escoger. Para gustos, proyecciones.

En este caso hemos optado por Eckert IV, una proyección pseudocilíndrica, para representar el mapamundi y la Albers Equal Area projection para Europa. Dibujaremos las líneas de costa, las fronteras entre países, los paralelos y meridianos y le daremos un toque de color a los continentes. Basemap, además de disponer de una base de datos con información para representar líneas costeras y fronteras políticas, permite utilizar una imagen como fondo para el mapa. Entre las opciones que ofrece Basemap podemos encontrar el Blue Marble de la NASA —m.bluemarble()—. Aquí hemos optado por una imagen shaded relief de tonos claros con m.shadedrelief().

def basic_world_map(ax=None, region='world'):
    if region=='world':
        m = Basemap(resolution='i',projection='eck4',
                    lat_0=0,lon_0=0)
        # draw parallels and meridians.
        m.drawparallels(np.arange(-90.,91.,30.))
        m.drawmeridians(np.arange(-180.,181.,30.))
    elif region=='europe':
        m = Basemap(width=4000000,height=4000000,
                    resolution='l',projection='aea',\
                    lat_1=40.,lat_2=60,lon_0=10,lat_0=50)
        # draw parallels and meridians.
        m.drawparallels(np.arange(-90.,91.,10.))
        m.drawmeridians(np.arange(-180.,181.,10.))
        m.shadedrelief(scale=0.5)
    m.drawcoastlines()
    m.drawcountries()
    m.fillcontinents(color='coral', alpha=0.3)
    return m

Creamos un subplot y le asignamos un título a la figura. En esa figura vamos a representar las localizaciones de los circuitos con puntos con un área que vendrá determinada por el número de carreras disputadas —tanto mayor será el círculo cuantas más carreras se hayan disputado.

Para que los circulos no sean demasiado grandes —en Monza se han celebrado 64 Grandes Premios— limitaremos el radio del círculo a entre 3 y 20 puntos. Le damos a los cículos algo de transparencia con alpha=0.7, añadimos una nota de texto y guardamos la figura.

maximum = data['Grands Prix held'].max()
minimum = data['Grands Prix held'].min()
f, ax = plt.subplots(figsize=(20, 8))
ax.set_title('Formula 1 Grand Prix Circuits since 1950\n(Radius by number of races held)')
m = basic_world_map(ax)
for cir, loc, num in zip(data['Circuit'].values, data['Location'].values, data['Grands Prix held'].values):
    lat, lng = address_to_latlng(cir + ', ' + loc)
    x, y = m(lat, lng)
    m.scatter(x, y, s=np.pi * (3 + (num-minimum)/(maximum-minimum)*17)**2, marker='o', c='red', alpha=0.7)
ax.annotate(u'\N{COPYRIGHT SIGN} 2014, Pablo Fernandez', (0, 0))
f.savefig('f1-circuits.png', dpi=72, transparent=False, bbox_inches='tight')

f1-circuits

Podemos ver una gran concentración de Grandes Premios en Europa, continente que vio nacer a la Fórmula 1 y base de operaciones de la mayoría de equipos que compiten en ella. Si centramos la imagen sobre europa, a la cual hemos añadido un fondo, podremos ver con mayor claridad la distribución de las que han sido sedes de algún Gran Premio por el viejo continente.

f, ax = plt.subplots(figsize=(20, 8))
ax.set_title('Formula 1 Grand Prix Circuits in Europe since 1950\n(Radius by number of races held)')
m = basic_world_map(ax, 'europe')
for cir, loc, num in zip(data['Circuit'].values, data['Location'].values, data['Grands Prix held'].values):
    lat, lng = address_to_latlng(cir + ', ' + loc)
    x, y = m(lat, lng)
    m.scatter(x, y, s=np.pi * (3 + (num-minimum)/(maximum-minimum)*17)**2, marker='o', c='red', alpha=0.7)
ax.annotate(u'\N{COPYRIGHT SIGN} 2014, Pablo Fernandez', (100000, 100000))
f.savefig('f1-circuits-europe.png', dpi=72, transparent=False, bbox_inches='tight')

f1-circuits-europe

Manual de introducción a matplotlib.pyplot (VII): Tipos de gráfico (IV)

Esto pretende ser un tutorial del módulo pyplot de la librería matplotlib. El tutorial lo dividiremos de la siguiente forma (que podrá ir cambiando a medida que vayamos avanzando).

  1. Primeros pasos
  2. Creando ventanas, manejando ventanas y configurando la sesión
  3. Configuración del gráfico
  4. Tipos de gráfico I
  5. Tipos de gráfico II
  6. Tipos de gráfico III
  7. Tipos de gráfico IV
  8. Texto y anotaciones (arrow, annotate, table, text...)
  9. Herramientas estadísticas (acorr, cohere, csd, psd, specgram, spy, xcorr, ...)
  10. Eventos e interactividad (connect, disconnect, ginput, waitforbuttonpress...)
  11. Miscelánea

[Para este tutorial se ha usado python 2.7.1, ipython 0.11, numpy 1.6.1, matplotlib 1.1.0, netcdf4-python 0.9.9 y Basemap 1.0.2]

[DISCLAIMER: Muchos de los gráficos que vamos a representar no tienen ningún sentido físico y los resultados solo pretenden mostrar el uso de la librería].

En todo momento supondremos que se ha iniciado la sesión y se ha hecho

import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
from mpl_toolkits.basemap import Basemap as Bm

Hasta ahora hemos visto como configurar las ventanas, manejo de las mismas, definir áreas de gráfico, algunos tipos de gráficos... Ahora vamos a ver un último capítulo sobre tipos de gráficos. En esta última entrada sobre los tipos de gráfico hemos metido gráficos que quizá no estén muy relacionados entre sí por lo que quizá este capítulo puede parecer un poco cajón desastre.

Continue reading

Ejemplo de uso de Basemap y NetCDF4

Continuando lo que enseñó Juanlu en la anterior entrada vamos a mostrar líneas de nivel y temperatura del aire en la superficie, en este caso la presión al nivel del mar del día 01 de enero de 2012 a las 00.00 UTC según los datos extraídos del reanálisis NCEP/NCAR, sobre un mapa con la ayuda de la librería Basemap.

Como los datos del reanálisis NCEP/NCAR vienen en formato netCDF usaremos la librería netcdf4-python. El formato netCDF es un estándar abierto y es ampliamente usado en temas de ciencias de la tierra, atmósfera, climatología, meteorología,... No es estrictamente necesario usar netcdf4-python para acceder a ficheros netCDF puesto que desde scipy tenéis esta funcionalidad. Pero bueno, yo uso esta por una serie de ventajas que veremos otro día.

En la presente entrada se ha usado python 2.7.2, numpy 1.6.1, matplotlib 1.1.0, netCDF4 0.9.7 y Basemap 1.0.2.

Primero de todo vamos a importar todo lo que necesitamos:

## Importamos las librerías que nos hacen falta
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits import basemap as bm

Los ficheros netCDF de presión al nivel del mar y de Temperatura del aire de la superficie los podéis descargar de aquí y aquí, respectivamente. Veréis un enlace que pone 'FTP a copy of the file', lo pincháis y guardáis en el mismo sitio donde tengáis el script que estamos haciendo en la presente entrada.

Una vez que tenemos los ficheros los podemos abrir usando la librería netCDF4-python:

## Abrimos los ficheros de datos,
## el nombre de los ficheros lo tendréis que cambiar
## con el nombre de los ficheros que os habéis descargado
slp = nc.Dataset('X83.34.8.250.104.4.18.19.nc') #slp por 'sea level pressure'
tsfc = nc.Dataset('X83.34.8.250.104.4.15.31.nc') #tsfc 'por temperature at surface'

Continue reading