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

Introduction to Cartopy

\n", "

Unidata AMS 2021 Student Conference

\n", "\n", "
\n", "
\n", "\n", "---\n", "\n", "
\"Cartopy
\n", "\n", "\n", "### Focuses\n", "* Create a basic figure using CartoPy\n", "* Add maps to the figure\n", "* Plot georeferenced data on the figure\n", "\n", "\n", "### Objectives\n", "1. [Create a Map](#1.-Make-a-simple-map)\n", "1. [Add Features to Map](#2.-Add-features-to-map)\n", "1. [Plotting Data on Map](#3.-Plotting-data)\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "from metpy.calc import wind_speed\n", "from metpy.units import units\n", "from metpy.plots import USCOUNTIES" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Make a simple map\n", "\n", "Let's start with a simple map - without adding any additional parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Works with matplotlib's built-in transform support\n", "fig = plt.figure(figsize=(10,4))\n", "ax = fig.add_subplot(111, projection=ccrs.Robinson())\n", "\n", "# Set extent to cover the entire globe\n", "ax.set_global()\n", "\n", "# Add a stock image to the map - the background\n", "ax.stock_img()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have a basic map, let's add some parameters and zoom in a bit" ] }, { "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", "\n", "---\n", "\n", "## 2. Add features to map\n", "\n", "CartoPy has a few helper methods for adding maps to the plot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Setup the figure and geoaxes\n", "fig = plt.figure(figsize=(10,8))\n", "ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())\n", "\n", "# Add the stock image\n", "ax.stock_img()\n", "\n", "# Add coastline contours to the map\n", "ax.add_feature(cfeature.COASTLINE)\n", "\n", "# Set the extent\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In addition to coastlines, there are a variety of natural earth features one can add to a map" ] }, { "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, 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", "# Set the extent\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": {}, "source": [ "Interested in other map features? Check out the documentation from the [Natural Earth Project](http://www.naturalearthdata.com/)\n", "\n", "### US County Boundaries\n", "\n", "MetPy also has US County boundaries available at 20m, 5m, and 500k resolutions. Checkout the example below to see the difference between different resolutions." ] }, { "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')\n", " axis.set_title(scale)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "## 3. Plotting data\n", "\n", "Let's try plotting a point on the map near Boulder, Colorado using `ax.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.add_feature(cfeature.COASTLINE)\n", "ax.add_feature(cfeature.BORDERS, linewidth=2)\n", "ax.add_feature(cfeature.STATES, linestyle='--', edgecolor='black')\n", "\n", "# Add a point to the map using longitude, lat, and use a circle as the marker\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", "# Set the projection of the data point such that it transforms the point from lon, lat to the projected coordinate system Lambert Conformal\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 (ex. `contour`):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create some synthetic gridded wind data\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)) * units.knots\n", "u = np.zeros_like(v) * units.knots\n", "speed = wind_speed(u, v)\n", "\n", "# Create arrays of longitude and latitude\n", "x = np.linspace(-120, -60, 5)\n", "y = np.linspace(30, 55, 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)\n", "\n", "# Plot wind barbs--CartoPy handles reprojecting the vectors properly for the\n", "# coordinate system\n", "ax.barbs(x, y, u.m, v.m, transform=ccrs.PlateCarree(), color='tab:blue')\n", "\n", "ax.set_extent([-130, -60, 20, 55])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Check out these resources" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Interested in learning more about CartoPy? Be sure to check out the [CartoPy Example Gallery](https://scitools.org.uk/cartopy/docs/latest/gallery/index.html)\n", "* Also be sure to checkout the [general guide to CartoPy](https://scitools.org.uk/cartopy/docs/latest/#getting-started) which includes instructions on downloading to your local machine \n", "\n", "### Related Notebooks\n", "* [Matplotlib: Basics notebook](https://nbviewer.jupyter.org/github/Unidata/pyaos-ams-2021/blob/master/notebooks/visualization/matplotlib-basics.ipynb)" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:pyaos-ams-2021]", "language": "python", "name": "conda-env-pyaos-ams-2021-py" }, "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.9.1" } }, "nbformat": 4, "nbformat_minor": 4 }