Open Geospatial Consortium (OGC) Services¶
This is an interactive jupyter notebook that shows users how to access OGC services using owslib.
OWSLib is a Python package for client programming with Open Geospatial Consortium (OGC) web service interface standards, including WMS, WFS, WMTS, WCS, and SOS (https://geopython.github.io/OWSLib/).
Web Map Service¶
Web Map Service (WMS) is an OGC (Open GIS Consortium) standard. There are three major operations are defined for creating and displaying map images, including:
1) GetCapabilities. It is used for retrieving metadata on the service level (Required).
2) GetMap. It is the key core operation of WMS. It is designed for retrieving a map image with the geospatial parameters and the size (Required).
3) GetFeatureInfo. It is designed for retrieving information about certain special features displayed on a map (Optional).
Setup library environment & Create Web Map Service object¶
1) import the owslib.wms library
2) create the Web Map Service Object. The avialable layers can be shown when the wms is valid.
# load owslib library
from owslib.wms import WebMapService
# Create your WebMapService object
wms = WebMapService('http://apps.ecmwf.int/wms/?token=public', version='1.1.1') # version 1.3.0 works as well
#Show the number of available layers
print (len(wms.contents))
Plot the WebMapService from OGC¶
3) In the get map function, we shoud define the map size, the projection information, the bounding box of the map, map format, and the parameter for transparent
%matplotlib inline
import os, sys
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
def getMap(layerName,bbox,filename):
wms.getOperationByName('GetMap').formatOptions
img = wms.getmap(layers=[layerName],
size=(600,300),
srs='EPSG:4326',
bbox=bbox,
format='image/png',
transparent=True)
tmpfile = open(filename,'wb')
tmpfile.write(img.read())
tmpfile.close()
getMap('foreground',(-180,-90,180,90), 'foreground.png')
getMap('background',(-180,-90,180,90), 'background.png')
getMap('composition_bbaod550',(-180,-90,180,90), 'bbaod550.png')
image_back = mpimg.imread('background.png')
image_compos = mpimg.imread('bbaod550.png')
image_fore = mpimg.imread('foreground.png')
fig = plt.figure(figsize=(12,7))
img_back = plt.imshow(image_back,extent=[-180,180,-90,90],aspect='auto')
img_compos = plt.imshow(image_compos,extent=[-180,180,-90,90],aspect='auto')
img_fore = plt.imshow(image_fore,extent=[-180,180,-90,90],aspect='auto')
Web Feature Service¶
Web Feature Service (WFS) is an OGC (Open Geospatial Consortium) which based on GML (Geography Markup Language).
1)import liarary
2) define WebFeatureService object to connect GeoServer WFS service
%matplotlib inline
from owslib.wfs import WebFeatureService
import pandas as pd
import geopandas as gpd
from geopandas import GeoSeries, GeoDataFrame
import requests
import geojson
3) get wfs layer from defined URL
4) print information for the first 10 layers, incluing layer id, title, bbox information
wfs_url = "http://data.nanoos.org/geoserver/ows"
wfs = WebFeatureService(wfs_url, version='1.0.0')
print (len(wfs.contents.keys()))
for layerID in list(wfs.contents.keys())[:6]:
layer = wfs[layerID]
print('Layer ID:', layerID)
print('Title:', layer.title)
#print('Boundaries:', layer.boundingBoxWGS84, '\n')
print ("BBOX:", layer.boundingBox, '\n')
5) Choose one layer id and show the title of the layer,
then get the vector data, save the feature data into json file
layer_id = 'oa:goaoninv'
meta = wfs.contents[layer_id]
print(meta.title)
# Get the vector data
data = wfs.getfeature(typename=[layer_id], bbox=(-179.5, -77.692, 179.20170000000002, 81.93299999999999), outputFormat='application/json')
# Save the to file
fn = 'wfsout.geojson'
with open(fn, 'w', encoding="utf-8") as fh:
fh.write(data.read())
6) Read the json file using geopandas
wfs_gdf = gpd.read_file(fn)
7) Choose the world map as backgroud map, show the world infomation
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.head(9)
Show the CRS inforamation for the world data set
print (world.crs)
Visualization of World Map¶
world.plot(figsize=(20, 20))
Visualization of World Map with Web Feature Service from OGC¶
Show the wfs layer information in the world map
wfs_gdf.plot(ax=world.plot(cmap='Set3', figsize=(20, 20)),
marker='o', color='purple', markersize=15);
Web Coverage Service¶
Web Coverage Service (WCS) is a standard that shares geospatial data as "Coverages" defined by OGC. The "Coverage" means that returns the data of any specified point in time-space domain and the form is easy to input the model. The WCS supports GetCapabilities, DescribeCoverage and GetCoverage operations. The three operations are required operations.
1) The GetCapabilities operation allows users to get service metadata from a WCS server and brief descriptions of the data collections;
2) The DescribeCoverage operation allows users to request full descriptions of one or more coverages;
3) The GetCoverage operation allows a client to request a coverage in common format.
This demo shows how to access WCS service.
1) import library, define WebCoverageService, define the url and version, print the wcs information
from owslib.wcs import WebCoverageService
import xarray as xr
import matplotlib.pyplot as plt
%matplotlib inline
wcs_url = 'http://geo.weather.gc.ca/geomet/?lang=en&service=WCS'
# Connect to service
wcs = WebCoverageService(wcs_url, version='2.0.1')
#Get the title of the service
print(wcs.identification.title)
# List the available first ten contents
sorted(list(wcs.contents.keys()))[:10]
2) Then get the layer information, and show the title, the BoundingBox and the formats.
layerid = 'OCEAN.GIOPS.3D_SALW_0000'
wcs_layer = wcs[layerid]
#Title
print('Layer title :', wcs_layer.title)
#bounding box
print('BoundingBox :', wcs_layer.boundingBoxWGS84)
# supported data formats - we'll use geotiff
print('Formats :', wcs_layer.supportedFormats)
#grid dimensions
print('Grid upper limits :', wcs_layer.grid.highlimits)
3) define the geographic projection, the bounding box, the resolution and format of the output, we can get the WCS result by using getCoverage.
format_wcs = 'image/netcdf'
bbox_wcs = wcs_layer.boundingboxes[0]['bbox'] # Get the entire domain
crs_wcs = wcs_layer.boundingboxes[0]['nativeSrs'] # Coordinate system
w = int(wcs_layer.grid.highlimits[0] )
h = int(wcs_layer.grid.highlimits[1])
print("Format:", format_wcs)
print("Bounding box:", bbox_wcs)
print("Projection:", crs_wcs)
print("Resolution: {} x {}".format(w, h))
output = wcs.getCoverage(identifier=[layerid, ], crs=crs_wcs, bbox=bbox_wcs, width=w, height=h, format=format_wcs)
4) Save the WCS as .nc file
wcsfn = layerid + '.nc'
with open(wcsfn, 'wb') as fh:
fh.write(output.read())
5) Read the .nc file and show the data in the map
wcsdt = xr.open_dataset(wcsfn)
print(wcsdt.data_vars)
plt.figure(figsize=(20,10))
wcsdt.Band1.plot()
plt.show()
Sensor Observation Service from OGC¶
The Sensor Observation Service (SOS) standard is applicable to use cases in which sensor data needs to be managed in an interoperable way. This standard defines a Web service interface which allows querying observations, sensor metadata, as well as representations of observed features. KVP binding and SOAP binding are specified in the SOS.
1) import libraries, list the amount of sensors
#SOS Demo
from owslib.sos import SensorObservationService
service = SensorObservationService('http://sensorweb.demo.52north.org/52n-sos-webapp/sos/kvp',version='2.0.0')
print (len(service.contents))
for content in sorted(service.contents):
print(content)
id = service.identification
print (id.title)
provider=service.provider
print (provider.name)
len(service.operations)
#get FOI
get_foi=service.get_operation_by_name('GetFeatureOfInterest')
try:
x = unicode('test')
for x in sorted(get_foi.parameters['featureOfInterest']['values']):
print(x.encode('utf8'))
except:
for x in sorted(get_foi.parameters['featureOfInterest']['values']):
print(x)
Web Map Tile Service¶
A Web Map Tile Service (WMTS) is a standard protocol for serving pre-rendered or run-time computed georeferenced map tiles over the Internet.
1) define the WebMapTileService object, define the url information
from owslib.wmts import WebMapTileService
wmts_url = "http://geodata.nationaalgeoregister.nl/tiles/service/wmts/ahn1?service=wmts&request=getcapabilities"
#wmts_url ="https://osmlab.github.io/wmts-osm/WMTSCapabilities.xml"
wmts = WebMapTileService(wmts_url)
2) show the WMTS information, inludeing the number of layes, version, title, the first layer name,
print (len(wmts.contents))
print (wmts.identification.type)
print (wmts.identification.version)
print (wmts.identification.title)
layer = sorted(list(wmts.contents))[0]
print (layer)
#print ((wmts.contents))
#print (wmts.contents[layer].styles)
print (wmts.contents[layer])
#print (list(wmts.contents))
3) Show operations in wmts, including getCapabilities, getTile and getfeatureinfo
for operation in wmts.operations:
print (operation.name)
4) fetech a tile, define the layername, tilematrixset, tilematrix, row, colum, and format
# Fetch a tile (using some defaults): EPSG:28992
tile = wmts.gettile(layer='opentopo',
tilematrixset='EPSG:28992', tilematrix='0',
row=0, column=0, format="image/jpeg")
5) save the tile as jpg
out = open('wmtsdemo.jpg', 'wb')
out.write(tile.read())
out.close()
6) plot the tile as a map
image_file = "wmtsdemo.jpg"
image = plt.imread(image_file)
fig, ax = plt.subplots()
ax.imshow(image)
ax.axis('off') # clear x-axis and y-axis
