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)
Out[13]:
Casi cualquier representación de las que hemos cisto anteriormente con matploltib es posible.
As an example of visualizing some more continuous geographic data, let's consider the "polar vortex" that hit the eastern half of the United States in January of 2014. A great source for any sort of climatic data is NASA's Goddard Institute for Space Studies. Here we'll use the GIS 250 temperature data, which we can download using shell commands (these commands may have to be modified on Windows machines). The data used here was downloaded on 6/12/2016, and the file size is approximately 9MB:
The data comes in NetCDF format, which can be read in Python by the netCDF4 library.
You can install this library as shown here
$ conda install netcdf4We read the data as follows:
In [14]:
# preserve
from netCDF4 import Dataset
from netCDF4 import date2index
from datetime import datetime
In [15]:
# preserve
data = Dataset('../data/gistemp250.nc')
The file contains many global temperature readings on a variety of dates; we need to select the index of the date we're interested in—in this case, January 15, 2014:
In [16]:
# preserve
timeindex = date2index(datetime(2014, 1, 15),
data.variables['time'])
Now we can load the latitude and longitude data, as well as the temperature anomaly for this index:
In [17]:
# preserve
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lon, lat = np.meshgrid(lon, lat)
temp_anomaly = data.variables['tempanomaly'][timeindex]
Finally, we'll use the pcolormesh() method to draw a color mesh of the data.
We'll look at North America, and use a shaded relief map in the background.
Note that for this data we specifically chose a divergent colormap, which has a neutral color at zero and two contrasting colors at negative and positive values.
We'll also lightly draw the coastlines over the colors for reference:
In [18]:
# preserve
fig = plt.figure(figsize=(8,4))
# Elegimos la proyección
ax = plt.axes(projection=ccrs.PlateCarree())
# Y lo que queremos representar en el mapa
coastline = NaturalEarthFeature(category='physical', name='coastline', scale='50m')
# ax.add_feature(land, color=cfeature.COLORS['land'])
ax.add_feature(coastline, facecolor=cfeature.COLORS['land'], edgecolor='k', alpha=0.5)
ax.pcolormesh(lon, lat, temp_anomaly, cmap='RdBu_r')
Out[18]:
The data paints a picture of the localized, extreme temperature anomalies that happened during that month. The eastern half of the United States was much colder than normal, while the western half and Alaska were much warmer. Regions with no recorded temperature show the map background.
En este notebook hemos aprendido:
cartopy es capaz de muchas más cosas, por ejemplo: