En ocasiones necesitaremos representar datos sobre un mapa. En este tipo de casos basemap es una buena alternativa dentro del ecosistema Python, pero pronto será sustituida por cartopy. Aunque seguirá teniendo mantenimiento hasta 2020 y cartopy todavía no ha incorporado todas las características de basemap, miraremos al futuro y haremos nuestros primeros ejemplos con la nueva biblioteca. Si aun así te interesa basemap puedes ver esta entrada en el blog de Jake Vanderplas o este notebook de su libro sobre data science.
En primer lugar, como siempre importamos la librería y el resto de cosas que necesitaremos:
In [1]:
%matplotlib notebook
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
Veremos en un primer ejemplo cómo crear un mapa con una proyección y añadiremos la información que nos interese:
In [2]:
# Inicializamos una figura con el tamaño que necesitemos
# si no la queremos por defecto
plt.figure()
# Creamos unos ejes con la proyección que queramos
# por ejemplo, Mercator
ax = plt.axes(projection=ccrs.Mercator())
# Y lo que queremos representar en el mapa
# Tierra
ax.add_feature(cfeature.LAND)
# Océanos
ax.add_feature(cfeature.OCEAN)
# Líneas de costa (podemos modificar el color)
ax.add_feature(cfeature.COASTLINE, edgecolor=(0.3, 0.3, 0.3))
# Fronteras
ax.add_feature(cfeature.BORDERS, edgecolor=(0.4, 0.4, 0.4))
# Ríos y lagos
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.LAKES)
# Por último, podemos pintar el grid, si nos interesa
ax.gridlines()
Out[2]:
Veremos ahora otro ejemplo usando una proyección distinta y colorearemos el mapa de forma diferente:
In [3]:
# Inicializamos una figura con el tamaño que necesitemos
# si no la queremos por defecto
plt.figure()
# Elegimos la proyección InterruptedGoodeHomolosine
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())
# Y lo que queremos representar en el mapa
ax.stock_img()
ax.gridlines()
Out[3]:
Puede interesarnos poner etiquetas a los ejes. Podemos utilizar entonces las herramientas dentro de : cartopy.mpl.gridliner
In [4]:
# Importamos los formatos de ejes para latitud y longitud
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
In [5]:
plt.figure()
# Elegimos la proyección PlateCarree
ax = plt.axes(projection=ccrs.PlateCarree())
# Y lo que queremos representar en el mapa
# Tierra
ax.add_feature(cfeature.LAND)
# Océanos
ax.add_feature(cfeature.OCEAN)
# Líneas de costa (podemos modificar el color)
ax.add_feature(cfeature.COASTLINE, edgecolor=(0.3, 0.3, 0.3))
# Fronteras
ax.add_feature(cfeature.BORDERS, edgecolor=(0.4, 0.4, 0.4))
# Ríos y lagos
ax.add_feature(cfeature.RIVERS)
ax.add_feature(cfeature.LAKES)
# Dentro de los ejes seleccionamos las lineas del grid y
# activamos la opción de mostrar etiquetas
gl = ax.gridlines(draw_labels=True)
# Sobre las líneas del grid, ajustamos el formato en x e y
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
Para las ocasiones en las que no queremos mostrar un mapa entero, sino que sólo necesitamos representar una determinada localización, vale todo lo anterior, simplemente tendremos espedificar el área a mostrar con el método set_extent
y tomar alguna precaución...
In [6]:
plt.figure()
# Elegimos la proyección
ax = plt.axes(projection=ccrs.PlateCarree())
# Fijar el punto y la extensión del mapa que queremos ver
ax.set_extent([3, -10, 35, 45], ccrs.Geodetic())
# Y lo que queremos representar en el mapa
ax.add_feature(cfeature.COASTLINE,
edgecolor=(0.3, 0.3, 0.3),
facecolor=cfeature.COLORS['land']
)
gl = ax.gridlines(draw_labels=True)
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
Como se ve en la figura enterior, la representación que obtenemos es demasiado burda. Esto se debe a que los datos por defecto se encuentran descargados a una escala poco detallada.
Cartopy permite acceder a datos propios almacenados en nuestro ordenador o descargarlos de algunas bases de datos reconocidas. En este caso accederemos a NaturalEarthFeature, que es la que hemos utilizado por defecto hasta ahora sin saberlo.
In [7]:
# Importando Natural Earth Feature
from cartopy.feature import NaturalEarthFeature
In [8]:
fig = plt.figure()
# Elegimos la proyección
ax = plt.axes(projection=ccrs.PlateCarree())
# Fijar el punto y la extensión del mapa que queremos ver
ax.set_extent([3, -10, 35, 45], ccrs.Geodetic())
# Y lo que queremos representar en el mapa
# Hasta ahora utilizábamos:
# ax.add_feature(cfeature.COASTLINE,
# edgecolor=(0.3, 0.3, 0.3),
# facecolor=cfeature.COLORS['land']
# )
# Pero ahora descargaremos primero la característica que
# queremos representar:
coastline = NaturalEarthFeature(category='physical', name='coastline', scale='10m')
# Y después la añadiremso al mapa con las propiedades que creamos convenientes.
ax.add_feature(coastline,
facecolor=cfeature.COLORS['land'],
edgecolor='k',
alpha=0.5)
gl = ax.gridlines(draw_labels=True)
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
Desde Naturale Earth Feature, no sólo podemos descargar caraterísticas físicas, sino que también podemos acceder a datasets demográficos
In [9]:
import pandas as pd
Habitualmente, no queremos sólo pintar un mapa sino que queremos representar datos sobre él. Estos datos pueden venir del dataset anterior de cualquier otra fuente.
En este ejemplo representaremos los datos de impactos de meteoritos que han caído en la tierra que están recogidos en el dataset: www.kaggle.como/nasa/meteorite-landings. En el link se pude encontrar toda la información.
In [10]:
# Leeremos el csv que ya tenemos descargado utilizando pandas
import pandas as pd
In [11]:
meteorite = pd.read_csv("../data/meteorite-landings.csv")
meteorite.head()
Out[11]:
In [12]:
import numpy as np
In [13]:
# Cremos un mapa sobre el que representar los datos:
fig = plt.figure(figsize=(8,4))
# Elegimos la proyección PlateCarree
ax = plt.axes(projection=ccrs.PlateCarree())
# Y representamos las líneas de costa
coastline = NaturalEarthFeature(category='physical', name='coastline', scale='50m')
ax.add_feature(coastline, facecolor=cfeature.COLORS['land'], edgecolor='k', alpha=0.5)
# Ahora podemos añadir sobre ese mapa los datos con un scatter
cl = ax.scatter(meteorite.reclong,
meteorite.reclat,
c=np.log10(meteorite.mass),
alpha=0.9,
cmap=plt.cm.magma_r)
plt.colorbar(cl, shrink=0.400)