{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "
\n", "\"Unidata\n", "
\n", "\n", "

Plotting on a Map with CartoPy

\n", "

Unidata Python Workshop

\n", "\n", "
\n", "
\n", "\n", "
\n", "\n", "
\"CartoPy\"
\n", "\n", "## Overview:\n", "\n", "* **Teaching:** 15 minutes\n", "* **Exercises:** 15 minutes\n", "\n", "### Questions\n", "1. How do we plot on a map in Python?\n", "1. How do I specify a map projection?\n", "1. How do I tell CartoPy how to reference my data?\n", "1. How do I add map features to a CartoPy plot?\n", "\n", "### Objectives\n", "1. Create a basic figure using CartoPy\n", "1. Add maps to the figure\n", "1. Plot georeferenced data on the figure" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "## 1. Basic CartoPy Plotting" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- High level API for dealing with maps\n", "- CartoPy allows you to plot data on a 2D map.\n", "- Support many different map projections\n", "- Support for shapefiles from the GIS world" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Set things up\n", "%matplotlib inline\n", "\n", "# Importing CartoPy\n", "import cartopy.crs as ccrs\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simplest plot we can make sets a projection with no parameters. The one below uses the Robinson projection:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Works with matplotlib's built-in transform support.\n", "fig = plt.figure(figsize=(10, 4))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())\n", "\n", "# Sets the extent to cover the whole globe\n", "ax.set_global()\n", "\n", "# Adds standard background map\n", "ax.stock_img()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We also have fine-tuned control over the globe used in the projection as well as lots of standard parameters, which depend on individual projections:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Set up a globe with a specific radius\n", "globe = ccrs.Globe(semimajor_axis=6371000.)\n", "\n", "# Set up a Lambert Conformal projection\n", "proj = ccrs.LambertConformal(standard_parallels=[25.0], globe=globe)\n", "\n", "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=proj)\n", "\n", "# Sets the extent using a lon/lat box\n", "ax.set_extent([-130, -60, 20, 55])\n", "\n", "ax.stock_img()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## 2. Adding maps to CartoPy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "CartoPy provides a couple helper methods for adding maps to the plot:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "ax.stock_img()\n", "\n", "# Specifies the detail level of the map.\n", "# Options are '110m' (default), '50m', and '10m'\n", "ax.coastlines(resolution='50m')\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Cartopy also has a lot of built-in support for a variety of map features:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# import cartopy's collection of map features\n", "import cartopy.feature as cfeature\n", "\n", "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "# Add variety of features\n", "ax.add_feature(cfeature.LAND)\n", "ax.add_feature(cfeature.OCEAN)\n", "ax.add_feature(cfeature.COASTLINE)\n", "\n", "# Can also supply matplotlib kwargs\n", "ax.add_feature(cfeature.BORDERS, linestyle=':')\n", "ax.add_feature(cfeature.STATES, linestyle=':')\n", "ax.add_feature(cfeature.LAKES, alpha=0.5)\n", "ax.add_feature(cfeature.RIVERS, edgecolor='tab:green')\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The map features are available at several different scales depending on how large the area you are covering is. The scales can be accessed using the `with_scale` method. Natural Earth features are available at 110m, 50m and 10m." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "# Add variety of features\n", "ax.add_feature(cfeature.LAND)\n", "ax.add_feature(cfeature.OCEAN)\n", "ax.add_feature(cfeature.COASTLINE)\n", "\n", "# Can also supply matplotlib kwargs\n", "ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle=':')\n", "ax.add_feature(cfeature.STATES.with_scale('50m'), linestyle=':')\n", "ax.add_feature(cfeature.LAKES.with_scale('50m'), alpha=0.5)\n", "ax.add_feature(cfeature.RIVERS.with_scale('50m'), edgecolor='tab:green')\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "You can also grab other features from the Natural Earth project: http://www.naturalearthdata.com/\n", "\n", "## US Counties\n", "MetPy has US Counties built in at the 20m, 5m, and 500k resolutions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from metpy.plots import USCOUNTIES" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "proj = ccrs.LambertConformal(central_longitude=-85.0, central_latitude=45.0)\n", "\n", "fig = plt.figure(figsize=(12, 9))\n", "ax1 = fig.add_subplot(1, 3, 1, projection=proj)\n", "ax2 = fig.add_subplot(1, 3, 2, projection=proj)\n", "ax3 = fig.add_subplot(1, 3, 3, projection=proj)\n", "\n", "for scale, axis in zip(['20m', '5m', '500k'], [ax1, ax2, ax3]):\n", " axis.set_extent([270.25, 270.9, 38.15, 38.75], ccrs.Geodetic())\n", " axis.add_feature(USCOUNTIES.with_scale(scale), edgecolor='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## 3. Plotting Data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "CartoPy supports all of the matplotlib plotting options you would expect on a map. It handles transforming your data between different coordinate systems transparently, provided you provide the correct information. (More on this later...). To start, let's put a marker at -105, 40:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "ax.add_feature(cfeature.COASTLINE)\n", "ax.add_feature(cfeature.BORDERS, linewidth=2)\n", "ax.add_feature(cfeature.STATES, linestyle='--', edgecolor='black')\n", "\n", "ax.plot(-105, 40, marker='o', color='tab:red')\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So that did not succeed at putting a marker at -105 longitude, 40 latitude (Boulder, CO). Instead, what actually happened is that it put the marker at (-105, 40) in the map projection coordinate system; in this case that's a Lambert Conformal projection, and x,y are assumed in meters relative to the origin of that coordinate system. To get CartoPy to treat it as longitude/latitude, we need to tell it that's what we're doing. We do this through the use of the `transform` argument to all of the plotting functions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "ax.add_feature(cfeature.COASTLINE)\n", "ax.add_feature(cfeature.BORDERS, linewidth=2)\n", "ax.add_feature(cfeature.STATES, linestyle='--', edgecolor='black')\n", "\n", "data_projection = ccrs.PlateCarree()\n", "ax.plot(-105, 40, marker='o', color='tab:red', transform=data_projection)\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This approach by CartoPy separates the data coordinate system from the coordinate system of the plot. It allows you to take data in any coordinate system (lon/lat, Lambert Conformal) and display it in any map you want. It also allows you to combine data from various coordinate systems seamlessly. This extends to all plot types, not just `plot`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create some synthetic gridded wind data\n", "import numpy as np\n", "from metpy.calc import wind_speed\n", "\n", "# Note that all of these winds have u = 0 -> south wind\n", "v = np.full((5, 5), 10, dtype=np.float64) + 10 * np.arange(5)\n", "u = np.zeros_like(v)\n", "speed = wind_speed(u, v)\n", "x, y = np.meshgrid(np.linspace(-120, -60, 5), np.linspace(25, 50, 5))\n", "\n", "# Plot as normal\n", "fig = plt.figure(figsize=(10, 8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "ax.add_feature(cfeature.COASTLINE)\n", "ax.add_feature(cfeature.BORDERS, linewidth=2)\n", "\n", "# Plot wind barbs--CartoPy handles reprojecting the vectors properly for the\n", "# coordinate system\n", "ax.barbs(x, y, u, v, transform=ccrs.PlateCarree(), color='tab:blue')\n", "\n", "# Contour the wind speed\n", "contours = np.arange(5, 50, 10)\n", "ax.contour(x, y, speed, contours, colors='black',\n", " linestyles='--', transform=ccrs.PlateCarree())\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "* Create a map, on a Mercator Projection, which at least has coastlines and country and state borders. Bonus points for putting on colored land and oceans, or other map features.\n", "* Plot our location correctly on the map.\n", "* Set the bounds of the map to zoom in mostly over our state/region." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# YOUR CODE GOES HERE\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# %load solutions/map.py\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }