Commit a28f142a authored by Jon Pye's avatar Jon Pye
Browse files

Added demo-related comments and alternate commands and some extra features to pretty thngs up.

parent 9933468b
# --------------
# Accessing GeoServer through Python's "owslib"
# --------------
# We can get data from our Geoserver using plain Python code.
# It can be used with other Python modules to generate a pretty good looking map.
# Imports to get the data out of GeoServer
from owslib.wfs import WebFeatureService
import geojson
import numpy as np
# Imports to plot the map
# Imports to plot the data onto a map
import as ccrs
from import MapQuestOpenAerial, MapQuestOSM, OSM
import cartopy.feature as feature
import matplotlib.pyplot as plt
from shapely.geometry import Point
def flip_geojson_coordinates(geo):
......@@ -24,44 +31,64 @@ def flip_geojson_coordinates(geo):
for k in geo:
## ======
## Connect to the GeoServer Web Feature Service. This will serve up sets of vectors, called layers
## They're sets of points usually in our case.
## ======
wfs = WebFeatureService('', version="2.0.0")
# Print the name of every layer available on this Geoserver
print list(wfs.contents)
# Define the bounding box for the map using the stations.
a = wfs.contents['otn:stations']
b = a.boundingBoxWGS84
# Get a list of all stations. There is a lot of metadata associated with this JSON data.
# But mostly we're interested in the coordinates.
json_response = wfs.getfeature(typename=['otn:stations'], propertyname=None, outputFormat='application/json').read()
geom = geojson.loads(json_response)
print geom.keys()
# geom contains a collection of features, and their default coordinate reference system is a stored property.
print geom['type'], geom['crs'] # check our coordinate reference system
geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
# Open Source Imagery from MapQuest (max zoom = 16?)
tiler = MapQuestOpenAerial()
ax = plt.axes(projection=ccrs.PlateCarree())
extent = (b[0]-0.1*dx,b[2]+0.1*dx,b[1]-0.1*dy,b[3]+0.1*dy)
ax.set_extent(extent, geodetic)
# Get a list of all features.
locs = geom.get("features", [])
lons, lats = zip(*[(x["geometry"]["coordinates"][0], x["geometry"]["coordinates"][1]) for x in locs])
# Can get other data here, too.
# Each feature has a lot of reference data associated with, as defined in Geoserver.
print locs[0]['properties'].keys()
print locs[0]['properties'].values()
# But we also have the coordinates, which is what gets us onto the map.
#lons, lats = zip(*[(x["geometry"]["coordinates"][0], x["geometry"]["coordinates"][1]) for x in locs])
# You could also subselect if that fits your project.
lons, lats = zip(*[(x["geometry"]["coordinates"][0], x["geometry"]["coordinates"][1]) for x in locs if x['properties']['collectioncode']=="CBS"])
# Plot the lons and lats, transforming them to the axis projection
# We get around a lot of projection-related unpleasantness by dealing exclusively in point data here.
ax.scatter(lons, lats, transform=ccrs.PlateCarree())
print 'done with geometry'
ax.stock_img() # fairly nice looking bathymetry provided by Cartopy
# put up some nice looking grid lines
gl.xlabels_top = False
gl.ylabels_right = False
# Add the land mass colorings
ax.add_feature(feature.NaturalEarthFeature('physical', 'land', '10m',
edgecolor='face', facecolor=feature.COLORS['land']))
# You can also define a feature and add it in a separate command
states_provinces = feature.NaturalEarthFeature('cultural', 'admin_1_states_provinces_lines', scale='50m', facecolor='none')
ax.add_feature(states_provinces, edgecolor='gray')
# There's also a very simple method to call to draw coastlines on your map.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment