python-awips: Working with the Maps and Topo Databases
\n",
"
Unidata AMS 2021 Student Conference
\n",
"\n",
"\n",
"
\n",
"\n",
"---\n",
"\n",
"\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
}