Primeros pasos con Cartopy: paquete de python para crear mapas

Cartopy es un paquete de python diseñado para el procesamiento de datos geográficos con el fin de producir mapas y hacer análisis de datos geoespaciales. En este artículo veremos la forma en que podemos crear mapas utilizando Cartopy.

Cartopy hace uso de las potentes bibliotecas PROJ, NumPy y Shapely e incluye una interfaz programática construida sobre Matplotlib, empleándose para la creación de mapas con calidad de publicación.

Las características clave de cartopy son sus definiciones de proyección orientadas a objetos y su capacidad para transformar puntos, líneas, vectores, polígonos e imágenes entre esas proyecciones.

Instalación de Python, Anaconda y jupyter notebook

Una de las formas más rápidas y sencillas de instalar Python, es hacerlo mediante Anaconda.

Al lanzar el cuaderno de jupyter se abrirá una nueva pestaña en nuestro navegador. Una vez abierto, crearemos un nuevo cuaderno con python 3 haciendo clic en New → Python 3:

Mi primer mapa con cartopy

Una vez abierto el nuevo cuaderno, lo primero que tenemos que hacer es importar las librerías de cartopy y matplotlib:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

Después de importar las librerías podemos empezar a construir nuestro primer mapa. El ejemplo más sencillo es el que se muestra a continuación:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

ax = plt.axes(projection=ccrs.Mollweide())
ax.stock_img()
plt.show()

La clase principal para integrar cartopy en matplotlib es GeoAxes, que es una subclase de un matplotlib Axes normal. La clase GeoAxes agrega funcionalidad adicional a un eje que es específico para dibujar mapas. Resumiendo, podemos decir que cartopy utiliza matplotlib para crear unos “ejes” o entorno para el mapa, que se enriquezcan para poder utilizar coordenadas geográficas y no solo cartesianas.

A continuación estamos indicando que utilice la proyección Mollweide. Esta proyección es pseudo-cilíndrica y de igual área y se usa comúnmente para mapas del mundo que a veces se interrumpe con varios meridianos centrales. Empleamos stock_img() para agregar una imagen subyacente al mapa.

Trabajando con proyecciones y capas

Cartopy dispone de una lista de proyecciones muy amplia listas para utilizar.

En el ejemplo anterior se utiliza una proyección Goode Homosiline que es una proyección compuesta de áreas iguales que enfatiza las características de la tierra o el océano.

Del mismo modo, cartopy dispone de una serie de capas por defecto (features) tomadas de Natural Earth. Son las capas siguientes:

• BORDERS: Es un mapa de países del mundo.
• COASTLINE: representa las líneas de costa incluyendo las islas más importantes.
• LAKES: incluye lagos naturales y artificiales.
• LAND: es como COASTILE, pero en lugar de representar una línea muestra un polígono de los continentes e islas.
• OCEAN: polígonos de los océanos
• RIVERS: ríos del mundo.
• STATES: para los límites de los estados de EEUU.

Para utilizar estas capas importamos primero cartopy.feature como se indica en el siguiente ejemplo.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cf

ax = plt.axes(projection=ccrs.Mercator())

ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.LAND)

plt.show()

En la imagen anterior vemos un mapa al que hemos añadido las capas de línea de costa (COASTLINE) y LAND.

Configurar la vista del mapa

Para configurar los límites del mapa en cartopy podemos emplear la función: set_extend(x0 ,x1 ,y0 ,y1) en donde x0 , x1 se corresponden con la x mínima y la x máxima del área representada. De igual forma y0 e y1 se corresponden con la y mínima y la y máxima.

Siguiendo este criterio podemos definir un área del mapa que incluya la península ibérica:

ax.set_extent([-10, 4, 34, 45])

Vemos a continuación un ejemplo en que aplicamos estos límites.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cf

ax = plt.axes(projection=ccrs.Mercator())

ax.set_extent([-10, 4, 34, 45])
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.LAND)
ax.add_feature(cf.RIVERS)
ax.add_feature(cf.BORDERS)
plt.show()

Como habrás observado en la imagen anterior, ahora vemos unas líneas para el eje x y otras para el eje y, así como etiquetas indicando la latitud y la longitud. Esto es el resultado de la línea que añadimos en el código:

ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)

Todos los estilos de las líneas, etiquetas… y, en general de la vista del mapa, pueden ser personalizados. Por ejemplo si modificamos el valor de x_inline a True, los valores de longitud se inscriben dentro del mapa.

