Request GSKY’s WMS from Python (bushfire example)

In this notebook,

  • request images through GSKY WMS endpoints
  • create a GIF using those images

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).

In this example, we extract data from the “Multi-sensor (Landsat and Sentinel 2) surface reflectance false colour (Beta)” (DEA products) to view a bushfire in NSW, Australia over two days, 12-13 of September, 2019.

You can view this example via Terria Map.

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

[11]:
from owslib.wms import WebMapService
from PIL import Image, ImageDraw, ImageFilter
import numpy
import matplotlib.pyplot as plt
%matplotlib inline

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

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

Using OWSLib, we can begin by inspecting the service metadata to find out the available data layers:

[13]:
wms = WebMapService(gsky_url, version='1.3.0')
for layer in list(wms.contents):
    print("Layer Name:", wms[layer].name)
    print("Title:", wms[layer].title, '\n')
Layer Name: blend_sentinel2_landsat_nbart_daily_false_colour
Title: Multi-sensor (Landsat and Sentinel 2) surface reflectance false colour (Beta)

Layer Name: blend_sentinel2_landsat_nbart_daily_true_colour
Title: Multi-sensor (Landsat and Sentinel 2) surface reflectance true colour (Beta)

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: sentinel2_nbart_daily
Title: Sentinel 2 Analysis Ready Data

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.

[14]:
layer = "blend_sentinel2_landsat_nbart_daily_false_colour"

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

[55]:
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[12081:12083], '\n')
CRS Options:
         ['CRS:84', 'EPSG:3857', 'EPSG:4326']

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

Time Positions:
         ['2019-09-12T00:00:00.000Z']

Now, we can construct the bounding box and time selection for our bushfire example in NSW Australia, over two days, 12-13 of September, 2019.

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

[84]:
subset_bbox = (152.3, -30.1, 152.45, -30.02)
#subset_bbox = (152, -31.2, 152.3, -30.5) You can modify this bounding box to look at different area.

OWSLib’s library can now be used to make the GetMap request and print out the images:

[87]:
# 2019-09-12 - 09-13
dates = time[12081:12083]

images = []

for date in dates:
    print(date)
    img = wms.getmap(layers=[layer],
                    srs=crs[2],
                    bbox=(subset_bbox[0], subset_bbox[1], subset_bbox[2], subset_bbox[3]),
                    size=(512, 256),
                    format='image/png',
                    time=date
                    )
    im = Image.open(img)
    images.append(im)
    im.save('gsky_getMap%.10s.png' % date)

    im = Image.open(img)
    fig = plt.figure
    imgplot = plt.imshow(im,extent=[subset_bbox[0], subset_bbox[2], subset_bbox[1], subset_bbox[3]])


images[0].save('bushfire_NSW.gif',save_all=True, append_images=images[1:],optimize=False,duration=2,loop=0)
2019-09-12T00:00:00.000Z
2019-09-13T00:00:00.000Z
../_images/_notebook_request_GSKY_WMS_sentinel2_bushfire_NSW_Sep2019_15_1.png

Another way to make GIF is using https://giphy.com/create/gifmaker, where you can upload all .png files, and create a fancier GIF. For example, you can add annotation.


For more information on the OGC WMS standard specifications and the Python OWSLib package: