<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 SC57 - Working with big, multi-dimensional geoscientific datasets in Python
Step1:  Basic data arrays in numpy
Step2:  numpy is a powerful but "low-level" array manipulation tool. Axis only have numbers and no names (it is easy to forget which axis is what, a common source of trivial bugs), arrays can't carry metadata (e.g. units), and the data is unstructured (i.e. the coordinates and/or other related arrays have to be handled separately
Step3:  Let's Do Some Math
Step4:  Notice xarray has changed the colormap according to the dataset (borrowing logic from Seaborn).
* With degrees C, the data passes through 0, so a diverging colormap is used
* With Kelvin, the default colormap is used.
Step5:  Adding Data to DataSets
Step6:  Let's add those kelvin and centigrade dataArrays to the dataset.
Step7:  3. Selecting data with named dimensions
In xarray there are many different ways for selecting and indexing data.
Positional indexing (old way)
This is the "old way", i.e. like numpy
Step8:  This selection implies prior knowledge about the structure of the data, and is therefore much less readable than the "xarray methods" presented below.
Selection by index
Selection based on the index of a coordinate
Step9:  Selection by value
Selection based on the value of a coordinate
Step10:  Selection by value works well for time, too
Step11:  Selecting a range of values
The syntax is similar, but you'll need to use a slice
Step12:  Nearest neighbor lookup
Step13:  4. Operations and computation
We can do arithmetic directly on Dataset and DataArray objects. 
Labels are preserved and dataArray dimensions automatically aligned.
Broadcasting
<img src="./figures/broadcast.png" width="50%" align="left">
Step14:  Alignment
<img src="./figures/align.png" width="50%" align="left">
Step15:  Aggregation
Step16:  Masking with .where()
Step17:  5. Groupby and "split-apply-combine"
Xarray implements the "split-apply-combine" paradigm with groupby. This works really well for calculating climatologies
Step18:  <img src="./figures/split_apply-combine.png" alt="split" style="width
Step19:  You can also do arithmetic with groupby objects, which repeats the arithmetic over each group
Step20:  Resample adjusts a time series to a new resolution
Step21:  6. Graphics
xarray plotting functions rely on matplotlib internally, but they make use of all available metadata to make the plotting operations more intuitive and interpretable.
1D plots
Step22:  2D plots
Step23:  Customizing 2d plots
Step24:  Dealing with Outliers
Step25:  Facet plots
Step26:  Plotting on maps
For plotting on maps, we rely on the excellent cartopy library.
Step27:  Facet plots on maps
Step28:  Seaborn is Cool
Statistical visualization with Seaborn
Step29:  7. Out-of-core computation
Here's a quick demo of how xarray can leverage dask to work with data that doesn't fit in memory. This lets xarray substitute for tools like cdo and nco.
Let's open 10 years of runoff data
xarraycan open multiple files at once using string pattern matching.  
In this case we open all the files that match our filestr, i.e. all the files for the 2080s. 
Each of these files (compressed) is approximately 80 MB.
PS - these files weren't available during the tutorial. The data we used was daily discharge hydrological data from the ISIMIP project (e.g. HadGEM2-ES / PCRGLOBWB / RCP2p6), which we cannot share here but is available for download.
Step30:  xarray even puts them in the right order for you.
Step31:  How big is all this data uncompressed? Will it fit into memory?
Step32:  Working with Big Data
This data is too big for our memory. 
That means we need to process it in chunks. 
We can do this chunking in xarray very easily.  
xarray computes data 'lazily'. That means that data is only loaded into memory when it is actually required. This also allows us to inspect datasets without loading all the data into memory.  
To do this xarray integrates with dask to support streaming computation on datasets that don’t fit into memory.
<img src="./figures/dask-array.png" width="40%" align="center">
<ASSISTANT_TASK:>
Python Code:
# standard imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
import warnings
%matplotlib inline
np.set_printoptions(precision=3, linewidth=80, edgeitems=1)  # make numpy less verbose
xr.set_options(display_width=70)
warnings.simplefilter('ignore') # filter some warning messages
Explanation: SC57 - Working with big, multi-dimensional geoscientific datasets in Python: a tutorial introduction to xarray
Original notebook by Stephan Hoyer, Rossbypalooza, 2016.
Modified by Edward Byers, Matthew Gidden and Fabien Maussion for EGU General Assembly 2017, Vienna, Austria
Thursday, 27th April, 15:30–17:00 / Room -2.91  
Convenors
* Dr Edward Byers    - International Institute for Applied Systems Analysis, Laxenburg, Austria
* Dr Matthew Gidden  - International Institute for Applied Systems Analysis, Laxenburg, Austria
* Dr Fabien Maussion - University of Innsbruck, Innsbruck, Austria
With
you can reach
Structure of this tutorial
Introduction to key features of xarray
Basic operations in xarray: opening, inspecting, selecting and indexing data
Selecting data with named dimensions
Operations and computation
Groupby and "split-apply-combine"
Graphics
Out-of-core computation
1. Key features of xarray
What is xarray?
xarray is an open source project and Python package
xarray has been designed to perform labelled data analysis on multi-dimensional arrays
the xarray approach adopts the Common Data Model for self-describing scientific data in widespread use in the Earth sciences
xarray.Dataset is an in-memory representation of a netCDF file.
xarray is built on top of the dataprocessing library Pandas (the best way to work with tabular data (e.g., CSV files) in Python)
Our data
<img src="./figures/dataset.png" width="50%" align="right"> 
numeric
multi-dimensional
labelled
(lots of) metadata
sometimes (very) large
What is xarray good for?
Gridded, multi-dimensional and large datasets, commonly used in earth sciences, but also increasingly finance, engineering (signal/image processing), and biological sciences
Integration with other data analysis packages such as Pandas 
I/O operations (NetCDF)
Plotting
Out of core computation and parallel processing
Extensions based on xarray
...
Where can I find more info?
For more information about xarray
Read the online documentation
Ask questions on StackOverflow
View the source code and file bug reports on GitHub
For more doing data analysis with Python:
Thomas Wiecki, A modern guide to getting started with Data Science and Python
Wes McKinney, Python for Data Analysis (book)
Packages building on xarray for the geophysical sciences
For analyzing GCM output:
xgcm by Ryan Abernathey
oogcm by Julien Le Sommer
MPAS xarray by Phil Wolfram
marc_analysis by Daniel Rothenberg
Other tools:
windspharm: wind spherical harmonics by Andrew Dawson
eofs: empirical orthogonal functions by Andrew Dawson
infinite-diff by Spencer Hill 
aospy by Spencer Hill and Spencer Clark
regionmask by Mathias Hauser
salem by Fabien Maussion
Resources for teaching and learning xarray in geosciences:
- Fabien's teaching repo: courses that combine teaching climatology and xarray
2. Basic operations in xarray
Import python packages
End of explanation
import numpy as np
a = np.array([[1, 3, 9], [2, 8, 4]])
a
a[1, 2]
a.mean(axis=0)
Explanation: Basic data arrays in numpy
End of explanation
ds = xr.tutorial.load_dataset('air_temperature')
ds
ds.air
ds.dims
ds.attrs
ds.air.values
type(ds.air.values)
ds.air.dims
ds.air.attrs
ds.air.attrs['tutorial-date'] = 27042017
ds.air.attrs
Explanation: numpy is a powerful but "low-level" array manipulation tool. Axis only have numbers and no names (it is easy to forget which axis is what, a common source of trivial bugs), arrays can't carry metadata (e.g. units), and the data is unstructured (i.e. the coordinates and/or other related arrays have to be handled separately: another source of bugs).
This is where xarray comes in!
Properties of xarray.Dataset and xarray.DataArray objects
We'll start with the "air_temperature" tutorial dataset. This tutorial comes with the xarray package. Other examples here.
End of explanation
kelvin = ds.air.mean(dim='time')
kelvin.plot();
centigrade = kelvin - 273.16
centigrade.plot();
Explanation: Let's Do Some Math
End of explanation
# ufuncs work too
np.sin(centigrade).plot();
Explanation: Notice xarray has changed the colormap according to the dataset (borrowing logic from Seaborn).
* With degrees C, the data passes through 0, so a diverging colormap is used
* With Kelvin, the default colormap is used.
End of explanation
ds
Explanation: Adding Data to DataSets
End of explanation
ds['centigrade'] = centigrade
ds['kelvin'] = kelvin
ds
ds.kelvin.attrs  # attrs are empty! Let's add some
ds.kelvin.attrs['Description'] = 'Mean air tempterature (through time) in kelvin.'
ds.kelvin
ds.to_netcdf('new file.nc')
Explanation: Let's add those kelvin and centigrade dataArrays to the dataset.
End of explanation
ds.air[:, 1, 2]  # note that the attributes, coordinates are preserved
ds.air[:, 1, 2].plot();
Explanation: 3. Selecting data with named dimensions
In xarray there are many different ways for selecting and indexing data.
Positional indexing (old way)
This is the "old way", i.e. like numpy:
End of explanation
ds.air.isel(time=0).plot();  # like above, but with a dimension name this time
Explanation: This selection implies prior knowledge about the structure of the data, and is therefore much less readable than the "xarray methods" presented below.
Selection by index
Selection based on the index of a coordinate:
End of explanation
ds.air.sel(lat=72.5, lon=205).plot(); 
Explanation: Selection by value
Selection based on the value of a coordinate:
End of explanation
ds.air.sel(time='2013-01-02').plot(); # Note that we will extract 4 time steps! 3d data is plotted as histogram
ds.air.sel(time='2013-01-02T06:00').plot();  # or look at a single timestep
Explanation: Selection by value works well for time, too
End of explanation
ds.air.sel(lat=slice(60, 50), lon=slice(200, 270), time='2013-01-02T06:00:00').plot();
Explanation: Selecting a range of values
The syntax is similar, but you'll need to use a slice:
End of explanation
ds.air.sel(lat=41.8781, lon=360-87.6298, method='nearest', tolerance=5).plot();
Explanation: Nearest neighbor lookup
End of explanation
a = xr.DataArray(np.arange(3), dims='time', 
                 coords={'time':np.arange(3)})
