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

python-awips: Working with the Maps and Topo Databases

\n", "

Unidata AMS 2021 Student Conference

\n", "\n", "
\n", "
\n", "\n", "---\n", "\n", "
\"Map
\n", "\n", "\n", "### Focuses\n", "\n", "* Use the AWIPS Maps Database to access GIS objects which are returned as Shapely geometries (*Polygon*, *Point*, *MultiLineString*, etc.) and can be easily plotted by Matplotlib, Cartopy, MetPy, and other packages. \n", "* Use **maps** and **topo** data types to obtain GIS data from the AWIPS databases.\n", "* Note the geometry definition of `the_geom` for each data type, which can be **Point**, **MultiPolygon**, or **MultiLineString**.\n", "* Step through how to plot several layers of data onto an image, including: county and state boundaries, CWA boundaries, interstates, cities, lakes, rivers, and topograpy.\n", "\n", "\n", "\n", "### Objectives\n", "\n", "1. [Define Map Data Request](#1.-Define-the-Map-Data-Request)\n", "1. [Define Map Properties](#2.-Define-Map-Properties)\n", "1. [Draw County and State Boundaries](#3.-Draw-County-and-State-Boundaries)\n", "1. [Draw CWA Boundary](#4.-Draw-CWA-Boundary)\n", "1. [Draw Interstates](#5.-Draw-Interstates)\n", "1. [Draw Nearby Cities](#6.-Draw-Nearby-Cities)\n", "1. [Draw Lakes](#7.-Draw-Lakes)\n", "1. [Draw Major Rivers](#8.-Draw-Major-Rivers)\n", "1. [Draw Topography](#9.-Draw-Topography)\n", "\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from __future__ import print_function\n", "from awips.dataaccess import DataAccessLayer\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", "import numpy as np\n", "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n", "from cartopy.feature import ShapelyFeature,NaturalEarthFeature\n", "from shapely.geometry import Polygon\n", "from shapely.ops import cascaded_union\n", "import numpy.ma as ma" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Define the Map Data Request\n", "\n", "If you read through the [python-awips: How to Access Data](https://nbviewer.jupyter.org/github/Unidata/pyaos-ams-2021/blob/master/notebooks/dataAccess/python-awips-HowToAccessData.ipynb) training, you will know that we need to set an EDEX url to access our server, and then we create a data request. In this example we use *maps* as the data type to define our request. We'll use Boulder as our location, so set the *Location Name* on the request to BDU. Then add several *Identifiers* for various fields of interest." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create EDEX data request\n", "DataAccessLayer.changeEDEXHost(\"edex-cloud.unidata.ucar.edu\")\n", "request = DataAccessLayer.newDataRequest('maps')\n", "request.addIdentifier('table', 'mapdata.county')\n", "\n", "# Define a WFO ID for location\n", "# tie this ID to the mapdata.county column \"cwa\" for filtering\n", "request.setLocationNames('BOU')\n", "request.addIdentifier('cwa', 'BOU')\n", "\n", "# enable location filtering (inLocation)\n", "# locationField is tied to the above cwa definition (BOU)\n", "request.addIdentifier('geomField', 'the_geom')\n", "request.addIdentifier('inLocation', 'true')\n", "request.addIdentifier('locationField', 'cwa')\n", "\n", "# take a look at the request\n", "print(request)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 2. Define Map Properties\n", "\n", "If more than one plot is drawn, then it's easiest to define common logic in a function. Here, a new function called **make_map** is defined. This function uses the [matplotlib.pyplot package (plt)](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.html) to create a figure and axis. The coastlines (continental boundaries) are added, along with lat/lon grids. \n", "\n", "> Note: We only use this function once in this notebook, but it's in here as an example of how to write a function and use it, for future reference." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Standard map plot\n", "def make_map(bbox, projection=ccrs.PlateCarree()):\n", " fig, ax = plt.subplots(figsize=(12,12),\n", " subplot_kw=dict(projection=projection))\n", " ax.set_extent(bbox)\n", " ax.coastlines(resolution='50m')\n", " gl = ax.gridlines(draw_labels=True)\n", " gl.top_labels = gl.right_labels = False\n", " gl.xformatter = LONGITUDE_FORMATTER\n", " gl.yformatter = LATITUDE_FORMATTER\n", " return fig, ax" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 3. Draw County and State Boundaries\n", "\n", "Start by getting the [GeometryData](http://unidata.github.io/python-awips/api/PyGeometryData.html) back from EDEX for the map request. We'll create a plot of the county boundaries for the WFO (in this case Boulder - BDU). Add in the " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get response and create dict of county geometries\n", "response = DataAccessLayer.getGeometryData(request, [])\n", "\n", "# Populate an array of the counties for BDU\n", "counties = np.array([])\n", "for ob in response:\n", " counties = np.append(counties,ob.getGeometry())\n", "\n", "\n", "# All WFO counties merged to a single Polygon\n", "merged_counties = cascaded_union(counties)\n", "envelope = merged_counties.buffer(2)\n", "boundaries=[merged_counties]\n", "\n", "# Get bounds of this merged Polygon to use as buffered map extent\n", "bounds = merged_counties.bounds\n", "bbox=[bounds[0]-1,bounds[2]+1,bounds[1]-1.5,bounds[3]+1.5]\n", "\n", "# Create the map using our defined function\n", "fig, ax = make_map(bbox=bbox)\n", "\n", "# Plot state boundaries handled by Cartopy\n", "states = NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none')\n", "ax.add_feature(states, linestyle='-', edgecolor='black',linewidth=2)\n", "\n", "# Plot CWA counties\n", "for i, geom in enumerate(counties):\n", " cbounds = Polygon(geom)\n", " intersection = cbounds.intersection\n", " geoms = (intersection(geom) for geom in counties if cbounds.intersects(geom))\n", " shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), facecolor='none', linestyle=\"-\",edgecolor='#86989B')\n", " ax.add_feature(shape_feature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 4. Draw CWA Boundary\n", "\n", "Use the single polygon that encompases all of the counties in the CWA for Boulder to draw that on top of the previous figure." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# use the merged polygon to draw the CWA boundary\n", "geom = boundaries[0]\n", "gbounds = Polygon(geom)\n", "intersection = gbounds.intersection\n", "geoms = (intersection(geom) for geom in boundaries if gbounds.intersects(geom))\n", "shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), facecolor='none', linestyle=\"-\",linewidth=3.,edgecolor='#cc5000')\n", "ax.add_feature(shape_feature)\n", "\n", "# Display the plot\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 5. Draw Interstates\n", "\n", "Using the previously-defined **envelope=merged_counties.buffer(2)** in **newDataRequest()** to request geometries which fall inside the buffered boundary. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create new request to get interstate polygons from the EDEX server\n", "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n", "request.addIdentifier('table', 'mapdata.interstate')\n", "request.addIdentifier('geomField', 'the_geom')\n", "request.setParameters('name')\n", "interstates = DataAccessLayer.getGeometryData(request, [])\n", "\n", "# Plot interstates\n", "for ob in interstates:\n", " shape_feature = ShapelyFeature(ob.getGeometry(),ccrs.PlateCarree(), facecolor='none', linestyle=\"-\",edgecolor='orange')\n", " ax.add_feature(shape_feature)\n", "\n", "# Display the plot\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 6. Draw Nearby Cities\n", "\n", "Request the city table and filter by population and progressive disclosure level:\n", "\n", "> **Warning**: the `prog_disc` field is not entirely understood and values appear to change significantly depending on WFO site. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create new request for local cities\n", "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n", "request.addIdentifier('table', 'mapdata.city')\n", "request.addIdentifier('geomField', 'the_geom')\n", "request.setParameters('name','population','prog_disc')\n", "cities = DataAccessLayer.getGeometryData(request, [])\n", "\n", "citylist = []\n", "cityname = []\n", "# For BOU, progressive disclosure values above 50 and pop above 5000 looks good\n", "for ob in cities:\n", " if ob.getString(\"population\"):\n", " if ob.getNumber(\"prog_disc\") > 50:\n", " if int(ob.getString(\"population\")) > 5000:\n", " citylist.append(ob.getGeometry())\n", " cityname.append(ob.getString(\"name\"))\n", "\n", "# Plot city markers\n", "ax.scatter([point.x for point in citylist], [point.y for point in citylist], transform=ccrs.PlateCarree(),marker=\"+\",facecolor='black')\n", "# Plot city names\n", "for i, txt in enumerate(cityname):\n", " ax.annotate(txt, (citylist[i].x,citylist[i].y), xytext=(3,3), textcoords=\"offset points\")\n", "\n", "# Display the plot\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 7. Draw Lakes\n", "\n", "Get lake geometry data from the map database from the EDEX server and add it to the plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Create a request for lakes\n", "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n", "request.addIdentifier('table', 'mapdata.lake')\n", "request.addIdentifier('geomField', 'the_geom')\n", "request.setParameters('name')\n", "\n", "# Get lake geometries\n", "response = DataAccessLayer.getGeometryData(request, [])\n", "lakes = np.array([])\n", "for ob in response:\n", " lakes = np.append(lakes,ob.getGeometry())\n", "\n", "# Plot lakes\n", "for i, geom in enumerate(lakes):\n", " cbounds = Polygon(geom)\n", " intersection = cbounds.intersection\n", " geoms = (intersection(geom) for geom in lakes if cbounds.intersects(geom))\n", " shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), facecolor='blue', linestyle=\"-\",edgecolor='#20B2AA')\n", " ax.add_feature(shape_feature)\n", "\n", "# Display the plot\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 8. Draw Major Rivers" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a new request for major rivers\n", "request = DataAccessLayer.newDataRequest('maps', envelope=envelope)\n", "request.addIdentifier('table', 'mapdata.majorrivers')\n", "request.addIdentifier('geomField', 'the_geom')\n", "request.setParameters('pname')\n", "rivers = DataAccessLayer.getGeometryData(request, [])\n", "\n", "# Plot rivers\n", "for ob in rivers:\n", " shape_feature = ShapelyFeature(ob.getGeometry(),ccrs.PlateCarree(), facecolor='none', linestyle=\":\", edgecolor='#20B2AA')\n", " ax.add_feature(shape_feature)\n", "\n", "# Display the plot \n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## 9. Draw Topography\n", "\n", "Spatial envelopes are required for topo requests, which can become slow to download and render for large (CONUS) maps." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Create new request for l\n", "request = DataAccessLayer.newDataRequest(\"topo\")\n", "request.addIdentifier(\"group\", \"/\")\n", "request.addIdentifier(\"dataset\", \"full\")\n", "request.setEnvelope(envelope)\n", "gridData = DataAccessLayer.getGridData(request)\n", "grid=gridData[0]\n", "topo=ma.masked_invalid(grid.getRawData()) \n", "lons, lats = grid.getLatLonCoords()\n", "# print(topo.min()) # minimum elevation in our domain (meters)\n", "# print(topo.max()) # maximum elevation in our domain (meters)\n", "\n", "# Plot topography\n", "cs = ax.contourf(lons, lats, topo, 80, cmap=plt.get_cmap('terrain'),alpha=0.1, extend='both')\n", "cbar = fig.colorbar(cs, shrink=0.5, orientation='horizontal')\n", "cbar.set_label(\"topography height in meters\")\n", "\n", "# Display the plot\n", "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Top\n", "\n", "---\n", "\n", "## See also\n", "\n", "Documentation for:\n", "\n", "* [AWIPS Maps Database Reference List](http://unidata.github.io/awips2/python/maps-database/#mapdatacwa)\n", "* [awips.DataAccessLayer](http://unidata.github.io/python-awips/api/DataAccessLayer.html)\n", "* [awips.PyGeometryData](http://unidata.github.io/python-awips/api/PyGeometryData.html)\n", "* [matplotlib.pyplot](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.html)\n", "* [matplotlib.pyplot.subplot](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.subplot.html)\n", "* [matplotlib.pyplot.contourf](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.contourf.html)\n", "\n", "\n", "### Related Notebooks:\n", "\n", "* [python-awips: How to Access Data](https://nbviewer.jupyter.org/github/Unidata/pyaos-ams-2021/blob/master/notebooks/dataAccess/python-awips-HowToAccessData.ipynb)\n", "\n", "\n", "Top" ] } ], "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": 1 }