Using GSKY’s WMS in Python

Requesting map images through NCI’s GSKY Data Server

NCI’s GSKY Data Server supports the Open Geospatial Consortium (OGC) Web Map Service (WMS), which is a standard protocol for serving geospatial data as images (e.g., PNG).

The following libraries will need to be imported for this example:

[1]:
from owslib.wms import WebMapService
from PIL import Image
import matplotlib.pyplot as plt
%matplotlib inline

To start, we will need the base GSKY server URL:

[2]:
gsky_url = 'http://gsky.nci.org.au/ows/dea'

Using OWSLib, we can begin by inspecting the service metadata:

[3]:
wms = WebMapService(gsky_url, version='1.3.0')

To find out the available data layers that can be requested:

[4]:
for layer in list(wms.contents):
    print("Layer Name:", wms[layer].name)
    print("Title:", wms[layer].title, '\n')
Layer Name: hltc_high
Title: DEA High Tide Composite 25m v2.0

Layer Name: hltc_low
Title: DEA Low Tide Composite 25m v2.0

Layer Name: item_relative
Title: DEA Intertidal Extents Model Relative Layer 25m v2.0

Layer Name: item_stddev
Title: DEA Intertidal Extents Model Confidence Layer 25m v2.0

Layer Name: landsat5_geomedian
Title: DEA Landsat 5 terrain corrected surface reflectance geometric median

Layer Name: landsat5_nbar_16day
Title: 16-day DEA Landsat 5 surface reflectance

Layer Name: landsat5_nbar_daily
Title: Daily DEA Landsat 5 surface reflectance

Layer Name: landsat5_nbart_16day
Title: 16-day DEA Landsat 5 terrain corrected surface reflectance

Layer Name: landsat5_nbart_daily
Title: Daily DEA Landsat 5 terrain corrected surface reflectance

Layer Name: landsat7_geomedian
Title: DEA Landsat 7 terrain corrected surface reflectance geometric median

Layer Name: landsat7_nbar_16day
Title: 16-day DEA Landsat 7 surface reflectance

Layer Name: landsat7_nbar_daily
Title: Daily DEA Landsat 7 surface reflectance

Layer Name: landsat7_nbart_16day
Title: 16-day DEA Landsat 7 terrain corrected surface reflectance

Layer Name: landsat7_nbart_daily
Title: Daily DEA Landsat 7 terrain corrected surface reflectance

Layer Name: landsat8_geomedian
Title: DEA Landsat 8 terrain corrected surface reflectance geometric median

Layer Name: landsat8_nbar_16day
Title: 16-day DEA Landsat 8 surface reflectance

Layer Name: landsat8_nbar_daily
Title: Daily DEA Landsat 8 surface reflectance

Layer Name: landsat8_nbart_16day
Title: 16-day DEA Landsat 8 terrain corrected surface reflectance

Layer Name: landsat8_nbart_daily
Title: Daily DEA Landsat 8 terrain corrected surface reflectance

Layer Name: wofs
Title: DEA Water Observation Feature Layer

We can also view metadata that is available about a selected layer. For example, you can view the abstract associated with that data layer.

[5]:
layer = "landsat8_nbart_16day"
[6]:
print(wms[layer].abstract)
This product has been corrected to remove the influences of the atmosphere, the time of year, terrain shadow and satellite view angles using the methods described in Li et al. 2012 https://doi.org/10.1016/j.rse.2012.06.018. Landsat 8 Operational Land Imager (OLI) data is available from March 2013 and onwards. More detailed information about the terrain corrected surface reflectance product suite produced using Digital Earth Australia including CCBY4.0 is available at http://dx.doi.org/10.4225/25/5a7a76d2e129e. This service provides access to Landsat 8 OLI terrain corrected surface reflectance data. The image composites are made from images acquired within a 16 day period, and may include clouds.

Or view the CRS options, bounding box, and time positions available (these details will be needed to construct the GetMap request):

[7]:
print("CRS Options: ")
crs = sorted(wms[layer].crsOptions)
print('\t', crs, '\n')

print("Bounding Box: ")
bbox = wms[layer].boundingBox
print('\t', bbox, '\n')

print("Time Positions: ")
time = wms[layer].timepositions
print('\t', time[:10], '\n')

print("Styles: ")
styles = wms[layer].styles
print('\t', styles, '\n')
CRS Options:
         ['CRS:84', 'EPSG:3857', 'EPSG:4326']

Bounding Box:
         (-180.0, -90.0, 180.0, 90.0, 'CRS:84')

Time Positions:
         ['2013-03-19T00:00:00.000Z', '2013-04-04T00:00:00.000Z', '2013-04-20T00:00:00.000Z', '2013-05-06T00:00:00.000Z', '2013-05-22T00:00:00.000Z', '2013-06-07T00:00:00.000Z', '2013-06-23T00:00:00.000Z', '2013-07-09T00:00:00.000Z', '2013-07-25T00:00:00.000Z', '2013-08-10T00:00:00.000Z']

Styles:
         {'fc': {'title': 'False colour'}, 'tc': {'title': 'True colour'}}

Now let’s use the information above to construct and make the GetMap request.

We’ll need to define a bounding box for our request:

[8]:
subset_bbox = (147, -37, 148, -35)

OWSLib’s library can now be used to make the GetMap request:

[9]:
output = wms.getmap(layers=[layer],
                        Styles='tc',
                        srs=crs[2],
                        bbox=(subset_bbox[0], subset_bbox[1], subset_bbox[2], subset_bbox[3]),
                        size=(256, 256),
                        format='image/png',
                        time=time[2]
                        )

To view the above constructed URL:

[10]:
print(output.geturl())
http://gsky.nci.org.au/ows/dea?SERVICE=WMS&service=WMS&version=1.3.0&request=GetMap&layers=landsat8_nbart_16day&styles=&width=256&height=256&crs=EPSG%3A4326&bbox=-37%2C147%2C-35%2C148&format=image%2Fpng&transparent=FALSE&bgcolor=0xFFFFFF&exceptions=XML&time=2013-04-20T00%3A00%3A00.000Z&Styles=tc

Lastly, we need to write the GetMap result to a file:

[11]:
pngfile = './output/gsky_getMap.png'
with open(pngfile, 'wb') as out:
    out.write(output.read())

And if we’d like to confirm the result:

[12]:
im = Image.open(pngfile)
plt.figure(figsize=(10, 10))
plt.imshow(im)
[12]:
<matplotlib.image.AxesImage at 0x11ad6ccf8>
../_images/_notebook_Notebook_GSKY_WMS_25_1.png

For more information on the OGC WMS standard specifications and the Python OWSLib package: http://www.opengeospatial.org/standards/wms https://geopython.github.io/OWSLib/#wms