{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Elaborate statistics features for Silvereye\n",
"\n",
"## Dependencies imports\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import xarray as xr\n",
"import os\n",
"import sys\n",
"import pandas as pd\n",
"from functools import wraps\n",
"import numpy as np\n",
"\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import seaborn as sns # noqa, pandas aware plotting library"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from datetime import date\n",
"from dateutil.relativedelta import relativedelta # $ pip install python-dateutil"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The last section of this notebook investigates ipyleaflet for visualisation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from ipyleaflet import Map, basemaps, basemap_to_tiles, ImageOverlay"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import PIL"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from io import StringIO, BytesIO\n",
"from base64 import b64encode"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if ('SP_SRC' in os.environ):\n",
" root_src_dir = os.environ['SP_SRC']\n",
"elif sys.platform == 'win32':\n",
" root_src_dir = r'C:\\src\\csiro\\stash\\silverpieces'\n",
"else:\n",
" root_src_dir = '/silverpieces'\n",
"\n",
"pkg_src_dir = root_src_dir\n",
"sys.path.append(pkg_src_dir)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from silverpieces import *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from silverpieces.functions import *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"if ('SP_DATA' in os.environ):\n",
" root_data_dir = os.environ['SP_DATA']\n",
"elif sys.platform == 'win32':\n",
" root_data_dir = r'C:\\data\\silverpieces'\n",
"else:\n",
" root_data_dir = '/silverpieces/notebooks/'\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'retina'\n",
"# the default cmap_sequential for xarray is viridis. 'RdBu' is divergent, but works better for wetness concepts\n",
"# # https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html\n",
"xr.set_options(cmap_sequential='bwr_r')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Can get tassie_silo_rain.nc data from https://cloudstor.aarnet.edu.au/plus/s/nj2RevvD1EUD77n\n",
"fn = os.path.join(root_data_dir, 'tassie_silo_rain.nc')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tassie = xr.open_mfdataset([fn])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tassie"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"daily_rain = tassie.daily_rain"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"daily_rain.isel(time=4300).plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Use case - inter-period statistical comparisons.\n",
"\n",
"We want to be able to compare a grid of statistics for a period compared to all periods of similar lengths.\n",
"The start and end of the period should be as arbitrary as possible. The sliding window could however be limited or fixed to a year: it is probably moot to compare windows with shifted seasonality. \n",
"\n",
"### How does the cumulated rainfall 2016-2018 over TAS compare with all 3 year periods over the whole record?\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"s = SpatialTemporalDataArrayStat()\n",
"start_time = pd.to_datetime('2016-01-01')\n",
"end_time = pd.to_datetime('2018-12-31')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"daily_rain.load()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"three_years_rains = s.periods_stat_yearly(daily_rain, start_time, end_time, func = np.sum)\n",
"three_years_rains.name = '3yrs cumulated rain'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"TIME_DIMNAME = 'time'\n",
"g_simple = three_years_rains.plot(x='lon', y='lat', col=TIME_DIMNAME, col_wrap=3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's define \"percentiles\" of interest as boundaries to classify against "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"q = np.array([.1, .5, .9])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y = s.quantile_over_time_dim(three_years_rains, q=q)\n",
"\n",
"y.name = '3-yr cumulated rain quantiles'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following color scheme may not be the best to see the map of quantile values, but should give an idea "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y.plot(x='lon', y='lat', col='quantile', col_wrap=3, cmap='gist_ncar')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So now we want a map that tells us where the last three years, for every grid cell, sits (which side of every quantile)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"last_three_years_cumrain = three_years_rains[-1,:,:]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat_q = s.searchsorted(y, last_three_years_cumrain)\n",
"cat_q.name = 'Quantile categories 10/50/90'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat_q.plot(cmap='bwr_r')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So, the three years 2016 to 2018 have been the wettest on most of the mountainous areas of the state compared to the last decade, except for the south west National Park which has been the driest comparatively.\n",
"\n",
"That said, the deviation from mean may still be quite small and it may not be a \"drought\" as such.\n",
"\n",
"Now let's look at inter-annual variability rather than 3-year moving windows."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat_q"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat_q.values.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearly_rain = s.periods_stat_yearly(daily_rain, '2016-01-01', '2016-12-31', func = np.sum) # Yes, xarray vanilla would suffice in this case.\n",
"yearly_rain.name = 'yearly rainfall'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearly_rain.plot(x='lon', y='lat', col=TIME_DIMNAME, col_wrap=3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearly_cat_q = yearly_rain.copy()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y = s.quantile_over_time_dim(yearly_rain, q=q)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for yr in range(len(yearly_rain.time)):\n",
" x = yearly_rain[yr,:,:]\n",
" cat_q = s.searchsorted(y, x)\n",
" yearly_cat_q[yr,:,:] = cat_q"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearly_cat_q.name = 'Quantile categories yearly rain 10/50/90'\n",
"yearly_cat_q.plot(x='lon', y='lat', col=TIME_DIMNAME, col_wrap=3, cmap='bwr_r')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ipyleaflet"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from silverpieces.vis import *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"bounds = make_bounds(cat_q)\n",
"imgurl = to_embedded_png(cat_q)\n",
"\n",
"io = ImageOverlay(url=imgurl, bounds=bounds)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"center = center_from_bounds(bounds)\n",
"zoom = 7\n",
"m = Map(center=center, zoom=zoom, interpolation='nearest')\n",
"m.layout.height = '600px'\n",
"\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.add_layer(io)\n",
"io.interact(opacity=(0.0,1.0,0.01))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try to add a colorscale legend."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from branca.colormap import linear "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Could not get something that displays in the widget from matplotlib colormaps\n",
"legend = linear.RdBu_09.scale(0,3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.cm.bwr_r "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"io = ImageOverlay(url=imgurl, bounds=bounds)\n",
"io.colormap=legend"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from ipywidgets.widgets import Output"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"out = Output(layout={'border': '1px solid black'})\n",
"with out:\n",
" display(legend)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from ipyleaflet import WidgetControl"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = Map(center=center, zoom=zoom, interpolation='nearest')\n",
"m.layout.height = '600px'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.add_layer(io)\n",
"io.interact(opacity=(0.0,1.0,0.01))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"widget_control = WidgetControl(widget=out, position='topright')\n",
"m.add_control(widget_control)\n",
"display(m)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Py3 (sv)",
"language": "python",
"name": "dcx"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}