# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Exploring Observation Statistics From TILDE, Application to DART Data

# ## &nbsp;Table of contents
# ### 1. Introduction
# ### 2. Building a Basic TILDE Query and Extracting the Key Statistics
# ### 3. Building a Query using startdate and enddate parameters
# ### 4. Building a Query using parameters from the last 7 days
# ### 5. Building a Query using Percentiles
# ### 6. Comparing query results from several stations

# ## &nbsp;1. Introduction

# In this tutorial we will be learning how to use Python to access the TILDE API `Statistics` endpoint. To highlight the different functionality and the statistics available we will be using the DART (Deep-ocean Assessment and Reporting of Tsunamis) dataset. GeoNet has an API (Application Programming Interface) to access its DART data. You do not need to know anything about APIs to use this tutorial. If you would like more info see https://tilde.geonet.org.nz/api-docs/.
#
# This tutorial assumes you have basic knowledge of Python.

# ###### About GeoNet DART data

# GeoNet uses the 12 DART Tsunameters deployed offshore New Zealand and around the Southwestern Pacific Ocean to monitor ocean height. When a change has been detected of a certain magnitude, a buoy will "trigger" and go into a heightened detection mode. The DARTs have two operational reporting modes; standard and event. When in standard reporting mode, the BPR (bottom pressure recorder) and buoy system send four six-hour bundles of 15 minute water height values. When in event reporting mode, BPR data are sampled at 15 second intervals and are sent more regularly. The buoy surface location (latitude and longitude) will also be sent daily. <br>
# TILDE provides access to the 15 minutes and 15 second sampled data.
#
# For more DART information see the GeoNet page: https://www.geonet.org.nz/tsunami/dart

# ## &nbsp;2. Building a Basic TILDE Query and Extracting the Key Statistics

# ###### Import required modules and set the source URL

import requests
import json
import pandas as pd
import matplotlib.pyplot as plt

source = 'https://tilde.geonet.org.nz'

# ### A basic query with the minimum required parameters
#
# The query is `https://tilde.geonet.org.nz/v1/stats/` with the source URL being stated above `https://tilde.geonet.org.nz`.
#
# The minimum required parameters are:
# - domain = `dart`
# - key = `NZE/water-height/40/WTZ`, this specifies a station `NZE`, a data name `water-height`, a sensor code `40`, and a data type `WTZ`. A data type of `WTZ` is 15 minute sampled data. A data type of `UTZ` is 15 second sampled (triggered) data.
# With this information we can now use a more specific URL.
#
# Optional parameters are:
# - startdate = `latest` or `2021-01-01`
# - enddate = `30d` or `2021-01-31`

url = source+'/v1/stats/dart/NZE/water-height/40/WTZ'

r = requests.get(url)
print (r)

# We use `requests.get` to retrieve the data from the URL. The response status code says whether we were successful in getting the data requested and if we were unsuccessful, why not :
# <ul>
# <li>200 -- everything went okay, and the result has been returned (if any)
# <li>301 -- the server is redirecting you to a different endpoint. This can happen when a company switches domain names, or an endpoint name is changed.
# <li>400 -- the server thinks you made a bad request. This can happen when you don't send along the right data, among other things.
# <li>404 -- the resource you tried to access wasn't found on the server.
# </ul>

# ###### Store the response as a nested dictionary and see the output

# A nested dictionary in python is a dictionary within a dictionary. It's a collection of dictionaries into one single directory. In TILDE, each domain has different entries which sit within that specific domain name, e.g. for example, the DART entries sit within the DART domain. By creating a nested dictionary we can then specify the domain and then split it by domain name. `summary` is our nested dictionary and when printed we can see that there are a lot of key:value pairs. 

summary = r.json()
type(summary)

summary

# ##### We can then specify the series
# <br>
# Extract the key statistics and display them in a more readable format. To do this, we select the 'series' key of the dictionary and assign it to a new dictionary called `my_series`. Printing the output of `my_series`, can also be compared to `summary` above to see what happens when we select this key.

my_series = summary['series']
my_series

# We can also add some other key statistics from our dictionary `summary`, by using the same method that we did for 'series'. To specify the statistics required, we need to make sure we are selected the correct values. For example, within `summary`, the dictionary `first` has multiple keys ('value', 'err', 'qc', and 'ts'). To specify the first observed time we would need `summary['first']['ts']`. For other keys within `summary` that don't have sub-keys (e.g. count, mean, median), we can just use `summary['count']`. 

my_series['first obs'] = summary['first']['ts']
my_series['last obs'] = summary['last']['ts']
my_series['maximum obs'] = summary['maximum']['value']
my_series['minimum obs'] = summary['minimum']['value']
my_series['number obs'] = summary['count']
my_series['mean obs'] = summary['mean']
my_series['median obs'] = summary['median']
my_series['population stdev'] = summary['stddevPopulation']

