\n",
"\n",
"\n",
"\n",
"\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
}