b = xr.DataArray(np.arange(4), dims='space', 
                 coords={'space':np.arange(4)})
a + b
Explanation: 4. Operations and computation
We can do arithmetic directly on Dataset and DataArray objects. 
Labels are preserved and dataArray dimensions automatically aligned.
Broadcasting
<img src="./figures/broadcast.png" width="50%" align="left">
End of explanation
atime = np.arange(3)
btime = np.arange(5) + 1
atime, btime
a = xr.DataArray(np.arange(3), dims='time', 
                 coords={'time':atime})
b = xr.DataArray(np.arange(5), dims='time', 
                 coords={'time':btime})
a + b
Explanation: Alignment
<img src="./figures/align.png" width="50%" align="left">
End of explanation
ds.max()
ds.air.median(dim=['lat', 'lon']).plot();
Explanation: Aggregation
End of explanation
means = ds.air.mean(dim=['time'])
means.where(means > 273.15).plot();
Explanation: Masking with .where()
End of explanation
ds.air.groupby('time.season').mean()
Explanation: 5. Groupby and "split-apply-combine"
Xarray implements the "split-apply-combine" paradigm with groupby. This works really well for calculating climatologies:
End of explanation
ds.air.groupby('time.month').mean('time')
clim = ds.air.groupby('time.month').mean('time')
Explanation: <img src="./figures/split_apply-combine.png" alt="split" style="width: 50%;"/>
End of explanation
anomalies = ds.air.groupby('time.month') - clim
anomalies
anomalies.plot();
anomalies.sel(time= '2013-02').plot();  # Find all the anomolous values for February
Explanation: You can also do arithmetic with groupby objects, which repeats the arithmetic over each group:
End of explanation
tmin = ds.air.resample('1D', dim='time', how='min')  # Resample to one day '1D
tmax = ds.air.resample('1D', dim='time', how='max')
(tmin.sel(time='2013-02-15') - 273.15).plot();
ds_extremes = xr.Dataset({'tmin': tmin, 'tmax': tmax})
ds_extremes
Explanation: Resample adjusts a time series to a new resolution:
End of explanation
zonal_t_average = ds.air.mean(dim=['lon', 'time']) - 273.15
zonal_t_average.plot();  # 1D arrays are plotted as line plots
Explanation: 6. Graphics
xarray plotting functions rely on matplotlib internally, but they make use of all available metadata to make the plotting operations more intuitive and interpretable.
1D plots
End of explanation
t_average = ds.air.mean(dim='time') - 273.15
t_average.plot();  # 2D arrays are plotted with pcolormesh
t_average.plot.contourf();  # but you can use contour(), contourf() or imshow() if you wish
Explanation: 2D plots
End of explanation
t_average.plot.contourf(cmap='BrBG_r', vmin=-15, vmax=15);
t_average.plot.contourf(cmap='BrBG_r', levels=22, center=False);
Explanation: Customizing 2d plots
End of explanation
air_outliers = ds.air.isel(time=0).copy()
air_outliers[0, 0] = 100
air_outliers[-1, -1] = 400
air_outliers.plot();  # outliers mess with the datarange and colorscale!
# Using `robust=True` uses the 2nd and 98th percentiles of the data to compute the color limits.
air_outliers.plot(robust=True);
Explanation: Dealing with Outliers
End of explanation
t_season = ds.air.groupby('time.season').mean(dim='time') - 273.15
# facet plot allows to do multiplot with the same color mappings
t_season.plot.contourf(x='lon', y='lat', col='season', col_wrap=2, levels=22);
Explanation: Facet plots
End of explanation
import cartopy.crs as ccrs
f = plt.figure(figsize=(8, 4))
# Define the map projection *on which* you want to plot
ax = plt.axes(projection=ccrs.Orthographic(-80, 35))
# ax is an empty plot. We now plot the variable t_average onto ax
# the keyword "transform" tells the function in which projection the air temp data is stored 
t_average.plot(ax=ax, transform=ccrs.PlateCarree())
# Add gridlines and coastlines to the plot
ax.coastlines(); ax.gridlines(); 
Explanation: Plotting on maps
For plotting on maps, we rely on the excellent cartopy library.
End of explanation
# this time we need to retrieve the plots to do things with the axes later on
p = t_season.plot(x='lon', y='lat', col='season', transform=ccrs.PlateCarree(),
                  subplot_kws={'projection': ccrs.Orthographic(-80, 35)})
