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.

In [1]:
# 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))
74

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

In [2]:
%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

In [3]:
%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

In [4]:
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')
    
40
Layer ID: nvs:chl_avg_w3m_1d
Title: Chlorophyll-a Concentration Daily Average, Upper 3 meters
BBOX: (-135.0, 40.0, -122.0, 50.0, urn:ogc:def:crs:EPSG::4326) 

Layer ID: nvs:chl_avg_w3m_1m
Title: Chlorophyll-a Concentration Monthly Average, Upper 3 meters
BBOX: (-135.0, 40.0, -122.0, 50.0, urn:ogc:def:crs:EPSG::4326) 

Layer ID: nvs:chl_avg_w3m_7d
Title: Chlorophyll-a Concentration Weekly Average, Upper 3 meters
BBOX: (-135.0, 40.0, -122.0, 50.0, urn:ogc:def:crs:EPSG::4326) 

Layer ID: czo:wbd_huc10_EelRiver
Title: Eel River Basin Boundary
BBOX: (-124.3494266195082, 39.24061694638239, -122.66799631065977, 40.72629319824295, urn:ogc:def:crs:EPSG::4326) 

Layer ID: oa:goaoninv
Title: GOA-ON Assets
BBOX: (-179.5, -77.692, 179.2017, 81.93299999999999, urn:ogc:def:crs:EPSG::4326) 

Layer ID: nvs:nvs_stations
Title: NANOOS Complete Station Inventory from NVS (All stations on NVS)
BBOX: (-145.089, 40.294, -121.9556, 50.116, urn:ogc:def:crs:EPSG::4326) 

5) Choose one layer id and show the title of the layer,

then get the vector data, save the feature data into json file

In [5]:
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())
GOA-ON Assets

6) Read the json file using geopandas

In [6]:
wfs_gdf = gpd.read_file(fn)

7) Choose the world map as backgroud map, show the world infomation

In [7]:
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.head(9)
Out[7]:
pop_est continent name iso_a3 gdp_md_est geometry
0 920938 Oceania Fiji FJI 8374.0 MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1 53950935 Africa Tanzania TZA 150600.0 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2 603253 Africa W. Sahara ESH 906.5 POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3 35623680 North America Canada CAN 1674000.0 MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4 326625791 North America United States of America USA 18560000.0 MULTIPOLYGON (((-122.84000 49.00000, -120.0000...
5 18556698 Asia Kazakhstan KAZ 460700.0 POLYGON ((87.35997 49.21498, 86.59878 48.54918...
6 29748859 Asia Uzbekistan UZB 202300.0 POLYGON ((55.96819 41.30864, 55.92892 44.99586...
7 6909701 Oceania Papua New Guinea PNG 28020.0 MULTIPOLYGON (((141.00021 -2.60015, 142.73525 ...
8 260580739 Asia Indonesia IDN 3028000.0 MULTIPOLYGON (((141.00021 -2.60015, 141.01706 ...

Show the CRS inforamation for the world data set

In [8]:
print (world.crs)
{'init': 'epsg:4326'}

Visualization of World Map

In [9]:
world.plot(figsize=(20, 20))
Out[9]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f4da3250dd8>

Visualization of World Map with Web Feature Service from OGC

Show the wfs layer information in the world map

In [10]:
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

In [11]:
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]
Meteorological Service of Canada Geospatial Web Services 2.8.0
Out[11]:
['GDPS.CONV_KINDEX.PT3H',
 'GDPS.CONV_KINDEX.PT6H',
 'GDPS.CONV_ML-LCL-HGT.3h',
 'GDPS.CONV_ML-LCL-HGT.6h',
 'GDPS.CONV_ML-LI.400.3h',
 'GDPS.CONV_ML-LI.400.6h',
 'GDPS.CONV_ML-LI.500.3h',
 'GDPS.CONV_ML-LI.500.6h',
 'GDPS.CONV_ML-LI.600.3h',
 'GDPS.CONV_ML-LI.600.6h']

2) Then get the layer information, and show the title, the BoundingBox and the formats.

In [12]:
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)
Layer title : None
BoundingBox : None
Formats : ['image/tiff', 'image/x-aaigrid', 'image/netcdf', 'image/png', 'image/jpeg', 'image/png; mode=8bit']
Grid upper limits : ['1799', '849']

3) define the geographic projection, the bounding box, the resolution and format of the output, we can get the WCS result by using getCoverage.

In [13]:
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)
Format: image/netcdf
Bounding box: (-80.1, -180.0, 89.9, 180.0)
Projection: http://www.opengis.net/def/crs/EPSG/0/4326
Resolution: 1799 x 849

4) Save the WCS as .nc file

In [14]:
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

In [15]:
wcsdt = xr.open_dataset(wcsfn)

print(wcsdt.data_vars)

plt.figure(figsize=(20,10))

wcsdt.Band1.plot()

plt.show()
Data variables:
    crs      |S1 ...
    Band1    (lat, lon) float32 ...

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

In [16]:
#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))
33
In [17]:
for content in sorted(service.contents):
    print(content)
adxl345
apu
b
http://www.52north.org/test/offering/1
http://www.52north.org/test/offering/10
http://www.52north.org/test/offering/2
http://www.52north.org/test/offering/3
http://www.52north.org/test/offering/4
http://www.52north.org/test/offering/5
http://www.52north.org/test/offering/6
http://www.52north.org/test/offering/7
http://www.52north.org/test/offering/8
http://www.52north.org/test/offering/developer
idwniu
iok
k
miyo
myNewOffering
offering
sensor10_offering_AirTemperatur
sensor10_offering_Windspeed
sensor11_offering_AirTemperatur
sensor13_offering
sensor1_offering
sensor1_offering2
sensor1_offering_2
temp_senor_data
temp_sensor_data
test100
test55
testArd
testOffering2
www.edisoft.org/java/offering/1
In [18]:
id = service.identification
print (id.title)
provider=service.provider
print (provider.name)
len(service.operations)
52N SOS
52North
Out[18]:
16
In [19]:
#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)
adxl345
building
feature_1
http://www.52north.org/test/featureOfInterest/1
http://www.52north.org/test/featureOfInterest/2
http://www.52north.org/test/featureOfInterest/3
http://www.52north.org/test/featureOfInterest/4
http://www.52north.org/test/featureOfInterest/5
http://www.52north.org/test/featureOfInterest/6
http://www.52north.org/test/featureOfInterest/7
http://www.52north.org/test/featureOfInterest/8
http://www.52north.org/test/featureOfInterest/Heiden
http://www.52north.org/test/featureOfInterest/M√ľnster/FE101
http://www.52north.org/test/featureOfInterest/Portland
http://www.52north.org/test/featureOfInterest/TODO
http://www.52north.org/test/featureOfInterest/world
sensor1_id
testArd
www.edisoft.org/java/featureOfInterest/1
www.edisoft.org/java/featureOfInterest/2

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

In [20]:
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,

In [21]:
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))
44
OGC WMTS
1.0.0
Web Map Tile Service
2016_ortho25
Layer Name: 2016_ortho25 Title: 2016_ortho25

3) Show operations in wmts, including getCapabilities, getTile and getfeatureinfo

In [22]:
for operation in wmts.operations:
    print (operation.name)
GetCapabilities
GetTile
GetFeatureInfo

4) fetech a tile, define the layername, tilematrixset, tilematrix, row, colum, and format

In [23]:
# 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

In [24]:
out = open('wmtsdemo.jpg', 'wb')
out.write(tile.read())
out.close()

6) plot the tile as a map

In [25]:
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
Out[25]:
(-0.5, 255.5, 255.5, -0.5)