{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# NCL_shapefiles_1.py\nThis script illustrates the following concepts:\n - Reading shapefiles\n - Plotting data from shapefiles\n - Using shapefile data to plot unemployment percentages in the U.S.\n - Drawing a custom colorbar on a map\n - Drawing filled polygons over a Lambert Conformal plot\n - Drawing the US with a Lambert Conformal projection\n - Zooming in on a particular area on a Lambert Conformal map\n - Centering the labels under the colorbar boxes\n\nSee following URLs to see the reproduced NCL plot & script:\n - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/shapefiles_1.ncl\n - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/shapefiles_1_lg.png\n\nNote:\n At the time of making this example, there isn't a good way to draw tick \n marks along with the latitude and longitude labels. We have chosen to draw\n gridlines to show exactly where the labels are pointing. The gridlines can\n be removed by calling ``gl.xlines = False`` and ``gl.ylines = False``\n after drawing the labels.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import packages:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nimport matplotlib.colors as colors\nimport matplotlib.cm as cm\nimport matplotlib.ticker as mticker\nimport shapefile as shp\nimport numpy as np\nimport cartopy.crs as ccrs\nimport cartopy.feature as cfeature\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import util as gvutil" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read in data:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Open all shapefiles and associated .dbf, .shp, and .prj files\nopen(gdf.get(\"shape_files/states.dbf\"), 'r')\nopen(gdf.get(\"shape_files/states.shp\"), 'r')\nopen(gdf.get(\"shape_files/states.shx\"), 'r')\nopen(gdf.get(\"shape_files/states.prj\"), 'r')\n\n# Open shapefiles\nshapefile = shp.Reader(gdf.get(\"shape_files/states.dbf\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set color map colors and bounds\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "colormap = colors.ListedColormap(['blue', 'lime', 'yellow', 'red'])\n\ncolorbounds = [0.5, 1.5, 2.5, 3.5, 4.5]\n\nnorm = colors.BoundaryNorm(colorbounds, colormap.N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Helper function to determine state color:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def color_assignment(record):\n population = record.PERSONS\n unempolyment = record.UNEMPLOY\n percent = unempolyment / population\n if (0.01 <= percent and percent < 0.02):\n return colormap.colors[0]\n elif (0.02 <= percent and percent < 0.03):\n return colormap.colors[1]\n elif (0.03 <= percent and percent < 0.04):\n return colormap.colors[2]\n elif (0.04 <= percent):\n return colormap.colors[3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(10, 8))\nax = plt.axes(projection=ccrs.LambertConformal(standard_parallels=(33, 45),\n central_longitude=-98))\nax.set_extent([-125, -74, 22, 50])\n\nax.add_feature(cfeature.LAND, color='silver', zorder=0)\nax.add_feature(cfeature.LAKES, color='white', zorder=1)\n\nfor i in range(0, len(shapefile.shapes())):\n shape = shapefile.shape(i)\n record = shapefile.record(i)\n color = color_assignment(record)\n # if a shape has multiple parts make each one a separate patch\n if len(shape.parts) > 1:\n for j in range(0, len(shape.parts)):\n start_index = shape.parts[j]\n # the last part uses the remaining points and doesn't require and end_index\n if (j is (len(shape.parts) - 1)):\n patch = mpatches.Polygon(shape.points[start_index:],\n facecolor=color,\n edgecolor='black',\n linewidth=0.5,\n transform=ccrs.PlateCarree(),\n zorder=2)\n else:\n end_index = shape.parts[j + 1]\n patch = mpatches.Polygon(shape.points[start_index:end_index],\n facecolor=color,\n edgecolor='black',\n linewidth=0.5,\n transform=ccrs.PlateCarree(),\n zorder=2)\n ax.add_patch(patch)\n else:\n patch = mpatches.Polygon(shape.points,\n facecolor=color,\n edgecolor='black',\n linewidth=0.5,\n transform=ccrs.PlateCarree(),\n zorder=2)\n ax.add_patch(patch)\n\n# Create colorbar\nplt.colorbar(cm.ScalarMappable(cmap=colormap, norm=norm),\n ax=ax,\n boundaries=colorbounds,\n orientation='horizontal',\n shrink=0.75,\n ticks=[1, 2, 3, 4],\n label='percent',\n aspect=30,\n pad=0.075)\n\n# Add latitude and longitude labels\ngl = ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)\ngl.xlocator = mticker.FixedLocator(np.linspace(-120, -80, 5))\ngl.ylocator = mticker.FixedLocator(np.linspace(25, 45, 5))\ngl.xlabel_style = {'rotation': 0}\ngl.ylabel_style = {'rotation': 0}\n\n# Use geocat.viz.util convenience function to set titles and labels\ngvutil.set_titles_and_labels(ax, maintitle='Percentage unemployment, by state')\n\nplt.show()" ] } ], "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.9" } }, "nbformat": 4, "nbformat_minor": 0 }