In [1]:
name = "2019-05-30-cartopy-map"
title = "Pleasing maps with cartopy"
tags = "matplotlib, cartopy, dataviz, maps"
author = "Callum Rollo"
In [2]:
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
In [3]:
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
In [4]:
import cartopy
In addition, we import cartopy's coordinate reference system submodule:
In [5]:
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
In [6]:
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.
In [7]:
# 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
In [8]:
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
In [9]:
# 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]
In [10]:
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
In [11]:
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
In [12]:
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
In [13]:
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
In [14]:
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);
In [15]:
HTML(html)
Out[15]: