"""
Loading and processing Australian Bureau of Meteorology data.
License: The code in this notebook is licensed under the Apache License,
Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0). Digital Earth
Australia data is licensed under the Creative Commons by Attribution 4.0
license (https://creativecommons.org/licenses/by/4.0/).
Contact: If you need assistance, please post a question on the Open Data
Cube Discord chat (https://discord.com/invite/4hhBQVas5U) or on the GIS Stack
Exchange (https://gis.stackexchange.com/questions/ask?tags=open-data-cube)
using the `open-data-cube` tag (you can view previously asked questions
here: https://gis.stackexchange.com/questions/tagged/open-data-cube).
If you would like to report an issue with this script, you can file one
on GitHub (https://github.com/GeoscienceAustralia/dea-notebooks/issues/new).
Last modified: March 2021
"""
import datetime
import pytz
import ciso8601
from types import SimpleNamespace
import requests
import lxml
import lxml.etree
[docs]
def get_stations(time=None,
observation='http://bom.gov.au/waterdata/services/parameters/Water Course Discharge',
url='http://www.bom.gov.au/waterdata/services'):
""" Get list of stations
:param time: tuple of datetime.datetime objects, or None to query from 1980-1-1 to Now
Returns
========
List of stations:
.name -- string, human readable station name
.pos -- Coordinate of the station or None
.url -- service url identifier
"""
data = tpl_get_stations.format(observation=observation,
**_fmt_time(time))
data = data.replace('\n', '')
rr = requests.post(url, data=data)
return _parse_station_data(rr.text)
[docs]
def get_station_data(station,
time=None,
observation='http://bom.gov.au/waterdata/services/parameters/Water Course Discharge',
url='http://www.bom.gov.au/waterdata/services'):
"""
Query Gauge Data.
:param station: One of the stations see get_stations
:param time: tuple of datetime.datetime objects, or None to query from 1980-1-1 to Now
Returns
========
Pandas dataframe with Timestamp(index), Value columns
"""
data = tpl_get_obs.format(station=station.url,
observation=observation,
**_fmt_time(time))
data = data.replace('\n', '')
rr = requests.post(url, data=data)
return _parse_get_data(rr.text)
[docs]
def mk_station_selector(on_select,
stations=None,
dst_map=None,
**kw):
"""
Add stations to the map and register on_click event.
:param on_select: Will be called when user selects station on the map `on_select(station)`
:param stations: List of stations as returned from get_stations
:param dst_map: Map to add stations markers to
Any other arguments are passed on to Map(..) constructor.
Returns
=======
(map, marker_cluster)
Passes through map=dst_map if not None, or returns newly constructed Map object.
"""
import ipyleaflet as L
if stations is None:
stations = get_stations()
stations = [st for st in stations if st.pos is not None]
pos2st = {st.pos: st for st in stations}
def on_click(event='', type='', coordinates=None):
pos = tuple(coordinates)
st = pos2st.get(pos)
if st is None:
# should probably log warning here
print("Can't map click to station")
return
on_select(st)
markers = [L.Marker(location=st.pos,
draggable=False,
title=st.name)
for st in stations]
cluster = L.MarkerCluster(markers=markers)
if dst_map is None:
dst_map = L.Map(**kw)
dst_map.add_layer(cluster)
cluster.on_click(on_click)
return dst_map, cluster
[docs]
def ui_select_station(stations,
zoom=3,
center=(-24, 138),
**kw):
"""
Create an interactive map for selecting river gauging stations.
"""
import ipywidgets as W
from IPython.display import display
import matplotlib.pyplot as plt
import ipyleaflet as L
from odc.ui import ui_poll
dbg_display = W.Output()
fig_display = W.Output()
btn_done = W.Button(description='Done')
scroll_wheel_zoom = kw.pop('scroll_wheel_zoom', True)
map_widget = L.Map(zoom=zoom,
center=center,
scroll_wheel_zoom=scroll_wheel_zoom,
**kw)
state = SimpleNamespace(pos=None,
gauge_data=None,
finished=False,
station=None)
plt_interactive_state = plt.isinteractive()
plt.interactive(False)
with fig_display:
fig, ax = plt.subplots(1, figsize=(14,4))
ax.set_visible(False)
display(fig)
def _on_select(station):
if state.finished:
print('Please re-run the cell')
return
state.station = station
state.pos = station.pos
state.gauge_data = None
print('Fetching data for: {}'.format(station.name))
try:
xx = get_station_data(station).dropna()
except Exception:
print('Failed to read data')
return
print('Got {} observations'.format(xx.shape[0]))
state.gauge_data = xx
with fig_display:
ax.clear()
ax.set_visible(True)
xx.plot(ax=ax)
ax.set_xlabel("Date")
ax.set_ylabel("Cubic meters per second")
ax.legend([station.name])
fig_display.clear_output(wait=True)
with fig_display:
display(fig)
def on_select(station):
with dbg_display:
_on_select(station)
def on_done(btn):
if state.finished:
with dbg_display:
print('Please re-run the cell')
return
state.finished = True
n_obs = 0 if state.gauge_data is None else state.gauge_data.shape[0]
with dbg_display:
print('''Finished
Station: {}
Number of Observations: {}'''.format(state.station.name, n_obs))
def on_poll():
with dbg_display:
if state.finished:
return state.gauge_data, state.station
return None
mk_station_selector(on_select,
stations=stations,
dst_map=map_widget)
## UI:
##
## MMMMMMMMMMMMM BBBBB
## MMMMMMMMMMMMM .....
## MMMMMMMMMMMMM .....
## MMMMMMMMMMMMM .....
## MMMMMMMMMMMMM .....
## FFFFFFFFFFFFFFFFFFF
## FFFFFFFFFFFFFFFFFFF
# M - Map F - Figure
# B - Button . - Debug output
btn_done.on_click(on_done)
r_panel = W.VBox([btn_done, dbg_display],
layout=W.Layout(width='30%'))
ui = W.VBox([W.HBox([map_widget, r_panel]),
fig_display])
display(ui)
result = ui_poll(on_poll, 1/20) # this will block until done is pressed
#restore interactive state
fig_display.clear_output(wait=True)
with fig_display:
plt.interactive(plt_interactive_state)
plt.show()
return result
def _fmt_time(time=None):
if time is None:
time = (datetime.datetime(1980, 1, 1),
datetime.datetime.now())
t_start, t_end = (t.isoformat() for t in time)
return dict(t_start=t_start, t_end=t_end)
def _parse_float(x):
if x is None:
return float('nan')
try:
return float(x)
except ValueError:
return float('nan')
def _parse_time(x):
t = ciso8601.parse_datetime(x).astimezone(pytz.utc)
return t.replace(tzinfo=None)
def _parse_get_data(text, raw=False):
root = lxml.etree.fromstring(text)
data = [[e.text for e in root.findall('.//{http://www.opengis.net/waterml/2.0}' + t)]
for t in ['time', 'value']]
dd = [(_parse_time(t),
_parse_float(v))
for t, v in zip(*data)]
if raw:
return dd
import pandas as pd
return pd.DataFrame(dd, columns=('Timestamp', 'Value')).set_index('Timestamp')
def _parse_station_data(text):
def parse_pos(pos):
if pos is None:
return None
return tuple(_parse_float(x)
for x in pos.split(' '))
root = lxml.etree.fromstring(text)
data = [[e.text for e in root.findall('.//{http://www.opengis.net/gml/3.2}' + t)]
for t in ['name', 'identifier', 'pos']]
return [SimpleNamespace(name=name, url=url, pos=parse_pos(pos))
for name, url, pos in zip(*data)]
# observation = 'http://bom.gov.au/waterdata/services/parameters/Water Course Discharge'
#
tpl_get_stations = '''
<soap12:Envelope xmlns:soap12="http://www.w3.org/2003/05/soap-envelope"
xmlns:sos="http://www.opengis.net/sos/2.0"
xmlns:wsa="http://www.w3.org/2005/08/addressing"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:ows="http://www.opengis.net/ows/1.1"
xmlns:fes="http://www.opengis.net/fes/2.0"
xmlns:gml="http://www.opengis.net/gml/3.2"
xmlns:swes="http://www.opengis.net/swes/2.0"
xsi:schemaLocation="http://www.w3.org/2003/05/soap-envelope http://www.w3.org/2003/05/soap-envelope/soap-envelope.xsd http://www.opengis.net/sos/2.0 http://schemas.opengis.net/sos/2.0/sos.xsd"
>
<soap12:Header>
<wsa:To>http://www.ogc.org/SOS</wsa:To>
<wsa:Action>http://www.opengis.net/def/serviceOperation/sos/foiRetrieval/2.0/GetFeatureOfInterest</wsa:Action>
<wsa:ReplyTo>
<wsa:Address>http://www.w3.org/2005/08/addressing/anonymous</wsa:Address>
</wsa:ReplyTo>
<wsa:MessageID>0</wsa:MessageID>
</soap12:Header>
<soap12:Body>
<sos:GetFeatureOfInterest service="SOS" version="2.0.0">
<sos:observedProperty>{observation}</sos:observedProperty>
<sos:temporalFilter>
<fes:During>
<fes:ValueReference>om:phenomenonTime</fes:ValueReference>
<gml:TimePeriod gml:id="tp1">
<gml:beginPosition>{t_start}</gml:beginPosition>
<gml:endPosition>{t_end}</gml:endPosition>
</gml:TimePeriod>
</fes:During>
</sos:temporalFilter>
</sos:GetFeatureOfInterest>
</soap12:Body>
</soap12:Envelope>
'''
# {station}, {observation}, {t_start}, {t_end}
tpl_get_obs = '''
<soap12:Envelope xmlns:soap12="http://www.w3.org/2003/05/soap-envelope"
xmlns:sos="http://www.opengis.net/sos/2.0"
xmlns:wsa="http://www.w3.org/2005/08/addressing"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:ows="http://www.opengis.net/ows/1.1"
xmlns:fes="http://www.opengis.net/fes/2.0"
xmlns:gml="http://www.opengis.net/gml/3.2"
xmlns:swes="http://www.opengis.net/swes/2.0"
xsi:schemaLocation="http://www.w3.org/2003/05/soap-envelope http://www.w3.org/2003/05/soap-envelope/soap-envelope.xsd http://www.opengis.net/sos/2.0 http://schemas.opengis.net/sos/2.0/sos.xsd"
>
<soap12:Header>
<wsa:To>http://www.ogc.org/SOS</wsa:To>
<wsa:Action>http://www.opengis.net/def/serviceOperation/sos/core/2.0/GetObservation</wsa:Action>
<wsa:ReplyTo>
<wsa:Address>http://www.w3.org/2005/08/addressing/anonymous</wsa:Address>
</wsa:ReplyTo>
<wsa:MessageID>0</wsa:MessageID>
</soap12:Header>
<soap12:Body>
<sos:GetObservation service="SOS" version="2.0.0">
<sos:procedure>http://bom.gov.au/waterdata/services/tstypes/Pat4_C_B_1_DailyMean</sos:procedure>
<sos:observedProperty>{observation}</sos:observedProperty>
<sos:featureOfInterest>{station}</sos:featureOfInterest>
<sos:temporalFilter>
<fes:During>
<fes:ValueReference>om:phenomenonTime</fes:ValueReference>
<gml:TimePeriod gml:id="tp1">
<gml:beginPosition>{t_start}</gml:beginPosition>
<gml:endPosition>{t_end}</gml:endPosition>
</gml:TimePeriod>
</fes:During>
</sos:temporalFilter>
</sos:GetObservation>
</soap12:Body>
</soap12:Envelope>
'''