# We can now print `my_series` in a nice, easily readable table by using tabulate (https://pypi.org/project/tabulate/), which allows to print items to a table using (my_series.item()). The type of table format has many different options, such as 'simple', 'grid', 'fancy_grid', and 'presto' available. 

# If `tabulate` is not installed (if using an older version of the `environment.yml` file), either run `conda env update -f environment.yml` after pulling the most recent environment.yml file from the data-tutorials github repository or run `conda install -c conda-forge tabulate` in a terminal within the GeoNet anaconda environment.

# +
from tabulate import tabulate

headers = ['Parameter', 'Value']
data = [(k,v) for k,v in my_series.items()]
print (tabulate(data, headers=headers, tablefmt='github'))
# -

# ## &nbsp;3. Building a Query using startdate and enddate parameters
#
# The minimum required parameters are:
# - domain = `dart`
# - key = `NZE/water-height/40/WTZ`, this specifies a station `NZE`, a data name `water-height`, a sensor code `40`, and a data type `WTZ`. A data type of `WTZ` is 15 minute sampled data. A data type of `UTZ` is 15 second sampled data.
#
# We will now be adding a start date and end date of:
# - startdate = `2021-01-01`
# - enddate = `2021-01-31`
#
# For this, we will begin by changing our URL and then generating a nested dictionary called `summary`, as shown previously. Printing the summary will show a very similar summary to above, but this time the first and last datapoints we be what we have specified in the URL. Feel free to test this by uncommenting `summary` to see what the output is. 

# +
url = source+'/v1/stats/dart/NZE/water-height/40/WTZ/2021-01-01/2021-01-31'

r = requests.get(url)
summary = r.json()
#summary
# -

# We will now use a similar method as we did before and select the `series` component of the dictionary as well as adding some other key statistics as we did above. 

my_series2 = summary['series']
my_series2

my_series2['first obs'] = summary['first']['ts']
my_series2['last obs'] = summary['last']['ts']
my_series2['maximum obs'] = summary['maximum']['value']
my_series2['minimum obs'] = summary['minimum']['value']
my_series2['number obs'] = summary['count']
my_series2['mean obs'] = summary['mean']
my_series2['median obs'] = summary['median']
my_series2['population stdev'] = summary['stddevPopulation']

# As before we will now use tabulate and print this into a table. In this table we can see that the `first obs` and the `last obs` are using the start and end dates that we specified in the URL.

headers = ['Parameter', 'Value']
data = [(k,v) for k,v in my_series2.items()]
print (tabulate(data, headers=headers, tablefmt='github'))

# ## &nbsp;4. Building a Query using parameters from the last 7 days
#
# The minimum required parameters are:
# - domain = `dart`
# - key = `NZE/water-height/40/WTZ`, this specifies a station `NZE`, a data name `water-height`, a sensor code `40`, and a data type `WTZ`. A data type of `WTZ` is 15 minute sampled data. A data type of `UTZ` is 15 second sampled data.
#
# We will now be adding:
# - `latest`
# - period = `7d`
#
# The latest `period` can be 30d (latest 30 days), 7d, 2d, 1d, or 6h (latest 6 hours). Other periods will produce an error.
#
# When calculating the start of the query, the current time is taken into account.
#
# We will begin again by changing our URL and then generating a nested dictionary called `summary`, as shown previously. Printing the summary will show a very similar summary to above, but this time we will be seeing the first and last data points for the data period that has been selected, in this case the last 7 days from the current time. Feel free to test this by uncommenting `summary` to see what the output is.  

# +
url = source+'/v1/stats/dart/NZE/water-height/40/WTZ/latest/7d'

r = requests.get(url)
summary = r.json()
# summary
# -

# We will now use a similar method as we have done before and select the `series` component of the dictionary as well as adding some other key statistics. 

#selecting the series component of the `summary` dictionary
my_series3 = summary['series']

#selecting key statistics
my_series3['first obs'] = summary['first']['ts']
my_series3['last obs'] = summary['last']['ts']
my_series3['maximum obs'] = summary['maximum']['value']
my_series3['minimum obs'] = summary['minimum']['value']
my_series3['number obs'] = summary['count']
my_series3['mean obs'] = summary['mean']
my_series3['median obs'] = summary['median']
my_series3['population stdev'] = summary['stddevPopulation']

# As before we will now use tabulate and print this into a table. In this table we can see that the `first obs` and the `last obs` are using the data period that we specified in the URL.

headers = ['Parameter', 'Value']
data = [(k,v) for k,v in my_series3.items()]
print (tabulate(data, headers=headers, tablefmt='github'))

