============================ Map Resources and Topography ============================ `Notebook `_ The python-awips package provides access to the entire AWIPS Maps Database for use in Python GIS applications. Map objects are returned as Shapely geometries (*Polygon*, *Point*, *MultiLineString*, etc.) and can be easily plotted by Matplotlib, Cartopy, MetPy, and other packages. Each map database table has a geometry field called ``the_geom``, which can be used to spatially select map resources for any column of type geometry, Notes ----- - This notebook requires: **python-awips, numpy, matplotplib, cartopy, shapely** - Use datatype **maps** and **addIdentifier('table', )** to define the map table: DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu") request = DataAccessLayer.newDataRequest('maps') request.addIdentifier('table', 'mapdata.county') - Use **request.setLocationNames()** and **request.addIdentifier()** to spatially filter a map resource. In the example below, WFO ID **BOU** (Boulder, Colorado) is used to query counties within the BOU county watch area (CWA) :: request.addIdentifier('geomField', 'the_geom') request.addIdentifier('inLocation', 'true') request.addIdentifier('locationField', 'cwa') request.setLocationNames('BOU') request.addIdentifier('cwa', 'BOU') See the Maps Database Reference Page for available database tables, column names, and types. Note the geometry definition of ``the_geom`` for each data type, which can be **Point**, **MultiPolygon**, or **MultiLineString**. Setup ----- .. code:: ipython3 from __future__ import print_function from awips.dataaccess import DataAccessLayer import matplotlib.pyplot as plt import cartopy.crs as ccrs import numpy as np from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER from cartopy.feature import ShapelyFeature,NaturalEarthFeature from shapely.geometry import Polygon from shapely.ops import cascaded_union # Standard map plot def make_map(bbox, projection=ccrs.PlateCarree()): fig, ax = plt.subplots(figsize=(12,12), subplot_kw=dict(projection=projection)) ax.set_extent(bbox) ax.coastlines(resolution='50m') gl = ax.gridlines(draw_labels=True) gl.top_labels = gl.right_labels = False gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER return fig, ax # Server, Data Request Type, and Database Table DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu") request = DataAccessLayer.newDataRequest('maps') request.addIdentifier('table', 'mapdata.county') Request County Boundaries for a WFO ----------------------------------- - Use **request.setParameters()** to define fields to be returned by the request. .. code:: ipython3 # Define a WFO ID for location # tie this ID to the mapdata.county column "cwa" for filtering request.setLocationNames('BOU') request.addIdentifier('cwa', 'BOU') # enable location filtering (inLocation) # locationField is tied to the above cwa definition (BOU) request.addIdentifier('geomField', 'the_geom') request.addIdentifier('inLocation', 'true') request.addIdentifier('locationField', 'cwa') # This is essentially the same as "'"select count(*) from mapdata.cwa where cwa='BOU';" (=1) # Get response and create dict of county geometries response = DataAccessLayer.getGeometryData(request, []) counties = np.array([]) for ob in response: counties = np.append(counties,ob.getGeometry()) print("Using " + str(len(counties)) + " county MultiPolygons") %matplotlib inline # All WFO counties merged to a single Polygon merged_counties = cascaded_union(counties) envelope = merged_counties.buffer(2) boundaries=[merged_counties] # Get bounds of this merged Polygon to use as buffered map extent bounds = merged_counties.bounds bbox=[bounds[0]-1,bounds[2]+1,bounds[1]-1.5,bounds[3]+1.5] fig, ax = make_map(bbox=bbox) # Plot political/state boundaries handled by Cartopy political_boundaries = NaturalEarthFeature(category='cultural', name='admin_0_boundary_lines_land', scale='50m', facecolor='none') states = NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none') ax.add_feature(political_boundaries, linestyle='-', edgecolor='black') ax.add_feature(states, linestyle='-', edgecolor='black',linewidth=2) # Plot CWA counties for i, geom in enumerate(counties): cbounds = Polygon(geom) intersection = cbounds.intersection geoms = (intersection(geom) for geom in counties if cbounds.intersects(geom)) shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), facecolor='none', linestyle="-",edgecolor='#86989B') ax.add_feature(shape_feature) .. parsed-literal:: Using 23 county MultiPolygons .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_4_1.png Create a merged CWA with cascaded\_union ---------------------------------------- .. code:: ipython3 # Plot CWA envelope for i, geom in enumerate(boundaries): gbounds = Polygon(geom) intersection = gbounds.intersection geoms = (intersection(geom) for geom in boundaries if gbounds.intersects(geom)) shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), facecolor='none', linestyle="-",linewidth=3.,edgecolor='#cc5000') ax.add_feature(shape_feature) fig .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_6_0.png WFO boundary spatial filter for interstates ------------------------------------------- Using the previously-defined **envelope=merged\_counties.buffer(2)** in **newDataRequest()** to request geometries which fall inside the buffered boundary. .. code:: ipython3 request = DataAccessLayer.newDataRequest('maps', envelope=envelope) request.addIdentifier('table', 'mapdata.interstate') request.addIdentifier('geomField', 'the_geom') request.setParameters('name') interstates = DataAccessLayer.getGeometryData(request, []) print("Using " + str(len(interstates)) + " interstate MultiLineStrings") # Plot interstates for ob in interstates: shape_feature = ShapelyFeature(ob.getGeometry(),ccrs.PlateCarree(), facecolor='none', linestyle="-",edgecolor='orange') ax.add_feature(shape_feature) fig .. parsed-literal:: Using 225 interstate MultiLineStrings .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_8_1.png Nearby cities ------------- Request the city table and filter by population and progressive disclosure level: **Warning**: the ``prog_disc`` field is not entirely understood and values appear to change significantly depending on WFO site. .. code:: ipython3 request = DataAccessLayer.newDataRequest('maps', envelope=envelope) request.addIdentifier('table', 'mapdata.city') request.addIdentifier('geomField', 'the_geom') request.setParameters('name','population','prog_disc') cities = DataAccessLayer.getGeometryData(request, []) print("Queried " + str(len(cities)) + " total cities") citylist = [] cityname = [] # For BOU, progressive disclosure values above 50 and pop above 5000 looks good for ob in cities: if ob.getString("population"): if ob.getNumber("prog_disc") > 50: if int(ob.getString("population")) > 5000: citylist.append(ob.getGeometry()) cityname.append(ob.getString("name")) print("Plotting " + str(len(cityname)) + " cities") # Plot city markers ax.scatter([point.x for point in citylist], [point.y for point in citylist], transform=ccrs.PlateCarree(),marker="+",facecolor='black') # Plot city names for i, txt in enumerate(cityname): ax.annotate(txt, (citylist[i].x,citylist[i].y), xytext=(3,3), textcoords="offset points") fig .. parsed-literal:: Queried 1203 total cities Plotting 57 cities .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_10_1.png Lakes ----- .. code:: ipython3 request = DataAccessLayer.newDataRequest('maps', envelope=envelope) request.addIdentifier('table', 'mapdata.lake') request.addIdentifier('geomField', 'the_geom') request.setParameters('name') # Get lake geometries response = DataAccessLayer.getGeometryData(request, []) lakes = np.array([]) for ob in response: lakes = np.append(lakes,ob.getGeometry()) print("Using " + str(len(lakes)) + " lake MultiPolygons") # Plot lakes for i, geom in enumerate(lakes): cbounds = Polygon(geom) intersection = cbounds.intersection geoms = (intersection(geom) for geom in lakes if cbounds.intersects(geom)) shape_feature = ShapelyFeature(geoms,ccrs.PlateCarree(), facecolor='blue', linestyle="-",edgecolor='#20B2AA') ax.add_feature(shape_feature) fig .. parsed-literal:: Using 208 lake MultiPolygons .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_12_1.png Major Rivers ------------ .. code:: ipython3 request = DataAccessLayer.newDataRequest('maps', envelope=envelope) request.addIdentifier('table', 'mapdata.majorrivers') request.addIdentifier('geomField', 'the_geom') request.setParameters('pname') rivers = DataAccessLayer.getGeometryData(request, []) print("Using " + str(len(rivers)) + " river MultiLineStrings") # Plot rivers for ob in rivers: shape_feature = ShapelyFeature(ob.getGeometry(),ccrs.PlateCarree(), facecolor='none', linestyle=":",edgecolor='#20B2AA') ax.add_feature(shape_feature) fig .. parsed-literal:: Using 1400 river MultiLineStrings .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_14_1.png Topography ---------- Spatial envelopes are required for topo requests, which can become slow to download and render for large (CONUS) maps. .. code:: ipython3 import numpy.ma as ma request = DataAccessLayer.newDataRequest("topo") request.addIdentifier("group", "/") request.addIdentifier("dataset", "full") request.setEnvelope(envelope) gridData = DataAccessLayer.getGridData(request) print(gridData) print("Number of grid records: " + str(len(gridData))) print("Sample grid data shape:\n" + str(gridData[0].getRawData().shape) + "\n") print("Sample grid data:\n" + str(gridData[0].getRawData()) + "\n") .. parsed-literal:: [] Number of grid records: 1 Sample grid data shape: (778, 1058) Sample grid data: [[1694. 1693. 1688. ... 757. 761. 762.] [1701. 1701. 1701. ... 758. 760. 762.] [1703. 1703. 1703. ... 760. 761. 762.] ... [1767. 1741. 1706. ... 769. 762. 768.] [1767. 1746. 1716. ... 775. 765. 761.] [1781. 1753. 1730. ... 766. 762. 759.]] .. code:: ipython3 grid=gridData[0] topo=ma.masked_invalid(grid.getRawData()) lons, lats = grid.getLatLonCoords() print(topo.min()) # minimum elevation in our domain (meters) print(topo.max()) # maximum elevation in our domain (meters) # Plot topography cs = ax.contourf(lons, lats, topo, 80, cmap=plt.get_cmap('terrain'),alpha=0.1, extend='both') cbar = fig.colorbar(cs, shrink=0.5, orientation='horizontal') cbar.set_label("topography height in meters") fig .. parsed-literal:: 623.0 4328.0 .. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_17_1.png