{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": [],
"toc_visible": true,
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\")
"
]
},
{
"cell_type": "markdown",
"source": [
"# Visualization with Matplotlib and Cartopy"
],
"metadata": {
"id": "6wgRMSYzJ2U8"
}
},
{
"cell_type": "markdown",
"source": [
"## Matplotlib"
],
"metadata": {
"id": "9hrcm3GDKY21"
}
},
{
"cell_type": "markdown",
"source": [
""
],
"metadata": {
"id": "p2PoyWtb69nv"
}
},
{
"cell_type": "markdown",
"source": [
"**Matplotlib**: Matplotlib is a comprehensive library for creating static, animated, and interactive visualizations in Python.\n",
"\n",
"Website: https://matplotlib.org/\n",
"\n",
"GitHub: https://github.com/matplotlib/matplotlib"
],
"metadata": {
"id": "lAD3BGBj7FhM"
}
},
{
"cell_type": "markdown",
"source": [
"In the previous notebook, we saw some basic examples of plotting and visualization in the context of learning `numpy`. In this notebook, we dive much deeper. The goal is to understand how `matplotlib` represents figures internally."
],
"metadata": {
"id": "nSAGGTtj7K7S"
}
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "xreaix9UpPl9"
},
"outputs": [],
"source": [
"from matplotlib import pyplot as plt\n",
"import numpy as np\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"source": [
"## Figure and Axes\n",
"The *figure* is the highest level of organization of `matplotlib` objects. If we want, we can create a figure explicitly."
],
"metadata": {
"id": "dc20UW3a7aTO"
}
},
{
"cell_type": "code",
"source": [
"fig = plt.figure()"
],
"metadata": {
"id": "OoCYwt3u7gpw"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig = plt.figure(figsize=(13, 5))"
],
"metadata": {
"id": "FBd1g03m7lk5"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig = plt.figure()\n",
"ax = fig.add_axes([0, 0, 1, 1])"
],
"metadata": {
"id": "FP4VUPXs8duF"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig = plt.figure()\n",
"ax = fig.add_axes([0, 0, 0.5, 1])"
],
"metadata": {
"id": "lYeUE0AC8h-Z"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig = plt.figure()\n",
"ax1 = fig.add_axes([0, 0, 0.5, 1])\n",
"ax2 = fig.add_axes([0.6, 0, 0.3, 0.5], facecolor='g')"
],
"metadata": {
"id": "okaAc9VQ8oPq"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Subplots\n",
"Subplot syntax is one way to specify the creation of multiple axes."
],
"metadata": {
"id": "L9X_blGe8wPS"
}
},
{
"cell_type": "code",
"source": [
"fig = plt.figure()\n",
"axes = fig.subplots(nrows=2, ncols=3)"
],
"metadata": {
"id": "NRoBDo7D8wyn"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig = plt.figure(figsize=(12, 6))\n",
"axes = fig.subplots(nrows=2, ncols=3)"
],
"metadata": {
"id": "jDAqJvF880ZA"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"axes"
],
"metadata": {
"id": "2bFPto179F3O"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"There is a shorthand for doing this all at once, **which is our recommended way to create new figures!**"
],
"metadata": {
"id": "y5EVIEJZ9JFr"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()"
],
"metadata": {
"id": "MMeoNrAj9Gun"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig, axes = plt.subplots(ncols=2, figsize=(8, 4), subplot_kw={'facecolor': 'g'})"
],
"metadata": {
"id": "_pt6ETt99NCL"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"axes"
],
"metadata": {
"id": "cQslx-qS9QA9"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Drawing into Axes\n",
"All plots are drawn into axes. It is easiest to understand how matplotlib works if you use the [object-oriented](https://matplotlib.org/faq/usage_faq.html#coding-styles) style."
],
"metadata": {
"id": "LgruyuT992q2"
}
},
{
"cell_type": "code",
"source": [
"# create some data to plot\n",
"import numpy as np\n",
"x = np.linspace(-np.pi, np.pi, 100)\n",
"y = np.cos(x)\n",
"z = np.sin(6*x)"
],
"metadata": {
"id": "lwztA08a93Ly"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(x, y)"
],
"metadata": {
"id": "nx6PAGZ69Qw3"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"This does the same thing as"
],
"metadata": {
"id": "mqrmDZ8J-TKV"
}
},
{
"cell_type": "code",
"source": [
"plt.plot(x, y)"
],
"metadata": {
"id": "AJUfAIRj-Ubq"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"This starts to matter when we have multiple axes to worry about."
],
"metadata": {
"id": "HVvEl6ry-ZWJ"
}
},
{
"cell_type": "code",
"source": [
"fig, axes = plt.subplots(figsize=(8, 4), ncols=2)\n",
"ax0, ax1 = axes\n",
"ax0.plot(x, y)\n",
"ax1.plot(x, z)"
],
"metadata": {
"id": "VVRTS3_Q-Z_A"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Labeling Plots"
],
"metadata": {
"id": "3z3x4GrJ_1dM"
}
},
{
"cell_type": "code",
"source": [
"fig, axes = plt.subplots(figsize=(8, 4), ncols=2)\n",
"ax0, ax1 = axes\n",
"\n",
"ax0.plot(x, y)\n",
"ax0.set_xlabel('x')\n",
"ax0.set_ylabel('y')\n",
"ax0.set_title('x vs. y')\n",
"\n",
"ax1.plot(x, z)\n",
"ax1.set_xlabel('x')\n",
"ax1.set_ylabel('z')\n",
"ax1.set_title('x vs. z')\n",
"\n",
"# squeeze everything in\n",
"plt.tight_layout()"
],
"metadata": {
"id": "AUkT0RkH_4_L"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Customizing Line Plots"
],
"metadata": {
"id": "xfCIldSk_837"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(x, y, x, z)"
],
"metadata": {
"id": "hcb8cjBC_9Yf"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# It’s simple to switch axes\n",
"fig, ax = plt.subplots()\n",
"ax.plot(y, x, z, x)"
],
"metadata": {
"id": "96wy0fxfAJd_"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"A “parametric” graph:"
],
"metadata": {
"id": "1bNkiVTnAS3d"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(y, z)"
],
"metadata": {
"id": "sE5uNA37AWRS"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Line Styles"
],
"metadata": {
"id": "wgcuuutlAYqB"
}
},
{
"cell_type": "code",
"source": [
"fig, axes = plt.subplots(figsize=(16, 5), ncols=3)\n",
"axes[0].plot(x, y, linestyle='dashed')\n",
"axes[0].plot(x, z, linestyle='--')\n",
"\n",
"axes[1].plot(x, y, linestyle='dotted')\n",
"axes[1].plot(x, z, linestyle=':')\n",
"\n",
"axes[2].plot(x, y, linestyle='dashdot', linewidth=5)\n",
"axes[2].plot(x, z, linestyle='-.', linewidth=0.5)"
],
"metadata": {
"id": "E2sopQz2AaXC"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Colors\n",
"As described in the [colors documentation](https://matplotlib.org/2.0.2/api/colors_api.html), there are some special codes for commonly used colors:\n",
"\n",
"* b: blue\n",
"* g: green\n",
"* r: red\n",
"* c: cyan\n",
"* m: magenta\n",
"* y: yellow\n",
"* k: black\n",
"* w: white"
],
"metadata": {
"id": "DjKpliLzAfyz"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(x, y, color='k')\n",
"ax.plot(x, z, color='r')"
],
"metadata": {
"id": "fDfQAo-6AiKA"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"# Other ways to specify colors:\n",
"\n",
"fig, axes = plt.subplots(figsize=(16, 5), ncols=3)\n",
"\n",
"# grayscale\n",
"axes[0].plot(x, y, color='0.8')\n",
"axes[0].plot(x, z, color='0.2')\n",
"\n",
"# RGB tuple\n",
"axes[1].plot(x, y, color=(1, 0, 0.7))\n",
"axes[1].plot(x, z, color=(0, 0.4, 0.3))\n",
"\n",
"# HTML hex code\n",
"axes[2].plot(x, y, color='#00dcba')\n",
"axes[2].plot(x, z, color='#b029ee')"
],
"metadata": {
"id": "oUf1EXXhAmf_"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"There is a default color cycle built into `matplotlib`."
],
"metadata": {
"id": "9-BGfYSRAlsI"
}
},
{
"cell_type": "code",
"source": [
"plt.rcParams['axes.prop_cycle']"
],
"metadata": {
"id": "MpcIyAo0Au9R"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12, 10))\n",
"for factor in np.linspace(0.2, 1, 11):\n",
" ax.plot(x, factor*y)"
],
"metadata": {
"id": "_KUyWyuoAzTH"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Markers\n",
"There are [lots of different markers](https://matplotlib.org/api/markers_api.html) availabile in matplotlib!"
],
"metadata": {
"id": "20PihEvkA03R"
}
},
{
"cell_type": "code",
"source": [
"fig, axes = plt.subplots(figsize=(12, 5), ncols=2)\n",
"\n",
"axes[0].plot(x[:20], y[:20], marker='.')\n",
"axes[0].plot(x[:20], z[:20], marker='o')\n",
"\n",
"axes[1].plot(x[:20], z[:20], marker='^',\n",
" markersize=10, markerfacecolor='r',\n",
" markeredgecolor='k')"
],
"metadata": {
"id": "EAdhXPJ9A4qO"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Label, Ticks, and Gridlines"
],
"metadata": {
"id": "MatcZn-SA6rv"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.plot(x, y)\n",
"\n",
"ax.set_xlabel('x')\n",
"ax.set_ylabel('y')\n",
"ax.set_title(r'A complicated math function: $f(x) = \\cos(x)$')\n",
"\n",
"ax.set_xticks(np.pi * np.array([-1, 0, 1]))\n",
"ax.set_xticklabels([r'$-\\pi$', '0', r'$\\pi$'])\n",
"ax.set_yticks([-1, 0, 1])\n",
"\n",
"ax.set_yticks(np.arange(-1, 1.1, 0.2), minor=True)\n",
"#ax.set_xticks(np.arange(-3, 3.1, 0.2), minor=True)\n",
"\n",
"ax.grid(which='minor', linestyle='--')\n",
"ax.grid(which='major', linewidth=2)"
],
"metadata": {
"id": "d2i2nad4A-1C"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Axis Limits"
],
"metadata": {
"id": "RQsvvs8NBChU"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(x, y, x, z)\n",
"ax.set_xlim(-5, 5)\n",
"ax.set_ylim(-3, 3)"
],
"metadata": {
"id": "uC3c4HoVBDGm"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Text Annotations"
],
"metadata": {
"id": "lIeB1w03BFsf"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(x, y)\n",
"ax.text(-3, 0.3, 'hello world')\n",
"ax.annotate('the maximum', xy=(0, 1),\n",
" xytext=(0, 0), arrowprops={'facecolor': 'k'})"
],
"metadata": {
"id": "ZQHNvu6bBGMv"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Scatter Plots"
],
"metadata": {
"id": "hAavOonBCk2r"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"splot = ax.scatter(y, z, c=x, s=(100*z**2 + 5))\n",
"fig.colorbar(splot)"
],
"metadata": {
"id": "mjx09x18CorB"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Bar Plots"
],
"metadata": {
"id": "2lNIhHysCxJF"
}
},
{
"cell_type": "code",
"source": [
"labels = ['first', 'second', 'third']\n",
"values = [10, 5, 30]\n",
"\n",
"fig, axes = plt.subplots(figsize=(10, 5), ncols=2)\n",
"axes[0].bar(labels, values)\n",
"axes[1].barh(labels, values)"
],
"metadata": {
"id": "vertirNpCx4Y"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"\n",
"## 2D Plotting Methods"
],
"metadata": {
"id": "1WfcckZRC3dL"
}
},
{
"cell_type": "markdown",
"source": [
"### `imshow`"
],
"metadata": {
"id": "xqpZz41ODAeu"
}
},
{
"cell_type": "code",
"source": [
"## imshow\n",
"x1d = np.linspace(-2*np.pi, 2*np.pi, 100)\n",
"y1d = np.linspace(-np.pi, np.pi, 50)\n",
"xx, yy = np.meshgrid(x1d, y1d)\n",
"f = np.cos(xx) * np.sin(yy)\n",
"print(f.shape)"
],
"metadata": {
"id": "nDF_XWSBCzzY"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12,4), ncols=2)\n",
"ax[0].imshow(f)\n",
"ax[1].imshow(f, origin='lower')"
],
"metadata": {
"id": "bMcztn89Ct_t"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### `pcolormesh`"
],
"metadata": {
"id": "TvPiUweqDDHt"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(ncols=2, figsize=(12, 5))\n",
"pc0 = ax[0].pcolormesh(x1d, y1d, f)\n",
"pc1 = ax[1].pcolormesh(xx, yy, f)\n",
"fig.colorbar(pc0, ax=ax[0])\n",
"fig.colorbar(pc1, ax=ax[1])"
],
"metadata": {
"id": "uw-NeHU5DFGj"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"x_sm, y_sm, f_sm = xx[:10, :10], yy[:10, :10], f[:10, :10]\n",
"\n",
"fig, ax = plt.subplots(figsize=(12,5), ncols=2)\n",
"\n",
"# last row and column ignored!\n",
"ax[0].pcolormesh(x_sm, y_sm, f_sm, edgecolors='k')\n",
"\n",
"# same!\n",
"ax[1].pcolormesh(x_sm, y_sm, f_sm[:-1, :-1], edgecolors='k')"
],
"metadata": {
"id": "XO--X1iXDIw7"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"y_distorted = y_sm*(1 + 0.1*np.cos(6*x_sm))\n",
"\n",
"plt.figure(figsize=(12,6))\n",
"plt.pcolormesh(x_sm, y_distorted, f_sm[:-1, :-1], edgecolors='w')\n",
"plt.scatter(x_sm, y_distorted, c='k')"
],
"metadata": {
"id": "XnUxk88BDLb-"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### `contour` / `contourf`"
],
"metadata": {
"id": "8Uqi_2msDhk0"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12, 5), ncols=2)\n",
"\n",
"# same thing!\n",
"ax[0].contour(x1d, y1d, f)\n",
"ax[1].contour(xx, yy, f)"
],
"metadata": {
"id": "2qPLGijEDi9G"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12, 5), ncols=2)\n",
"\n",
"c0 = ax[0].contour(xx, yy, f, 5)\n",
"c1 = ax[1].contour(xx, yy, f, 20)\n",
"\n",
"plt.clabel(c0, fmt='%2.1f')\n",
"plt.colorbar(c1, ax=ax[1])"
],
"metadata": {
"id": "W8xugK8XDjeI"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12, 5), ncols=2)\n",
"\n",
"clevels = np.arange(-1, 1, 0.2) + 0.1\n",
"\n",
"cf0 = ax[0].contourf(xx, yy, f, clevels, cmap='RdBu_r', extend='both')\n",
"cf1 = ax[1].contourf(xx, yy, f, clevels, cmap='inferno', extend='both')\n",
"\n",
"fig.colorbar(cf0, ax=ax[0])\n",
"fig.colorbar(cf1, ax=ax[1])"
],
"metadata": {
"id": "VqiskMXBDML3"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### `quiver`"
],
"metadata": {
"id": "MakDGbD2DnYR"
}
},
{
"cell_type": "code",
"source": [
"u = -np.cos(xx) * np.cos(yy)\n",
"v = -np.sin(xx) * np.sin(yy)\n",
"\n",
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.contour(xx, yy, f, clevels, cmap='RdBu_r', extend='both', zorder=0)\n",
"ax.quiver(xx[::4, ::4], yy[::4, ::4],\n",
" u[::4, ::4], v[::4, ::4], zorder=1)"
],
"metadata": {
"id": "TpJz3La8DqXs"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### streamplot"
],
"metadata": {
"id": "K-E4nNL8Dvil"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(figsize=(12, 7))\n",
"ax.streamplot(xx, yy, u, v, density=2, color=(u**2 + v**2))"
],
"metadata": {
"id": "PpPgXXLpDtAM"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Cartopy"
],
"metadata": {
"id": "7Or4dfNIWjxa"
}
},
{
"cell_type": "markdown",
"source": [
""
],
"metadata": {
"id": "BvTLSujHWemf"
}
},
{
"cell_type": "markdown",
"source": [
"Making maps is a fundamental part of geoscience research. Maps differ from regular figures in the following principle ways:\n",
"\n",
"* Maps require a *projection* of geographic coordinates on the 3D Earth to the 2D space of your figure.\n",
"\n",
"* Maps often include extra decorations besides just our data (e.g. continents, country borders, etc.)\n",
"\n",
"Mapping is a notoriously hard and complicated problem, mostly due to the complexities of projection.\n",
"\n",
"In this lecture, we will learn about [Cartopy](https://scitools.org.uk/cartopy/docs/latest/), one of the most common packages for making maps within python. Another popular and powerful library is [Basemap](https://matplotlib.org/basemap/); however, Basemap is [going away](https://matplotlib.org/basemap/users/intro.html#cartopy-new-management-and-eol-announcement) and being replaced with Cartopy in the near future. For this reason, new python learners are recommended to learn Cartopy."
],
"metadata": {
"id": "xu9BmV3EWtsu"
}
},
{
"cell_type": "markdown",
"source": [
"A lot of the material in this lesson was adopted from [Phil Elson](https://pelson.github.io/)’s excellent [Cartopy Tutorial](https://github.com/SciTools/cartopy-tutorial). Phil is the creator of Cartopy and published his tutorial under an [open license](http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/)."
],
"metadata": {
"id": "afQu9zlQWxc1"
}
},
{
"cell_type": "markdown",
"source": [
"## Background: Projections\n",
"\n",
"Most of our media for visualization *are* flat. Our two most common media are flat:\n",
"\n",
"* Paper\n",
"\n",
"* Screen"
],
"metadata": {
"id": "Jt651s2fW1aP"
}
},
{
"cell_type": "markdown",
"source": [
""
],
"metadata": {
"id": "Kgy-Wpo2W9qB"
}
},
{
"cell_type": "markdown",
"source": [
"## Introducing Cartopy\n",
"[Cartopy](https://scitools.org.uk/cartopy/docs/latest/) makes use of the powerful [PROJ](https://proj.org/), [numpy](https://numpy.org/doc/) and [shapely](https://shapely.readthedocs.io/en/stable/manual.html) libraries and includes a programatic interface built on top of [matplotlib](https://matplotlib.org/stable/) for the creation of publication quality maps.\n",
"\n",
"Key features of cartopy are its object oriented projection definitions, and its ability to transform points, lines, vectors, polygons and images between those projections."
],
"metadata": {
"id": "z1v0fGDmXH7s"
}
},
{
"cell_type": "markdown",
"source": [
"### Cartopy Projections and Other Reference Systems\n",
"In Cartopy, each projection is a class. Most classes of projection can be configured in projection-specific ways, although Cartopy takes an opinionated stance on sensible defaults.\n",
"\n",
"Let’s create a `Plate Carree` projection instance.\n",
"\n",
"To do so, we need cartopy’s `crs` module. This is typically imported as `ccrs` (Cartopy Coordinate Reference Systems).\n",
"\n",
"First, let's install `cartopy` to use it in this Google Colab notebook. 💻 We have to [first install `shapely` without binary](https://github.com/SciTools/cartopy/issues/1490), which hopefully should prevent your Google Colab notebook from crashing."
],
"metadata": {
"id": "uyrQ2GxvXNWb"
}
},
{
"cell_type": "code",
"source": [
"!pip install --no-binary 'shapely==1.6.4' 'shapely==1.6.4' --force\n",
"!pip install cartopy # Install latest version of cartopy"
],
"metadata": {
"id": "4VBdFV0xXFQz"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"import cartopy.crs as ccrs\n",
"import cartopy"
],
"metadata": {
"id": "mJFz7chgWue8"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"Cartopy’s projection list tells us that the `Plate Carree` [projection](https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html) is available with the `ccrs.PlateCarree` class:\n",
"\n",
"**Note**: we need to instantiate the class in order to do anything related to projections with it!"
],
"metadata": {
"id": "e_ESykiZXkyb"
}
},
{
"cell_type": "code",
"source": [
"ccrs.PlateCarree()"
],
"metadata": {
"id": "goIQ5Q5gXljq"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Drawing a map\n",
"Cartopy optionally depends upon matplotlib, and each projection knows how to create a matplotlib Axes (or AxesSubplot) that can represent itself.\n",
"\n",
"The Axes that the projection creates is a [cartopy.mpl.geoaxes.GeoAxes](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html). This Axes subclass overrides some of matplotlib’s existing methods, and adds a number of extremely useful ones for drawing maps.\n",
"\n",
"We’ll go back and look at those methods shortly, but first, let’s actually see the cartopy+matplotlib dance in action:"
],
"metadata": {
"id": "jLbKA4_fYO_E"
}
},
{
"cell_type": "code",
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"plt.axes(projection=ccrs.PlateCarree())"
],
"metadata": {
"id": "us1ATjnrYRLv"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"That was a little underwhelming, but we can see that the Axes created is indeed one of those GeoAxes[Subplot] instances.\n",
"\n",
"One of the most useful methods that this class adds on top of the standard matplotlib Axes class is the `coastlines` method. With no arguments, it will add the Natural Earth `1:110,000,000` scale coastline data to the map."
],
"metadata": {
"id": "E9JuDLL5YUQI"
}
},
{
"cell_type": "code",
"source": [
"plt.figure()\n",
"ax = plt.axes(projection=ccrs.PlateCarree())\n",
"ax.coastlines()"
],
"metadata": {
"id": "n-fkt30WYU6P"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"We could have also created a `matplotlib` subplot with one of the many approaches that exist. For example, the `plt.subplots` function could be used:"
],
"metadata": {
"id": "3jvbD-Q-YYIp"
}
},
{
"cell_type": "code",
"source": [
"fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})\n",
"ax.coastlines()"
],
"metadata": {
"id": "06GXXFY4YZUA"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"Projection classes have options we can use to customize the map"
],
"metadata": {
"id": "LWV3JxjQYbIw"
}
},
{
"cell_type": "code",
"source": [
"ccrs.PlateCarree?"
],
"metadata": {
"id": "xBSpLBo9YdQx"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))\n",
"ax.coastlines()"
],
"metadata": {
"id": "2_4CASQIYflw"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Useful Methods of a GeoAxes"
],
"metadata": {
"id": "ySthvmLuYiZJ"
}
},
{
"cell_type": "markdown",
"source": [
"The [cartopy.mpl.geoaxes.GeoAxes](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html) class adds a number of useful methods.\n",
"\n",
"Let’s take a look at:\n",
"\n",
"* [set_global](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.set_global) - zoom the map out as much as possible\n",
"\n",
"* [set_extent](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.set_extent) - zoom the map to the given bounding box\n",
"\n",
"* [gridlines](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.gridlines) - add gridlines (and optionally labels) to the axes\n",
"\n",
"* [coastlines](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.coastlines) - add Natural Earth coastlines to the axes\n",
"\n",
"* [stock_img](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.stock_img) - add a low-resolution Natural Earth background image to the axes\n",
"\n",
"* [imshow](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.imshow) - add an image (numpy array) to the axes\n",
"\n",
"* [add_geometries](https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html?highlight=set_global#cartopy.mpl.geoaxes.GeoAxes.add_geometries) - add a collection of geometries (Shapely) to the axes"
],
"metadata": {
"id": "I1N4_d2qYlaD"
}
},
{
"cell_type": "markdown",
"source": [
"### Some More Examples of Different Global Projections"
],
"metadata": {
"id": "YfVrmmbnYqBt"
}
},
{
"cell_type": "code",
"source": [
"projections = [ccrs.PlateCarree(),\n",
" ccrs.Robinson(),\n",
" ccrs.Mercator(),\n",
" ccrs.Orthographic(),\n",
" ccrs.InterruptedGoodeHomolosine()\n",
" ]\n",
"\n",
"\n",
"for proj in projections:\n",
" plt.figure()\n",
" ax = plt.axes(projection=proj)\n",
" ax.stock_img()\n",
" ax.coastlines()\n",
" ax.set_title(f'{type(proj)}')"
],
"metadata": {
"id": "lBgffMSjYsZD"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Regional Maps\n",
"To create a regional map, we use the `set_extent` method of GeoAxis to limit the size of the region."
],
"metadata": {
"id": "Lp5y7cNZYot5"
}
},
{
"cell_type": "code",
"source": [
"central_lon, central_lat = -10, 45\n",
"extent = [-40, 20, 30, 60]\n",
"ax = plt.axes(projection=ccrs.Orthographic(central_lon, central_lat))\n",
"ax.set_extent(extent)\n",
"ax.gridlines()\n",
"ax.coastlines(resolution='50m')"
],
"metadata": {
"id": "Lk5udJLiYwQ8"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Adding Features to the Map\n",
"To give our map more styles and details, we add `cartopy.feature` objects. Many useful features are built in. These “default features” are at coarse (110m) resolution.\n",
"\n",
"| **Name** | **Description** |\n",
"| ------------------------- | --------------- |\n",
"| `cartopy.feature.BORDERS` | Country boundaries |\n",
"| `cartopy.feature.COASTLINE` | Coastline, including major islands |\n",
"| `cartopy.feature.LAKES` | Natural and artificial lakes |\n",
"| `cartopy.feature.LAND` | Land polygons, including major islands |\n",
"| `cartopy.feature.OCEAN` | Ocean polygons |\n",
"| `cartopy.feature.RIVERS` | Single-line drainages, including lake centerlines |"
],
"metadata": {
"id": "n_9hozbrY2vV"
}
},
{
"cell_type": "markdown",
"source": [
"Below we illustrate these features in a customized map of North America."
],
"metadata": {
"id": "ZReBYSH9Y-Ux"
}
},
{
"cell_type": "code",
"source": [
"import cartopy.feature as cfeature\n",
"import numpy as np\n",
"\n",
"central_lat = 37.5\n",
"central_lon = -96\n",
"extent = [-120, -70, 24, 50.5]\n",
"central_lon = np.mean(extent[:2])\n",
"central_lat = np.mean(extent[2:])\n",
"\n",
"plt.figure(figsize=(12, 6))\n",
"ax = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))\n",
"ax.set_extent(extent)\n",
"\n",
"ax.add_feature(cartopy.feature.OCEAN)\n",
"ax.add_feature(cartopy.feature.LAND, edgecolor='black')\n",
"ax.add_feature(cartopy.feature.LAKES, edgecolor='black')\n",
"ax.add_feature(cartopy.feature.RIVERS)\n",
"ax.gridlines()"
],
"metadata": {
"id": "m_GVXZUkY8Fn"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"If we want higher-resolution features, Cartopy can automatically download and create them from the [Natural Earth Data database](http://www.naturalearthdata.com/) or the [GSHHS dataset](https://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html) database."
],
"metadata": {
"id": "OMdy7CNGZBdg"
}
},
{
"cell_type": "code",
"source": [
"rivers_50m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '50m')\n",
"\n",
"plt.figure(figsize=(12, 6))\n",
"ax = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))\n",
"ax.set_extent(extent)\n",
"\n",
"ax.add_feature(cartopy.feature.OCEAN)\n",
"ax.add_feature(cartopy.feature.LAND, edgecolor='black')\n",
"ax.add_feature(rivers_50m, facecolor='None', edgecolor='b')\n",
"ax.gridlines()"
],
"metadata": {
"id": "4u18kP_HZCCU"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Adding Data to the Map\n",
"\n",
"Now that we know how to create a map, let’s add our data to it! That’s the whole point.\n",
"\n",
"Because our map is a matplotlib axis, we can use all the familiar maptplotlib commands to make plots. By default, the map extent will be adjusted to match the data. We can override this with the `.set_global` or `.set_extent` commands."
],
"metadata": {
"id": "9NU4FlBnZFk9"
}
},
{
"cell_type": "code",
"source": [
"# create some test data\n",
"new_york = dict(lon=-74.0060, lat=40.7128)\n",
"honolulu = dict(lon=-157.8583, lat=21.3069)\n",
"lons = [new_york['lon'], honolulu['lon']]\n",
"lats = [new_york['lat'], honolulu['lat']]"
],
"metadata": {
"id": "YwCGYVsDZKQ0"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"**Key point**: the data also have to be transformed to the projection space.\n",
"\n",
"This is done via the `transform=` keyword in the plotting method. The argument is another `cartopy.crs` object. If you don’t specify a transform, Cartopy assume that the data is using the same projection as the underlying GeoAxis.\n",
"\n",
"From the [Cartopy Documentation](https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html):\n",
"\n",
"\"The core concept is that the projection of your axes is independent of the coordinate system your data is defined in. The `projection` argument is used when creating plots and determines the projection of the resulting plot (i.e. what the plot looks like). The `transform` argument to plotting functions tells Cartopy what coordinate system your data are defined in.\""
],
"metadata": {
"id": "mKY1JpaMZMYw"
}
},
{
"cell_type": "code",
"source": [
"ax = plt.axes(projection=ccrs.PlateCarree())\n",
"ax.plot(lons, lats, label='Equirectangular straight line')\n",
"ax.plot(lons, lats, label='Great Circle', transform=ccrs.Geodetic())\n",
"ax.coastlines()\n",
"ax.legend()\n",
"ax.set_global()"
],
"metadata": {
"id": "ZizwsnAxZP3S"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plotting 2D (Raster) Data\n",
"The same principles apply to 2D data. Below we create some example data defined in regular lat / lon coordinates."
],
"metadata": {
"id": "TujjC3c5ZTCL"
}
},
{
"cell_type": "code",
"source": [
"import numpy as np\n",
"lon = np.linspace(-80, 80, 25)\n",
"lat = np.linspace(30, 70, 25)\n",
"lon2d, lat2d = np.meshgrid(lon, lat)\n",
"data = np.cos(np.deg2rad(lat2d) * 4) + np.sin(np.deg2rad(lon2d) * 4)\n",
"plt.contourf(lon2d, lat2d, data)"
],
"metadata": {
"id": "6ZVIdqG6ZYXc"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"Now we create a `PlateCarree` projection and plot the data on it without any `transform` keyword. This happens to work because `PlateCarree` is the simplest projection of lat / lon data."
],
"metadata": {
"id": "mIXnMCUOZXE1"
}
},
{
"cell_type": "code",
"source": [
"ax = plt.axes(projection=ccrs.PlateCarree())\n",
"ax.set_global()\n",
"ax.coastlines()\n",
"ax.contourf(lon, lat, data)"
],
"metadata": {
"id": "bB5NYSg5ZbTM"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"However, if we try the same thing with a different projection, we get the wrong result."
],
"metadata": {
"id": "XZbEEdt0ZfGj"
}
},
{
"cell_type": "code",
"source": [
"projection = ccrs.RotatedPole(pole_longitude=-177.5, pole_latitude=37.5)\n",
"ax = plt.axes(projection=projection)\n",
"ax.set_global()\n",
"ax.coastlines()\n",
"ax.contourf(lon, lat, data)"
],
"metadata": {
"id": "Oqr9tFLuZg3v"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"To fix this, we need to pass the correct transform argument to `contourf`:"
],
"metadata": {
"id": "yxEdZVX7Zi5n"
}
},
{
"cell_type": "code",
"source": [
"projection = ccrs.RotatedPole(pole_longitude=-177.5, pole_latitude=37.5)\n",
"ax = plt.axes(projection=projection)\n",
"ax.set_global()\n",
"ax.coastlines()\n",
"ax.contourf(lon, lat, data, transform=ccrs.PlateCarree())"
],
"metadata": {
"id": "uuDkL-fVZjeC"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Showing Images"
],
"metadata": {
"id": "AOZTH4MwZnbw"
}
},
{
"cell_type": "markdown",
"source": [
"We can plot a satellite image easily on a map if we know its extent\n"
],
"metadata": {
"id": "T1YKQz5-ZxMS"
}
},
{
"cell_type": "code",
"source": [
"! wget https://lance-modis.eosdis.nasa.gov/imagery/gallery/2012270-0926/Miriam.A2012270.2050.2km.jpg"
],
"metadata": {
"id": "K_NPeEy1ZwGb"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"fig = plt.figure(figsize=(8, 12))\n",
"\n",
"# this is from the cartopy docs\n",
"fname = 'Miriam.A2012270.2050.2km.jpg'\n",
"img_extent = (-120.67660000000001, -106.32104523100001, 13.2301484511245, 30.766899999999502)\n",
"img = plt.imread(fname)\n",
"\n",
"ax = plt.axes(projection=ccrs.PlateCarree())\n",
"\n",
"# set a margin around the data\n",
"ax.set_xmargin(0.05)\n",
"ax.set_ymargin(0.10)\n",
"\n",
"# add the image. Because this image was a tif, the \"origin\" of the image is in the\n",
"# upper left corner\n",
"ax.imshow(img, origin='upper', extent=img_extent, transform=ccrs.PlateCarree())\n",
"ax.coastlines(resolution='50m', color='black', linewidth=1)\n",
"\n",
"# mark a known place to help us geo-locate ourselves\n",
"ax.plot(-117.1625, 32.715, 'bo', markersize=7, transform=ccrs.Geodetic())\n",
"ax.text(-117, 33, 'San Diego', transform=ccrs.Geodetic())"
],
"metadata": {
"id": "JTtPyLyUZn26"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Bonus: Xarray Integration\n",
"Cartopy transforms can be passed to `xarray`! This creates a very quick path for creating professional looking maps from netCDF data."
],
"metadata": {
"id": "hLExt3imZ5np"
}
},
{
"cell_type": "code",
"source": [
"import pooch\n",
"import urllib.request\n",
"datafile = pooch.retrieve('https://unils-my.sharepoint.com/:u:/g/personal/tom_beucler_unil_ch/EdY_7ttKxltLp-u6r1MpDggBvqGb355_VecnfPRLEIsMog?download=1',\n",
" known_hash='af8ff05bfeb8da2ec763773bfcc06112f20cb7167613359c98a2b80313c45b73')"
],
"metadata": {
"id": "muhhGm7mcDO4"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"import xarray as xr\n",
"ds = xr.open_dataset(datafile, drop_variables=['time_bnds'])\n",
"ds"
],
"metadata": {
"id": "JM7SeDgPZ-2T"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sst = ds.sst.sel(time='2000-01-01', method='nearest')\n",
"fig = plt.figure(figsize=(9,6))\n",
"ax = plt.axes(projection=ccrs.Robinson())\n",
"ax.coastlines()\n",
"ax.gridlines()\n",
"sst.plot(ax=ax, transform=ccrs.PlateCarree(),\n",
" vmin=2, vmax=30, cbar_kwargs={'shrink': 0.4})"
],
"metadata": {
"id": "G03F6GOAccMQ"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## Doing More\n",
"Browse the [Cartopy Gallery](https://scitools.org.uk/cartopy/docs/latest/gallery/index.html) to learn about all the different types of data and plotting methods available!"
],
"metadata": {
"id": "YAstBVw7cetj"
}
}
]
}