# ## &nbsp;5. Building a Query using Percentiles
#
# The minimum required parameters are:
# - domain = `dart`
# - key = `NZE/water-height/40/WTZ`, this specifies a station `NZE`, a data name `water-height`, a sensor code `40`, and a data type `WTZ`. A data type of `WTZ` is 15 minute sampled data. A data type of `UTZ` is 15 second sampled data.
#
# We will now be adding:
# - `?percentiles=10,75`, 10th and 75th percentiles
#
# We will begin again by changing our URL and then generating a nested dictionary called `summary`, as shown previously. Printing the summary will show a very similar summary to above, but this time we will be seeing the 10th and 75th percentiles Feel free to test this by uncommenting `summary` to see what the output is. In our previous `summary` dictionaries we had `percentile=None`, now we are seeing a specific result. To change to different percentiles, e.g. 25 and 70, just change the URL `?percentiles=25,70`.

# +
url = source+'/v1/stats/dart/NZE/water-height/40/WTZ/?percentiles=10,75'

r = requests.get(url)
summary = r.json()
#summary
# -

# #### Box-and-whisker plot.
# If we specified a low and a high percentile, we can summarise the distribution of observations, with a box-and-whisker style plot (https://en.wikipedia.org/wiki/Box_plot). We begin by specifying the median, minimum and maximum values which are avialable directly from the summary dictionary. 

medv = summary['median']
minv = summary['minimum']['value']
maxv = summary['maximum']['value']
print (minv,medv,maxv)

# From the `percentiles` values in the summary dictionary we can generate more values for the box and whisker plot. Where `lopcp` and `hipcp` are the low and high percentage values that we set in the URL (in this example, 10 and 75), `lopcv` and `hipcv` are the low and high percentage values, `first` and `last` are the first and last timestamp values from the summary dictionary. We can also print all of these variables and see that they match up with the `summary` information above. 

# +
pcs = summary['percentiles']

lopcp = pcs[0].get('perc')
hipcp = pcs[1].get('perc')

lopcv = pcs[0].get('value')
hipcv = pcs[1].get('value')

first = summary['first']['ts']
last = summary['last']['ts']

print (pcs)
print (lopcp,hipcp,lopcv,hipcv,first,last)
# -

# We will now plot this using a box-and-whisker style by using our above parameters. 

# +
fig,ax = plt.subplots()
#draw box from lowpcv to hipcv
ax.bar(0, hipcv-lopcv, bottom=lopcv, color='#41b0f0', label='precentiles '+str(lopcp)+'-'+str(hipcp))
# vertical line from minv to maxv
ax.vlines(0, minv, maxv, color='black', label='min-max')
#horizontal line at medv
ax.hlines(medv, -0.4, 0.4, linestyles='dashed', colors='red', label='median')

#make more like a box-and-whisker
ax.set_xlim(-1,1)
ax.ticklabel_format(useOffset=False, style='plain')
plt.xticks([])

plt.legend(loc='best')

plt.suptitle('Data summary '+my_series['station']+'-'+my_series['name']+'-'+my_series['sensorCode']+'-'+my_series['type'])
plt.title(first+' - '+last)
# -

# ## &nbsp;6. Comparing query results from several stations
#
# We now make a query for several stations and extract the relevant data from each query, then use this data to demonstrate how to calculate the water level range for several sites over the time period queried. This is not completely comprehensive as we are restricting our query to a sensor code of `40`.
#
# We need to create a list of stations using a data summary query, and then loop through that list, making an API call for each, and accumulating the results.
#
# Our query is for the time period from `2021-09-01` to `2021-09-10`.

# +
#make a query to get a list of all possible stations to query
url = 'https://tilde.geonet.org.nz/v1/dataSummary/dart'
r = requests.get(url)
summary = r.json()

stations = summary['domain']['dart']['stations']
stacodes = list(stations.keys())
# -

# The next sub-codes uses a python module to create the station list with corresponding height range. These are the inputs used in the subsequent plot.

# +
#make a dataframe for later use
df = pd.DataFrame()

#dates for query
d1 = '2021-09-01'
d2 = '2021-09-10'

for code in stacodes:
    url = source+'/v1/stats/dart/'+code+'/water-height/40/WTZ/'+d1+'/'+d2
    r = requests.get(url)
    summary = r.json()
    
    #if maximum and minimum have values, calculate range and append to dataframe
    if ((summary['maximum'])and(summary['minimum'])):
        range = summary['maximum']['value'] - summary['minimum']['value']
        #make a dictionary of the extracted data
        d = {'range':range, 'station':code}
        #append this to the dataframe
        df = df.append(d, ignore_index=True)

#make station the dataframe index
df.set_index('station', inplace=True)

# +
fig,ax = plt.subplots(figsize=(15,5))
#plot in metres, not mm
ax.plot(df['range']/1000, color='#41b0f0', marker='o', linestyle='none')
    
ax.set_ylabel('Height range (m)')
ax.set_xlabel('DART')

plt.suptitle('Height Range at DART Buoys', fontsize=16)
plt.title(d1+' - '+d2)

StarcoderdataJupyterScriptsDedupFiltered