name = "2019-05-30-cartopy-map"
title = "Pleasing maps with cartopy"
tags = "matplotlib, cartopy, dataviz, maps"
author = "Callum Rollo"
from nb_tools import connect_notebook_to_post
from IPython.core.display import HTML
html = connect_notebook_to_post(name, title, tags, author)
To make a pretty, publication grade map for your study area look no further than cartopy.
In this tutorial we will walk through generating a basemap with:
This code can be generalised to any region you wish to map
First we import some modules for manipulating and plotting data
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import xarray as xr
from pathlib import Path
Then we import cartopy itself
import cartopy
In addition, we import cartopy's coordinate reference system submodule:
import cartopy.crs as ccrs
A few other modules and functions which we will use later to add cool stuff to our plots. Also updating font sizes for improved readability
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
plt.rcParams.update({"font.size": 20})
SMALL_SIZE = 22
MEDIUM_SIZE = 22
LARGE_SIZE = 26
plt.rc("font", size=SMALL_SIZE)
plt.rc("xtick", labelsize=SMALL_SIZE)
plt.rc("ytick", labelsize=SMALL_SIZE)
plt.rc("axes", titlesize=SMALL_SIZE)
plt.rc("legend", fontsize=SMALL_SIZE)
To save space and time I have subset the bathymetry plotted in this example. If you wish to map a different area you will need to download the GEBCO topography data found here.
You can find a notebook intro to using xarray for netcdf here on the UEA python website. Or go to Callum's github for a worked example using GEBCO data.
# Open prepared bathymetry dataset using pathlib to sepcify the relative path
bathy_file_path = Path('../data/bathy.nc')
bathy_ds = xr.open_dataset(bathy_file_path)
bathy_lon, bathy_lat, bathy_h = bathy_ds.bathymetry.longitude, bathy_ds.bathymetry.latitude, bathy_ds.bathymetry.values
We're just interested in bathy here, so set any height values greater than 0 to to 0 and set contour levels to plot later
bathy_h[bathy_h > 0] = 0
bathy_conts = np.arange(-9000, 500, 500)
Here we load some scatter data from a two column csv for plotting later
# Load some scatter data of smaple locations near South Georgia
data = pd.read_csv("../data/scatter_coords.csv")
lons = data.Longitude.values
lats = data.Latitude.values
# Subset of sampling locations
sample_lon = lons[[0, 2, 7]]
sample_lat = lats[[0, 2, 7]]
Now to make the map itself. First we define our coordinate system. Here we are using a Plate Carrée projection, which is one of equidistant cylindrical projections.
A full list of Cartopy projections is available at http://scitools.org.uk/cartopy/docs/latest/crs/projections.html.
Then we create figure and axes instances and set the plotting extent in degrees [West, East, South, North]
coord = ccrs.PlateCarree()
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111, projection=coord)
ax.set_extent([-42, -23, -60, -50], crs=coord);
Now we contour the bathymetry data
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111, projection=coord)
ax.set_extent([-42, -23, -60, -50], crs=coord)
bathy = ax.contourf(bathy_lon, bathy_lat, bathy_h, bathy_conts, transform=coord, cmap="Blues_r")
A good start. To make it more map like we add gridlines, formatted labels and a colorbar
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111, projection=coord)
ax.set_extent([-42, -23, -60, -50], crs=coord)
bathy = ax.contourf(bathy_lon, bathy_lat, bathy_h, bathy_conts, transform=coord, cmap="Blues_r")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color="k", alpha=0.5, linestyle="--")
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylines = True
gl.xlines = True
fig.colorbar(bathy, ax=ax, orientation="horizontal", label="Bathymetry (m)", shrink=0.7, pad=0.08, aspect=40);
Now to add a few more features. First coastlines from cartopy's natural features toolbox. Then scatters of the samples we imported earlier
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111, projection=coord)
ax.set_extent([-42, -23, -60, -50], crs=coord)
bathy = ax.contourf(bathy_lon, bathy_lat, bathy_h, bathy_conts, transform=coord, cmap="Blues_r")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color="k", alpha=0.5, linestyle="--")
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylines = True
gl.xlines = True
fig.colorbar(bathy, ax=ax, orientation="horizontal", label="Bathymetry (m)", shrink=0.7, pad=0.08, aspect=40)
feature = cartopy.feature.NaturalEarthFeature(
name="coastline", category="physical", scale="50m", edgecolor="0.5", facecolor="0.8"
)
ax.add_feature(feature)
ax.scatter(lons, lats, zorder=5, color="red", label="Samples collected")
ax.scatter(sample_lon, sample_lat, zorder=10, color="k", marker="D", s=50, label="Samples sequenced");
To finish off the map we add a legend for the scatter plot, an inset map showing the area at a larger scale and some text identifying the islands
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(111, projection=coord)
ax.set_extent([-42, -23, -60, -50], crs=coord)
bathy = ax.contourf(bathy_lon, bathy_lat, bathy_h, bathy_conts, transform=coord, cmap="Blues_r")
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color="k", alpha=0.5, linestyle="--")
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylines = True
gl.xlines = True
fig.colorbar(bathy, ax=ax, orientation="horizontal", label="Bathymetry (m)", shrink=0.7, pad=0.08, aspect=40)
ax.add_feature(feature)
ax.scatter(lons, lats, zorder=5, color="red", label="Samples collected")
ax.scatter(sample_lon, sample_lat, zorder=10, color="k", marker="D", s=50, label="Samples sequenced")
fig.legend(bbox_to_anchor=(0.12, 0.2), loc="lower left")
tr2 = ccrs.Stereographic(central_latitude=-55, central_longitude=-35)
sub_ax = plt.axes(
[0.63, 0.65, 0.2, 0.2], projection=ccrs.Stereographic(central_latitude=-55, central_longitude=-35)
)
sub_ax.set_extent([-70, -15, -75, 10])
x_co = [-42, -42, -23, -23, -42]
y_co = [-60, -50, -50, -60, -60]
sub_ax.add_feature(feature)
sub_ax.plot(x_co, y_co, transform=coord, zorder=10, color="red")
ax.text(-38.5, -54.9, "South\nGeorgia", fontsize=14)
ax.text(-26.8, -58.2, "South\nSandwich\nIslands", fontsize=14);
HTML(html)
This post was written as an IPython (Jupyter) notebook. You can view or download it using nbviewer.