Requesting National Geophysical Compilations data subsets through NCI’s GSKY Data Server¶
NCI’s GSKY Data Server supports the Open Geospatial Consortium (OGC) Web Coverage Service (WCS), which is a standard protocol for serving geospatial data in common formats such as NetCDF and GeoTIFF.
Constructing WCS Requests
GSKY’s Web Coverage Service (WCS) allows users to request data or subsets of data in either NetCDF3 or GeoTIFF format. The request is made by constructing a GetCoverage URL, which is then used within a web browser to communicate to the GSKY Data Server.
For example, the GetCoverage request for the aster layer takes the following form:
http://gsky.nci.org.au/ows/national_geophysical_compilations? service = WCS &version = 1.0.0 &request = GetCoverage &coverage = value &format = value &bbox =value &time =value &crs =value
GetCoverage parameters:
| Parameter | Required/Optional | Input |
|---|---|---|
| service | Required | WCS |
| version | Required | 1.0.0 |
| request | Required | GetCoverage |
| coverage | Required | <variable> |
| format | Required | GeoTIFF, GeoTIFF_Float, NetCDF3 |
bbox* |
Required/Optional | <xmin,ymin,xmax,ymax> |
time* |
Required/Optional | <time_value> |
| srs, or crs | Optional | <srs_value> or <crs_value> |
*For large files and/or files with a time dimension, these might be necessary. If bbox is not defined the entire spatial domain will be returned (if server limits allow) and if time is not specified, either the first or sometimes last timestep is returned.
WCS GetCapabilities and DescribeCoverage
In order to contruct the GetCoverage URL, a GetCapabilities request can be made to the server. This requests returns a xml describing the available WCS parameters (metadata, services, and data) made available by NCI’s GSKY server. Additional metadata information can also be requested about a specific coverage layer by making a DescribeCoverage request.
GetCapabilities example:
DescribeCoverage example: > http://gsky.nci.org.au/ows/national_geophysical_compilations?service=WCS&version=1.0.0&coverage=magmap2019_grid_tmi_awags_mag_2019&request=DescribeCoverage
GetCoverage Request Using the information returned from the GetCapabilities and DescribeCoverage requests, a GetMap URL can be constructed and then entered into the address bar of any web browser.
Example GetCoverage (NetCDF format):
Example GetCoverage (GeoTIFF format):
Using GSKY’s WCS in Python¶
Many tools are available to perform the above steps less manually. In particular, Python’s OWSLib library can be used with GSKY’s WCS.
The following libraries will need to be imported for the below example
[1]:
from owslib.wcs import WebCoverageService
from PIL import Image
%matplotlib inline
To start, we will need the base GSKY server URL:
[2]:
gsky_url = 'http://gsky.nci.org.au/ows/national_geophysical_compilations?service=WCS&version=1.0.0&request=GetCapabilities'
Now using OWSLib, you can begin by inspecting the service metadata:
[3]:
wcs = WebCoverageService(gsky_url)
Find out the available data layers that can be requested:
[4]:
for layer in list(wcs.contents):
print("Layer Name:", layer)
print("Title:", wcs[layer].title, '\n')
Layer Name: gravmap2016_grid_grv_cscba
Title: Australia gravity grid 2016 (complete spherical cap Bouguer anomaly)
Layer Name: gravmap2016_grid_grv_ir
Title: Australia gravity grid 2016 (isostatic residual anomaly)
Layer Name: gravmap2016_grid_grv_scba
Title: Australia gravity grid 2016 (spherical cap Bouguer anomaly)
Layer Name: magmap2019_grid_tmi_1vd_awags_mag_2019
Title: Total Magnetic Intensity Grid of Australia 2019 - First Vertical Derivative (1VD)
Layer Name: magmap2019_grid_tmi_awags_mag_2019
Title: Total Magnetic Intensity (TMI) Grid of Australia 2019 - seventh edition - 80 m cell size
Layer Name: magmap2019_grid_tmi_cellsize40m_awags_mag_2019
Title: Total Magnetic Intensity (TMI) Grid of Australia 2019 - seventh edition - 40 m cell size
Layer Name: magmap2019_grid_tmi_rtp_awags_mag_2019
Title: Total Magnetic Intensity (TMI) Grid of Australia with Variable Reduction to Pole (VRTP) 2019 - seventh edition
Layer Name: radmap2019_grid_dose_terr_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 unfiltered terrestrial dose rate
Layer Name: radmap2019_grid_dose_terr_filtered_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 filtered terrestrial dose rate
Layer Name: radmap2019_grid_k_conc_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 unfiltered pct potassium
Layer Name: radmap2019_grid_k_conc_filtered_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 filtered pct potassium grid
Layer Name: radmap2019_grid_th_conc_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 unfiltered ppm thorium
Layer Name: radmap2019_grid_th_conc_filtered_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 filtered ppm thorium
Layer Name: radmap2019_grid_thk_ratio_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 ratio thorium over potassium
Layer Name: radmap2019_grid_u2th_ratio_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 ratio uranium squared over thorium
Layer Name: radmap2019_grid_u_conc_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 unfiltered ppm uranium
Layer Name: radmap2019_grid_u_conc_filtered_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 filtered ppm uranium
Layer Name: radmap2019_grid_uk_ratio_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 ratio uranium over potassium
Layer Name: radmap2019_grid_uth_ratio_awags_rad_2019
Title: Radiometric Grid of Australia (Radmap) v4 2019 ratio uranium over thorium
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 = "magmap2019_grid_tmi_1vd_awags_mag_2019"
[6]:
print(wcs[layer].abstract)
Total magnetic intensity (TMI) data measures variations in the intensity of the Earth magnetic field caused by the contrasting content of rock-forming minerals in the Earth crust. Magnetic anomalies can be either positive (field stronger than normal) or negative (field weaker) depending on the susceptibility of the rock.The first vertical derivative (1VD) grid is derived from the 2019 Total magnetic Intensity (TMI) grid of Australia which has a grid cell size of ~3 seconds of arc (approximately 80 m). As the vertical derivative filter is essentially a high-pass filter, longer wavelengths are suppressed, and shorter wavelengths emphasized. The magnetic units of the data are in nT per metre.
NCI Data Catalogue: https://dx.doi.org/10.25914/5f75640aae140
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 = wcs[layer].supportedCRS
print('\t', crs, '\n')
print("Bounding Box: ")
bbox = wcs[layer].boundingBoxWGS84
print('\t', bbox, '\n')
print("Time Positions: ")
time = wcs[layer].timepositions
print('\t', time[:10], '\n')
CRS Options:
[urn:ogc:def:crs:EPSG::4326]
Bounding Box:
(-180.0, -90.0, 180.0, 90.0)
Time Positions:
['2019-08-29T00:00:00.000Z']
Now let’s use the information above to construct and make GetCoverage requests The below sections will demonstrate both a request in GeoTIFF and NetCDF formats.
We’ll need to define a bounding box for our request:
[8]:
subset_bbox = (112, -44, 154, -10)
OWSLib’s library can now be used to make the GetCoverage request:
[9]:
output = wcs.getCoverage(identifier=layer,
time=[wcs[layer].timepositions[0]],
bbox=subset_bbox,format='GeoTIFF',
crs='EPSG:4326', width=256, height=256)
To view the above constructed URL:
[10]:
print(output.geturl())
http://gsky.nci.org.au/ows/national_geophysical_compilations?version=1.0.0&request=GetCoverage&service=WCS&Coverage=magmap2019_grid_tmi_1vd_awags_mag_2019&BBox=112%2C-44%2C154%2C-10&time=2019-08-29T00%3A00%3A00.000Z&crs=EPSG%3A4326&format=GeoTIFF&width=256&height=256
Write the result to a file:
[11]:
filename = './gsky_wcs.tiff'
with open(filename, 'wb') as f:
f.write(output.read())
And if we’d like to confirm the result, we can open and view the GeoTIFF with the Python GDAL library for example:
[12]:
import gdal
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
[13]:
ds = gdal.Open(filename)
band = ds.GetRasterBand(1).ReadAsArray()
fig = plt.figure(figsize=(15,15))
plt.imshow(band, cmap='binary')
plt.clim(-0.03,0.1)
To request a coverage returned as in the NetCDF format, a similar GetCoverage request is constructed with the format parameter now specifying the NetCDF option.
[14]:
layer = "gravmap2016_grid_grv_ir"
[15]:
output = wcs.getCoverage(identifier=layer,
time=[wcs[layer].timepositions[0]],
bbox=subset_bbox,format='NetCDF',
crs='EPSG:4326', width=256, height=256)
[16]:
print(output.geturl())
http://gsky.nci.org.au/ows/national_geophysical_compilations?version=1.0.0&request=GetCoverage&service=WCS&Coverage=gravmap2016_grid_grv_ir&BBox=112%2C-44%2C154%2C-10&time=2019-08-29T00%3A00%3A00.000Z&crs=EPSG%3A4326&format=NetCDF&width=256&height=256
Again, write the output to a file to save:
[17]:
filename = './gsky_wcs.nc'
with open(filename, 'wb') as f:
f.write(output.read())
To confirm or inspect the contents of the NetCDF file, libraries such as NetCDF4 Python or GDAL can be used.
[18]:
from netCDF4 import Dataset
[19]:
with Dataset(filename) as ds:
band = ds['Band1'][::-1]
fig = plt.figure(figsize=(10,10))
plt.imshow(band, cmap='afmhot')
For more information on the OGC WCS standard specifications and the Python OWSLib package: http://www.opengeospatial.org/standards/wcs