for ax in p.axes.flat:
    ax.coastlines()
Explanation: Facet plots on maps
End of explanation
import seaborn as sns
data = (ds_extremes
        .sel_points(lat=[41.8781, 37.7749], lon=[360-87.6298, 360-122.4194],
                    method='nearest', tolerance=3,
                    dim=xr.DataArray(['Chicago', 'San Francisco'],
                                     name='location', dims='location'))
        .to_dataframe()
        .reset_index()
        .assign(month=lambda x: x.time.dt.month))
plt.figure(figsize=(10, 5))
sns.violinplot('month', 'tmax', 'location', data=data, split=True, inner=None);
Explanation: Seaborn is Cool
Statistical visualization with Seaborn:
End of explanation
from glob import glob
files = glob('data/*dis*.nc')
runoff = xr.open_mfdataset(files)
runoff
Explanation: 7. Out-of-core computation
Here's a quick demo of how xarray can leverage dask to work with data that doesn't fit in memory. This lets xarray substitute for tools like cdo and nco.
Let's open 10 years of runoff data
xarraycan open multiple files at once using string pattern matching.  
In this case we open all the files that match our filestr, i.e. all the files for the 2080s. 
Each of these files (compressed) is approximately 80 MB.
PS - these files weren't available during the tutorial. The data we used was daily discharge hydrological data from the ISIMIP project (e.g. HadGEM2-ES / PCRGLOBWB / RCP2p6), which we cannot share here but is available for download.
End of explanation
runoff.time
Explanation: xarray even puts them in the right order for you.
End of explanation
runoff.nbytes / 1e9  # Convert to gigiabytes
Explanation: How big is all this data uncompressed? Will it fit into memory?
End of explanation
runoff = runoff.chunk({'lat': 60})
runoff.chunks
%time ro_seasonal = runoff.groupby('time.season').mean('time')
import dask
from multiprocessing.pool import ThreadPool
dask.set_options(pool=ThreadPool(1))
%time ro_seasonal.compute()
dask.set_options(pool=ThreadPool(4))
%time ro_seasonal = runoff.groupby('time.season').mean('time')
%time result = ro_seasonal.compute()
brazil = dict(lat=slice(10.75, -40.75), lon=slice(-100.25, -25.25))
result.dis.sel(**brazil).plot(col='season', size=4, cmap='Spectral_r')
Explanation: Working with Big Data
This data is too big for our memory. 
That means we need to process it in chunks. 
We can do this chunking in xarray very easily.  
xarray computes data 'lazily'. That means that data is only loaded into memory when it is actually required. This also allows us to inspect datasets without loading all the data into memory.  
To do this xarray integrates with dask to support streaming computation on datasets that don’t fit into memory.
<img src="./figures/dask-array.png" width="40%" align="center">
End of explanation
<END_TASK>