mirror of
https://github.com/Unidata/python-awips.git
synced 2025-02-24 14:57:57 -05:00
513 lines
17 KiB
ReStructuredText
513 lines
17 KiB
ReStructuredText
============================
|
||
Map Resources and Topography
|
||
============================
|
||
`Notebook <http://nbviewer.ipython.org/github/Unidata/python-awips/blob/master/examples/notebooks/Map_Resources_and_Topography.ipynb>`_
|
||
Python-AWIPS Tutorial Notebook
|
||
|
||
--------------
|
||
|
||
Objectives
|
||
==========
|
||
|
||
- Use python-awips to connect to an edex server
|
||
- Define data request object specifically for the maps database
|
||
- Manipulate request object for various different map resources
|
||
- Plot map resources in combination with one another for geographical
|
||
context
|
||
|
||
--------------
|
||
|
||
Table of Contents
|
||
-----------------
|
||
|
||
| `1
|
||
Imports <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#imports>`__\
|
||
| `2 Connect to
|
||
EDEX <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#connect-to-edex>`__\
|
||
| `3 Function:
|
||
make_map() <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#function-make-map>`__\
|
||
| `4 Create Initial Map From
|
||
CWA <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#create-initial-map-from-cwa>`__\
|
||
| `5 Draw Merged
|
||
CWA <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#draw-merged-cwa>`__\
|
||
| `6 Draw Interstates using Boundary
|
||
Filter <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#draw-interstates-using-boundary-filter>`__\
|
||
| `7 Draw Nearby
|
||
Cities <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#draw-nearby-cities>`__\
|
||
| `8 Draw Nearby
|
||
Lakes <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#draw-nearby-lakes>`__\
|
||
| `9 Draw Major
|
||
Rivers <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#draw-major-rivers>`__\
|
||
| `10 Draw
|
||
Topography <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#draw-topography>`__\
|
||
| `11 See
|
||
Also <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#see-also>`__\
|
||
| `11.1 Additional
|
||
Documentation <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html#additional-documentation>`__\
|
||
|
||
1 Imports
|
||
---------
|
||
|
||
The imports below are used throughout the notebook. Note the first
|
||
import is coming directly from python-awips and allows us to connect to
|
||
an EDEX server. The subsequent imports are for data manipulation and
|
||
visualization.
|
||
|
||
.. code:: ipython3
|
||
|
||
from awips.dataaccess import DataAccessLayer
|
||
import matplotlib.pyplot as plt
|
||
import cartopy.crs as ccrs
|
||
import numpy.ma as ma
|
||
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
|
||
from cartopy.feature import ShapelyFeature,NaturalEarthFeature
|
||
from shapely.ops import cascaded_union
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
2 Connect to EDEX
|
||
-----------------
|
||
|
||
First we establish a connection to Unidata’s public EDEX server. With
|
||
that connection made, we can create a `new data request
|
||
object <http://unidata.github.io/python-awips/api/IDataRequest.html>`__
|
||
and set the data type to **maps**.
|
||
|
||
.. code:: ipython3
|
||
|
||
# Server, Data Request Type, and Database Table
|
||
DataAccessLayer.changeEDEXHost("edex-cloud.unidata.ucar.edu")
|
||
request = DataAccessLayer.newDataRequest('maps')
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
3 Function: make_map()
|
||
----------------------
|
||
|
||
In many of our notebooks we end up plotting map images, and this logic
|
||
below is the same from those other notebooks. Typically, functions are
|
||
defined when they are called multiple times throughout a notebook. In
|
||
this case, we only use it in one code block cell, but because it is a
|
||
common function from several of our notebooks, it’s nice to keep the
|
||
logic neatly defined for consistency.
|
||
|
||
.. code:: ipython3
|
||
|
||
# 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
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
4 Create Initial Map From CWA
|
||
-----------------------------
|
||
|
||
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 and can be easily plotted by many Python 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.
|
||
|
||
..
|
||
|
||
**Tip**: Note the geometry definition of ``the_geom`` for each data
|
||
type, which can be **Point**, **MultiPolygon**, or
|
||
**MultiLineString**.
|
||
|
||
Here we’ll be using Boulder (BOU) as our example for plotting the County
|
||
Warning Area (CWA). We’ll query our EDEX server to get all counties in
|
||
the CWA for BOU, and then plot those counties along withe the state
|
||
boundaries and lines of longitude and latitude. In order to get this
|
||
information from EDEX, we’ll need to set several characteristics on our
|
||
data request object. We will use
|
||
`request.setParameters() <http://unidata.github.io/python-awips/api/IDataRequest.html#awips.dataaccess.IDataRequest.setParameters>`__
|
||
to refine our query to EDEX.
|
||
|
||
.. code:: ipython3
|
||
|
||
# Specify the necessary identifiers for requesting the Boulder CWA
|
||
request.addIdentifier('table', 'mapdata.county')
|
||
# 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')
|
||
|
||
# Get response and create dict of county geometries
|
||
response = DataAccessLayer.getGeometryData(request)
|
||
counties = []
|
||
for ob in response:
|
||
counties.append(ob.getGeometry())
|
||
print("Using " + str(len(counties)) + " county MultiPolygons")
|
||
|
||
|
||
# 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]
|
||
|
||
# Create the map we'll use for the rest of this notebook based on the
|
||
# boundaries of the CWA
|
||
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
|
||
shape_feature = ShapelyFeature(counties,ccrs.PlateCarree(),
|
||
facecolor='none', linestyle="-",edgecolor='#86989B')
|
||
ax.add_feature(shape_feature)
|
||
|
||
|
||
.. parsed-literal::
|
||
|
||
Using 22 county MultiPolygons
|
||
|
||
|
||
|
||
|
||
.. parsed-literal::
|
||
|
||
<cartopy.mpl.feature_artist.FeatureArtist at 0x11568f6d0>
|
||
|
||
|
||
|
||
|
||
.. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_13_2.png
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
5 Draw Merged CWA
|
||
-----------------
|
||
|
||
In the previous section we created a merged polygon with the applicable
|
||
counties. Here, we draw this new shape on top of our existing map in a
|
||
burnt orange color.
|
||
|
||
.. code:: ipython3
|
||
|
||
# Plot CWA envelope
|
||
shape_feature = ShapelyFeature(boundaries,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_16_0.png
|
||
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
6 Draw Interstates using Boundary Filter
|
||
----------------------------------------
|
||
|
||
Now, we’ll use the previously-defined
|
||
**envelope=merged_counties.buffer(2)** in a **newDataRequest()** to
|
||
request interstate geometries which fall inside the buffered boundary.
|
||
|
||
.. code:: ipython3
|
||
|
||
# Define the request for the interstate query
|
||
request = DataAccessLayer.newDataRequest('maps', envelope=envelope)
|
||
request.addIdentifier('table', 'mapdata.interstate')
|
||
request.addIdentifier('geomField', 'the_geom')
|
||
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_19_1.png
|
||
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
7 Draw Nearby Cities
|
||
--------------------
|
||
|
||
Request the city table based using the **envelope** attribute 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
|
||
|
||
# Define the request for the city query
|
||
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")
|
||
|
||
# Set aside two arrays - one for the geometry of the cities and one for their names
|
||
citylist = []
|
||
cityname = []
|
||
# For BOU, progressive disclosure values above 50 and pop above 5000 looks good
|
||
for ob in cities:
|
||
if ob.getString("population") != 'None':
|
||
if ob.getNumber("prog_disc") > 50 and 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 1205 total cities
|
||
Plotting 58 cities
|
||
|
||
|
||
|
||
|
||
.. image:: Map_Resources_and_Topography_files/Map_Resources_and_Topography_22_1.png
|
||
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
8 Draw Nearby Lakes
|
||
-------------------
|
||
|
||
Again, use the **envelope** attribute to define a new data requst for
|
||
the nearby lakes.
|
||
|
||
.. code:: ipython3
|
||
|
||
# Define request for lakes
|
||
request = DataAccessLayer.newDataRequest('maps', envelope=envelope)
|
||
request.addIdentifier('table', 'mapdata.lake')
|
||
request.addIdentifier('geomField', 'the_geom')
|
||
|
||
# Get lake geometries
|
||
response = DataAccessLayer.getGeometryData(request)
|
||
print("Using " + str(len(response)) + " lake MultiPolygons")
|
||
|
||
# Plot lakes
|
||
shape_feature = ShapelyFeature([lake.getGeometry() for lake in response],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_25_1.png
|
||
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
9 Draw Major Rivers
|
||
-------------------
|
||
|
||
.. code:: ipython3
|
||
|
||
# Define request for rivers
|
||
request = DataAccessLayer.newDataRequest('maps', envelope=envelope)
|
||
request.addIdentifier('table', 'mapdata.majorrivers')
|
||
request.addIdentifier('geomField', 'the_geom')
|
||
rivers = DataAccessLayer.getGeometryData(request)
|
||
print("Using " + str(len(rivers)) + " river MultiLineStrings")
|
||
|
||
# Plot rivers
|
||
shape_feature = ShapelyFeature([river.getGeometry() for river in rivers],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_28_1.png
|
||
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
10 Draw Topography
|
||
------------------
|
||
|
||
Spatial envelopes are required for topo requests, which can become slow
|
||
to download and render for large (CONUS) maps.
|
||
|
||
.. code:: ipython3
|
||
|
||
# Define topography request
|
||
request = DataAccessLayer.newDataRequest("topo", envelope=envelope)
|
||
request.addIdentifier("group", "/")
|
||
request.addIdentifier("dataset", "full")
|
||
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::
|
||
|
||
[<awips.dataaccess.PyGridData.PyGridData object at 0x115a20370>]
|
||
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_32_1.png
|
||
|
||
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|
||
|
||
11 See Also
|
||
-----------
|
||
|
||
11.1 Additional Documentation
|
||
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||
|
||
- This notebook requires: **python-awips, numpy, matplotplib, cartopy,
|
||
shapely**
|
||
- Use datatype **maps** and **addIdentifier(‘table’, <postgres maps
|
||
schema>)** 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 <http://unidata.github.io/awips2/python/maps-database/#mapdatacw>`__
|
||
for available database tables, column names, and types.
|
||
|
||
`Top <https://unidata.github.io/python-awips/examples/generated/Map_Resources_and_Topography.html>`__
|
||
|
||
--------------
|