GSKY Click and Ship

The following libraries will need to be imported for the below example.

[1]:
from owslib.wms import WebMapService
from owslib.wcs import WebCoverageService
import ipywidgets as widgets
from IPython.display import display
from PIL import Image
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import numpy as np
%matplotlib inline

Preview available data

Let’s preview the available data first with WMS

[2]:
gsky_url = 'http://gsky.nci.org.au/ows/dea'
[3]:
wms = WebMapService(gsky_url, version='1.3.0')

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.
[7]:
print(wms[layer].styles)
{'fc': {'title': 'False colour'}, 'tc': {'title': 'True colour'}}

Let’s request a preview of the dataset through WMS

WMS function

[8]:
def get_map(layer, bbox, time, Styles):
    output = wms.getmap(layers=[layer],
                        srs=wms[layer].crsOptions[0],
                        bbox=bbox,
                        size=(256, 256),
                        format='image/png',
                        time=time,
                        Styles=style
                        )

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

    return pngfile

Make a request

[9]:
subset_bbox = (148, -37, 151, -36)
[10]:
print("Time Positions: ")
time = wms[layer].timepositions
print('\t', time[:10], '\n')
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']

[11]:
print(wms[layer].styles)
{'fc': {'title': 'False colour'}, 'tc': {'title': 'True colour'}}
[12]:
style = 'tc'
[13]:
pngfile = get_map(layer, subset_bbox, time[5], style)

And if we’d like to confirm the result:

[14]:
im = Image.open(pngfile)
plt.figure(figsize=(10, 10))
plt.imshow(im, extent=[subset_bbox[0], subset_bbox[1], subset_bbox[1], subset_bbox[0]])
[14]:
<matplotlib.image.AxesImage at 0x1202f4940>
../_images/_notebook_Notebook_GSKY_ClicknShip_24_1.png
[16]:
def get_coverage(layer, Styles, bbox, time):
    wcs = WebCoverageService(gsky_url, version='1.0.0')
    output = wcs.getCoverage(identifier=layer,
                             Styles=style,
                             time=time,
                             bbox=subset_bbox,
                             format='NetCDF',
                             crs='EPSG:4326',
                             width=256,
                             height=256)


    filename = './output/gsky_wcs.nc'
    with open(filename, 'wb') as f:
        f.write(output.read())

    return filename

Select/click region

[17]:
lon_range = widgets.IntRangeSlider(
    value=[135, 145],
    min=110,
    max=155,
    step=1,
    description='Longitude:',
    disabled=False,
    continuous_update=False,
    orientation='horizontal',
    readout=True,
    readout_format='d',
)

lat_range = widgets.IntRangeSlider(
    value=[-30, -25],
    min=-45,
    max=-10,
    step=1,
    description='Latitude:',
    disabled=False,
    continuous_update=False,
    orientation='horizontal',
    readout=True,
    readout_format='d',
)
[18]:
display(lon_range)
display(lat_range)
[19]:
print('latitude = {0}'.format(lat_range.value))
print('longitude = {0}'.format(lon_range.value))
latitude = (-30, -25)
longitude = (135, 145)
[20]:
subset = (lon_range.value[0], lat_range.value[0], lon_range.value[1], lon_range.value[1])
[21]:
style = 'fc'

And ship!

[22]:
wcs_output = get_coverage(layer, style, subset, time[5])
[23]:
with Dataset(wcs_output) as ds:

    band = ds['Band1']
    lon = ds['lon']
    lat = ds['lat']


    # Set figure size
    plt.figure(figsize=(12,12))

    # Plot image
    plt.imshow(band, extent=[lon[0], lon[-1], lat[-1], lat[0]])

    # Add figure title and labels
    # We can make use of the defined variable attributes to do this
    plt.title(layer.title(), fontsize=20)
    plt.xlabel(lon.long_name+' ('+lon.units+') ', fontsize=16)
    plt.ylabel(lat.long_name+' ('+lat.units+') ', fontsize=16)


    # Adjust tick mark size
    plt.tick_params(labelsize=16)
../_images/_notebook_Notebook_GSKY_ClicknShip_34_0.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