diff --git a/.gitignore b/.gitignore index 82509949e..84ae88357 100644 --- a/.gitignore +++ b/.gitignore @@ -135,3 +135,4 @@ content/pages/communications.md content/pages/contributing.md content/pages/code_of_conduct.md content/notebooks_gallery/ +.jupyter_cache/ diff --git a/.jupyter_cache/global.db b/.jupyter_cache/global.db new file mode 100644 index 000000000..27e8c3dbd Binary files /dev/null and b/.jupyter_cache/global.db differ diff --git a/_config.yml b/_config.yml index a9af7c3f4..88acb41f9 100644 --- a/_config.yml +++ b/_config.yml @@ -5,10 +5,15 @@ title: Project Pythia Foundations author: The Project Pythia community logo: images/ProjectPythia_Logo_Final-01-Blue.png -# Force re-execution of notebooks on each build. -# See https://jupyterbook.org/content/execute.html execute: - execute_notebooks: force + # See https://jupyterbook.org/content/execute.html + execute_notebooks: cache + # "auto" should only execute those notebooks that don't have output in all cells. + # "force" : force execution of all notebooks + # "cache": Cache output of notebooks on each build. + timeout: 600 + # Do not abort notebook execution for errors, since they may be intentional. + allow_errors: True # Define the name of the latex output file for PDF builds latex: diff --git a/_toc.yml b/_toc.yml index 8385aec8e..5fcaab0c6 100644 --- a/_toc.yml +++ b/_toc.yml @@ -36,6 +36,8 @@ - file: core/numpy - file: core/matplotlib - file: core/cartopy + sections: + - file: core/cartopy/01_Cartopy_Intro - file: core/datetime - file: core/pandas - file: core/data-formats diff --git a/core/cartopy.md b/core/cartopy.md index efd3afeae..59889c1b6 100644 --- a/core/cartopy.md +++ b/core/cartopy.md @@ -4,4 +4,19 @@ This content is under construction! ``` -This section will contain tutorials on plotting map with [cartopy](https://scitools.org.uk/cartopy/docs/latest/). It will be cross-referenced with tutorials on [xarray](xarray) and [matplotlib](matplotlib). +This section contains tutorials on plotting maps with [cartopy](https://scitools.org.uk/cartopy/docs/latest/). +It will be cross-referenced with tutorials on [xarray](xarray) and [matplotlib](matplotlib). + +--- + +From the [Cartopy website](https://scitools.org.uk/cartopy/docs/latest): +Cartopy is a Python package designed for geospatial data processing in order to +produce maps and other geospatial data analyses. + +Cartopy makes use of the powerful PROJ.4, NumPy and Shapely libraries and includes a programmatic interface +built on top of Matplotlib for the creation of publication quality maps. + +You should have a basic familiarity with Matplotlib prior to working through the Cartopy notebooks presented here. + +Key features of cartopy are its object-oriented [projection definitions](https://scitools.org.uk/cartopy/docs/latest/crs/projections.html#cartopy-projections), +and its ability to transform points, lines, vectors, polygons and images between those projections. diff --git a/core/cartopy/01_Cartopy_Intro.ipynb b/core/cartopy/01_Cartopy_Intro.ipynb new file mode 100644 index 000000000..39c91ff08 --- /dev/null +++ b/core/cartopy/01_Cartopy_Intro.ipynb @@ -0,0 +1,578 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 01_Cartopy_Introduction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## In this notebook, we'll cover the following:\n", + "1. Extend Matplotlib's `axes` into georeferenced `GeoAxes`\n", + "2. Create a map with a specified projection and add various cartographic features to it.\n", + "3. Explore some of Cartopy's map projections\n", + "4. Create regional maps\n", + "\n", + "### References: \n", + "1. [Cartopy Documentation](https://scitools.org.uk/cartopy/docs/latest/)\n", + "2. [*Maps with Cartopy*(Ryan Abernathey)](https://rabernat.github.io/research_computing_2018/maps-with-cartopy.html)\n", + "3. [Map Projections (GeoCAT)](https://geocat-examples.readthedocs.io/en/latest/gallery/index.html#map-projections)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 0) Preliminaries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from cartopy import crs as ccrs, feature as cfeature" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1) Extend Matplotlib's `axes` into georeferenced `GeoAxes`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Recall that in Matplotlib, what we might tradtionally term a *figure* consists of two key components: a `figure` and an associated subplot `axes` instance.\n", + "\n", + "#### By virtue of importing Cartopy, we can now convert the `axes` into a `GeoAxes` by specifying a projection that we have imported from *Cartopy's Coordinate Reference System* class as `ccrs`. This will effectively *georeference* the subplot." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2) Create a map with a specified projection and add various cartographic features to it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Full list of projections in Cartopy: https://scitools.org.uk/cartopy/docs/latest/crs/projections.html " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(11, 8.5))\n", + "# Create a GeoAxes that uses the PlateCarree projection\n", + "# (basically a global lat-lon map projection, which translates from French to \"flat square\" in English,\n", + "# where each point is equally spaced in terms of degrees)\n", + "ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=-75))\n", + "ax.set_title(\"A Geo-referenced subplot, Plate Carree projection\")\n", + "# adding a trailing semicolon to the last line prevents output from the deployed axes method from appearing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Although the figure seems empty, it has in fact been georeferenced, using one of Cartopy's map projections that is provided by Cartopy's `crs` (coordinate reference system) class. We can now add in cartographic features, in the form of *shapefiles*, to our subplot. One of them is `coastlines`, which is a callable `GeoAxes` method that can be plotted directly on our subplot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ax.coastlines()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
When using the `%matplotlib inline` notebook magic, to get the figure to display again with the features that we've added since the original display, just type the name of the Figure object in its own cell.
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Cartopy provides other cartographic features via its `features` class, which we've imported as `cfeature`. These are also shapefiles, downloaded on initial request from http://www.naturalearthdata.com/ . Once downloaded, they \"live\" in your `~/.local/share/cartopy` directory." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### We add them to our subplot via the `add_feature` method. We can define attributes for them in a manner similar to Matplotlib's `plot` method. A list of the various Natural Earth shapefiles appears in https://scitools.org.uk/cartopy/docs/latest/matplotlib/feature_interface.html ." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='black')\n", + "ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once again, referencing the `Figure` object will re-render the figure in the notebook, now including the two features" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3) Explore some of Cartopy's map projections" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
You can find a list of supported projections in Cartopy, with examples, at https://scitools.org.uk/cartopy/docs/latest/crs/projections.html
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Mollweide Projection (often used with global satellite mosaics)\n", + "#### This time, we'll define an object to store our projection definition. Any time we wish to use this particular projection later in the notebook, we can use the object name rather than repeating the same call to `ccrs`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(11, 8.5))\n", + "projMoll = ccrs.Mollweide(central_longitude=0)\n", + "ax = plt.subplot(1, 1, 1, projection=projMoll)\n", + "ax.set_title(\"Mollweide Projection\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Add in the cartographic shapefiles\n", + "ax.coastlines()\n", + "ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue')\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Add a fancy background image to the map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ax.stock_img()\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Lambert Azimuthal Equal Area Projection" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(11, 8.5))\n", + "projLae = ccrs.LambertAzimuthalEqualArea(central_longitude=0.0, central_latitude=0.0)\n", + "ax = plt.subplot(1, 1, 1, projection=projLae)\n", + "ax.set_title(\"Lambert Azimuthal Equal Area Projection\")\n", + "ax.coastlines()\n", + "ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4) Create regional maps" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Now, let's go back to PlateCarree, but let's use Cartopy's `set_extent` method to restrict the map coverage to a North American view. Let's also choose a lower resolution for coastlines, just to illustrate how one can specify that. Plot lat/lon lines as well." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Reference for Natural Earth's three resolutions (10m, 50m, 110m; higher is coarser): http://www.naturalearthdata.com/downloads/ " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "projPC = ccrs.PlateCarree()\n", + "lonW = -140\n", + "lonE = -40\n", + "latS = 15\n", + "latN = 65\n", + "cLat = (latN + latS) / 2\n", + "cLon = (lonW + lonE) / 2\n", + "res = '110m'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(11, 8.5))\n", + "ax = plt.subplot(1, 1, 1, projection=projPC)\n", + "ax.set_title('Plate Carree')\n", + "gl = ax.gridlines(\n", + " draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'\n", + ")\n", + "ax.set_extent([lonW, lonE, latS, latN], crs=projPC)\n", + "ax.coastlines(resolution=res, color='black')\n", + "ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')\n", + "ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Did you notice the output above the map about *cartopy.mpl.feature_artist*? We can eliminate that in one of two ways:\n", + "1. Use a generic object name, `_`, to assign that last feature to.\n", + "2. Add a semicolon at the end of the last line in the cell." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
Note the in the `set_extent` call, we specified PlateCarree. This ensures that the values we passed into `set_extent` will be transformed from degrees into the values appropriate for the projection we use for the map.
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### The PlateCarree projection exaggerates the spatial extent of regions closer to the poles. Let's try a couple different projections. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "projStr = ccrs.Stereographic(central_longitude=cLon, central_latitude=cLat)\n", + "fig = plt.figure(figsize=(11, 8.5))\n", + "ax = plt.subplot(1, 1, 1, projection=projStr)\n", + "ax.set_title('Stereographic')\n", + "gl = ax.gridlines(\n", + " draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'\n", + ")\n", + "ax.set_extent([lonW, lonE, latS, latN], crs=projPC)\n", + "ax.coastlines(resolution=res, color='black')\n", + "ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')\n", + "# Use generic object name to suppress text output to the screen\n", + "_ = ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "projLcc = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)\n", + "fig = plt.figure(figsize=(11, 8.5))\n", + "ax = plt.subplot(1, 1, 1, projection=projLcc)\n", + "ax.set_title('Lambert Conformal')\n", + "gl = ax.gridlines(\n", + " draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'\n", + ")\n", + "ax.set_extent([lonW, lonE, latS, latN], crs=projPC)\n", + "ax.coastlines(resolution='110m', color='black')\n", + "ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')\n", + "# End last line with a semicolon to suppress text output to the screen\n", + "ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
Lat/lon labeling for projections other than Mercator and PlateCarree is a recent addition to Cartopy. As you can see, work still needs to be done to improve the placement of labels.
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create a regional map, centered over New York State and add in some more Natural Earth geographic features.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Set the domain for defining the plot region. We will use this in the `set_extent` line below. Since these coordinates are expressed in degrees, they correspond to the PlateCarree projection." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
Be patient: with a limited regional extent as specified here, the highest resolution (10m) shapefiles are used; as a result (as with any `GeoAxes` object that must be transformed from one coordinate system to another), this will take some time to plot (could be several minutes), and even longer if you haven't previously retrieved these features from the Natural Earth shapefile server.
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "latN = 45.2\n", + "latS = 40.2\n", + "lonW = -80.0\n", + "lonE = -71.5\n", + "cLat = (latN + latS) / 2\n", + "cLon = (lonW + lonE) / 2\n", + "projLccNY = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Some pre-defined Features exist as `cartopy.feature` constants. The resolution of these pre-defined Features will depend on the areal extent of your map, which you specify via `set_extent`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15, 10))\n", + "ax = plt.subplot(1, 1, 1, projection=projLccNY)\n", + "ax.set_extent([lonW, lonE, latS, latN], crs=projPC)\n", + "ax.add_feature(cfeature.LAND)\n", + "ax.add_feature(cfeature.OCEAN)\n", + "ax.add_feature(cfeature.COASTLINE)\n", + "ax.add_feature(cfeature.BORDERS, linestyle='--')\n", + "ax.add_feature(cfeature.LAKES, alpha=0.5)\n", + "ax.add_feature(cfeature.STATES)\n", + "ax.add_feature(cfeature.RIVERS)\n", + "ax.set_title('New York and Vicinity');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Let's create a new map, but this time use lower-resolution shapefiles from Natural Earth, and also eliminate plotting the country borders.\n", + "### Notice this is a bit more involved. First we create objects for our lower-resolution shapefiles via the `NaturalEarthFeature` method from Cartopy's `feature` class, and then we add them to the map with `add_feature`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15, 10))\n", + "ax = plt.subplot(1, 1, 1, projection=projLccNY)\n", + "ax.set_extent((lonW, lonE, latS, latN), crs=projPC)\n", + "\n", + "# The features with names such as cfeature.LAND, cfeature.OCEAN, are higher-resolution (10m)\n", + "# shapefiles from the Naturalearth repository. Lower resolution shapefiles (50m, 110m) can be\n", + "# used by using the cfeature.NaturalEarthFeature method as illustrated below.\n", + "\n", + "resolution = '110m'\n", + "\n", + "land_mask = cfeature.NaturalEarthFeature(\n", + " 'physical',\n", + " 'land',\n", + " scale=resolution,\n", + " edgecolor='face',\n", + " facecolor=cfeature.COLORS['land'],\n", + ")\n", + "sea_mask = cfeature.NaturalEarthFeature(\n", + " 'physical',\n", + " 'ocean',\n", + " scale=resolution,\n", + " edgecolor='face',\n", + " facecolor=cfeature.COLORS['water'],\n", + ")\n", + "lake_mask = cfeature.NaturalEarthFeature(\n", + " 'physical',\n", + " 'lakes',\n", + " scale=resolution,\n", + " edgecolor='face',\n", + " facecolor=cfeature.COLORS['water'],\n", + ")\n", + "state_borders = cfeature.NaturalEarthFeature(\n", + " category='cultural',\n", + " name='admin_1_states_provinces_lakes',\n", + " scale=resolution,\n", + " facecolor='none',\n", + ")\n", + "\n", + "ax.add_feature(land_mask)\n", + "ax.add_feature(sea_mask)\n", + "ax.add_feature(lake_mask)\n", + "ax.add_feature(state_borders, linestyle='solid', edgecolor='black')\n", + "ax.set_title('New York and Vicinity; lower resolution');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Finally, let's create a figure with two subplots. On one, we'll repeat our hi-res NYS map; on the second, we'll plot over a different part of the world." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create the figure object\n", + "fig = plt.figure(\n", + " figsize=(30, 24)\n", + ") # Notice we need a bigger \"canvas\" so these two maps will be of a decent size\n", + "\n", + "# First subplot\n", + "ax = plt.subplot(2, 1, 1, projection=projLccNY)\n", + "ax.set_extent([lonW, lonE, latS, latN], crs=projPC)\n", + "ax.add_feature(cfeature.LAND)\n", + "ax.add_feature(cfeature.OCEAN)\n", + "ax.add_feature(cfeature.COASTLINE)\n", + "ax.add_feature(cfeature.BORDERS, linestyle='--')\n", + "ax.add_feature(cfeature.LAKES, alpha=0.5)\n", + "ax.add_feature(cfeature.STATES)\n", + "ax.set_title('New York and Vicinity')\n", + "\n", + "# Set the domain for defining the second plot region.\n", + "latN = 70\n", + "latS = 30.2\n", + "lonW = -10\n", + "lonE = 50\n", + "cLat = (latN + latS) / 2\n", + "cLon = (lonW + lonE) / 2\n", + "\n", + "projLccEur = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)\n", + "\n", + "# Second subplot\n", + "ax2 = plt.subplot(2, 1, 2, projection=projLccEur)\n", + "ax2.set_extent([lonW, lonE, latS, latN], crs=projPC)\n", + "ax2.add_feature(cfeature.LAND)\n", + "ax2.add_feature(cfeature.OCEAN)\n", + "ax2.add_feature(cfeature.COASTLINE)\n", + "ax2.add_feature(cfeature.BORDERS, linestyle='--')\n", + "ax2.add_feature(cfeature.LAKES, alpha=0.5)\n", + "ax2.add_feature(cfeature.STATES)\n", + "ax2.set_title('Europe');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.8.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/environment.yml b/environment.yml index f36fb0544..e632eff59 100644 --- a/environment.yml +++ b/environment.yml @@ -7,5 +7,6 @@ dependencies: - matplotlib - numpy - pre-commit - - xarray + - sqlalchemy<1.4 - cartopy + - xarray diff --git a/pyproject.toml b/pyproject.toml index b05d384bc..d28feaaf1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,6 +10,7 @@ black = "pyproject.toml" [tool.nbqa.mutate] black = 1 pyupgrade = 1 +isort = 1 [tool.nbqa.addopts] pyupgrade = ["--py36-plus"]