Trabajando con proyecciones y sistemas de coordenadas

Hasta ahora habíamos visto que Cartopy permite utilizar muchas proyecciones como PlateCarree, AlbersEqualArea, AzimuthalEquidistant, EquidistantConic, LambertConformal, LambertCylindrical, Mercator y muchas otras. Ahora lo que pretendemos es poder añadir datos sobre ese mapa.

El dato más sencillo que podemos emplear es añadir un punto a nuestro mapa. Sobre el ejemplo en que construíamos un mapa de la península ibérica vamos a añadir un punto:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cf

ax = plt.axes(projection=ccrs.Mercator())

ax.set_extent([-10, 4, 34, 45])
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.LAND)
ax.add_feature(cf.RIVERS)
ax.add_feature(cf.BORDERS)

xCoordenadas, yCoordenadas =-5.6641123, 40.9650084

plt.plot(xCoordenadas, yCoordenadas, marker= 'x',transform=ccrs.PlateCarree())
plt.show()

En el código anterior hemos creado un punto mediante sus coordenadas geográficas (xCoordenadas, yCoordenadas). Ese dato lo incorporamos al mapa mediante matplotlib, definiendo un símbolo para el marcador en forma de ‘x’.

xCoordenadas, yCoordenadas =-5.6641123, 40.9650084

Estas coordenadas se corresponden con EPSG 4326 que utiliza una proyección PlateCarree, sin embargo habíamos creado el mapa en una proyección Mercator. Para que el punto se pueda ver en la posición correcta necesitamos indicarle el sistema de coordenadas en el que se encuentra. Es lo que hacemos utilizando: transform.

Cartopy maneja dos conceptos fundamentales:

• Por un lado esta projection que es el que se encarga de crear el gráfico del mapa.
• Y por otro transform que es el que se encarga de indicar el sistema de coordenadas en que se encuentran los datos. Si no utilizamos transform, cartopy interpreta que los datos están en el mismo sistema de coordenadas que la proyección.

plt.plot(xCoordenadas, yCoordenadas, marker= 'x',transform=ccrs.PlateCarree())

Importar datos geográficos en formato shp

Cartopy nos permite utilizar fuentes de datos geográficos como shapefile. Cartopy proporciona un lector de archivos de forma orientado a objetos basado en el módulo pyshp. El envoltorio de pyshp de Cartopy tiene la ventaja de ser Python puro y, por lo tanto, es fácil de instalar y extremadamente portátil.

Aunque en este ejemplo estamos utilizando Cartopy para la entrada de datos, es recomendable  el empleo de Fiona y GeoPandas.

Veamos un ejemplo en el que utilizamos un shapefile proporcionado por Natural earth:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

from cartopy.feature import ShapelyFeature

shpfilename = shpreader.natural_earth(resolution='110m',
                                      category='cultural',
                                      name='admin_0_countries')
ax = plt.axes(projection=ccrs.PlateCarree())

reader = shpreader.Reader(shpfilename)
shape_feature = ShapelyFeature(reader.geometries(),ccrs.PlateCarree(),facecolor='none')

ax.add_feature(shape_feature)
plt.show()

El ejemplo anterior es un caso especial, ya que Cartopy proporciona una interfaz para acceder a los datos de el sitio web NaturalEarthData.

Mediante shapereader.natural_earth obtenemos la ruta al archivo de natural_earth y además lo descarga y descomprime.

shapreader.BasicReader o su alias Reader proporciona una interfaz para acceder al contenido del shapefile:

En un caso más general, en lugar de enlazar al recurso de Natural Earth podemos utilizar un archivo shp. Puedes descargar el mismo archivo que hemos utilizado antes y ahora llamarlo desde el script de cartopy:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

from cartopy.feature import ShapelyFeature

shpfilename = 'ne_110m_admin_0_countries.shp'
ax = plt.axes(projection=ccrs.PlateCarree())

reader = shpreader.Reader(shpfilename)
shape_feature = ShapelyFeature(reader.geometries(),ccrs.PlateCarree(),facecolor='none')

ax.add_feature(shape_feature)
plt.show()

Como puedes observar, ahora estamos utilizando el archivo shp descargado, siendo el resultado el mismo que en el mapa anterior.

En conclusión podemos decir que cartopy es una librería de python fácil de utilizar y que nos permite crear mapas a partir de nuestros datos geográficos, permitiendo el procesado de los datos y la personalización de los estilos del mapa.