In [1]:
from notebook.services.config import ConfigManager
cm = ConfigManager()
cm.update('livereveal',
{'theme': 'serif',
'transition': 'slide',
'start_slideshow_at': 'selected',
'width': 1024,
'height': 768});
In [2]:
import warnings
from IPython.display import HTML
warnings.filterwarnings("ignore")
with open('./custom.css') as f:
css = f.read()
ids = """<li>
<a href="https://github.com/ocefpaf">
<span class="icon github">
<svg version="1.1" class="github-icon-svg" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px"
viewBox="0 0 16 16" enable-background="new 0 0 16 16" xml:space="preserve">
<path fill-rule="evenodd" clip-rule="evenodd" fill="#C2C2C2" d="M7.999,0.431c-4.285,0-7.76,3.474-7.76,7.761
c0,3.428,2.223,6.337,5.307,7.363c0.388,0.071,0.53-0.168,0.53-0.374c0-0.184-0.007-0.672-0.01-1.32
c-2.159,0.469-2.614-1.04-2.614-1.04c-0.353-0.896-0.862-1.135-0.862-1.135c-0.705-0.481,0.053-0.472,0.053-0.472
c0.779,0.055,1.189,0.8,1.189,0.8c0.692,1.186,1.816,0.843,2.258,0.645c0.071-0.502,0.271-0.843,0.493-1.037
C4.86,11.425,3.049,10.76,3.049,7.786c0-0.847,0.302-1.54,0.799-2.082C3.768,5.507,3.501,4.718,3.924,3.65
c0,0,0.652-0.209,2.134,0.796C6.677,4.273,7.34,4.187,8,4.184c0.659,0.003,1.323,0.089,1.943,0.261
c1.482-1.004,2.132-0.796,2.132-0.796c0.423,1.068,0.157,1.857,0.077,2.054c0.497,0.542,0.798,1.235,0.798,2.082
c0,2.981-1.814,3.637-3.543,3.829c0.279,0.24,0.527,0.713,0.527,1.437c0,1.037-0.01,1.874-0.01,2.129
c0,0.208,0.14,0.449,0.534,0.373c3.081-1.028,5.302-3.935,5.302-7.362C15.76,3.906,12.285,0.431,7.999,0.431z"/>
</svg>
</span>
<span class="username">ocefpaf</span>
</a>
</li>
<li>
<a href="https://twitter.com/ocefpaf">
<span class="icon twitter">
<svg version="1.1" class="twitter-icon-svg" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" x="0px" y="0px"
viewBox="0 0 16 16" enable-background="new 0 0 16 16" xml:space="preserve">
<path fill="#C2C2C2" d="M15.969,3.058c-0.586,0.26-1.217,0.436-1.878,0.515c0.675-0.405,1.194-1.045,1.438-1.809
c-0.632,0.375-1.332,0.647-2.076,0.793c-0.596-0.636-1.446-1.033-2.387-1.033c-1.806,0-3.27,1.464-3.27,3.27
c0,0.256,0.029,0.506,0.085,0.745C5.163,5.404,2.753,4.102,1.14,2.124C0.859,2.607,0.698,3.168,0.698,3.767
c0,1.134,0.577,2.135,1.455,2.722C1.616,6.472,1.112,6.325,0.671,6.08c0,0.014,0,0.027,0,0.041c0,1.584,1.127,2.906,2.623,3.206
C3.02,9.402,2.731,9.442,2.433,9.442c-0.211,0-0.416-0.021-0.615-0.059c0.416,1.299,1.624,2.245,3.055,2.271
c-1.119,0.877-2.529,1.4-4.061,1.4c-0.264,0-0.524-0.015-0.78-0.046c1.447,0.928,3.166,1.469,5.013,1.469
c6.015,0,9.304-4.983,9.304-9.304c0-0.142-0.003-0.283-0.009-0.423C14.976,4.29,15.531,3.714,15.969,3.058z"/>
</svg>
</span>
<span class="username">ocefpaf</span>
</a>
</li>
<li>
<a href="mailto:ocefpaf@gmail.com">ocefpaf@gmail.com</a>
</li>
"""
In [3]:
HTML(ids)
Out[3]:
One answer to all questions is to use the Anaconda Python distribution:
url=http://bit.ly/miniconda
wget $url -O miniconda.sh
bash miniconda.sh -b
export PATH=$HOME/miniconda/bin:$PATH
conda update --yes --all
url=http://bit.ly/conda-env
wget $url -O environment.yml
conda env create environment.yml
Reading some data
In [4]:
!head -n 5 ./data/dados_pirata.csv
We need to skip the header and the date-time columns
In [5]:
import numpy as np
data = np.loadtxt("./data/dados_pirata.csv",
skiprows=1,
usecols=range(2, 16),
delimiter=',')
data[data == -99999.] = np.NaN
data.shape, data.dtype
Out[5]:
But the depth information was in header!
We need to type that in...
In [6]:
%matplotlib inline
import matplotlib.pyplot as plt
z = [1, 10, 100, 120, 13, 140, 180, 20,
300, 40,5, 500, 60, 80]
fig, ax = plt.subplots()
# Points because z was not ordered.
ax.plot(data[42, :], z, 'ko')
ax.invert_yaxis()
Note that almost everything we did before, and some more, is now done at the loading time
In [7]:
import pandas as pd
df = pd.read_csv('./data/dados_pirata.csv',
index_col='datahora',
parse_dates=True,
na_values=-99999)
df.drop('Unnamed: 0', axis=1, inplace=True)
df.columns = ['{0:0>3}'.format(col.split('_')[1]) for
col in df.columns]
df.sort(axis=1, inplace=True)
df.head(3)
Out[7]:
Pandas behave like masked-arrays,
but with some extra convenience functions beyond nan{min,max,mean,etc}
In [8]:
df.dropna(how='all', axis=1, inplace=True)
df.head(3)
Out[8]:
Quick way to explore the data
In [9]:
desc = df.describe()
desc
Out[9]:
High level code! No need for comments to describe what is being computed below:
In [10]:
desc.ix['std'] ** 2
Out[10]:
What levels are we plotting?
In [11]:
ax = df[['001', '500']].plot(figsize=(11, 3))
More convenience functions:
In [12]:
df['001'].interpolate().plot(figsize=(11, 3))
ax = df['001'].plot()
Full control of the interpolation
In [13]:
kw = dict(method='time', limit=20)
df['001'].interpolate(**kw).plot(figsize=(11, 3))
ax = df['001'].plot()
Pandas loves time-series
In [14]:
key = lambda x: x.month
grouped = df.groupby(key)
monthly = grouped.mean()
fig, ax = plt.subplots(figsize=(9, 5))
ax = monthly.plot(ax=ax)
GeoPandas read GeoJSONs, Shapefiles, etc as pandas DataFrames
In [15]:
import rasterio
import geopandas as gpd
fname = "./data/2013-04-29-Running.geojson"
df = gpd.read_file(fname)
ax = df.plot()
The geometries are essentially shapely objects
In [16]:
shape = df['geometry'][0]
msg = """
Length is {} km.
The center is locate at {!r}.
The track bounds are {}.
""".format
print(msg(shape.length * 111, shape.centroid.xy,
shape.bounds))
ptyhon-ctd loads ctd, xbt, and many more ocean profilers directly into pandas
In [17]:
from ctd import DataFrame
cast = DataFrame.from_cnv('./data/CTD_001.cnv.gz', compression='gzip')
cast.head()
Out[17]:
The augmented DataFrame has some extra methods to process CTD data
In [18]:
downcast, upcast = cast['t090C'].split()
fig, ax = downcast.plot(figsize=(3, 7))
xray: netCDF4 + pandas for >2D tables
In [19]:
import xray
url = ('http://geoport.whoi.edu/thredds/dodsC/clay/usgs/'
'users/jcwarner/Projects/Sandy/triple_nest/'
'00_dir_NYB05.ncml')
ds = xray.open_dataset(url)
ds
Out[19]:
Slice the variables using the var name (like netCDF4)
In [20]:
temp = ds['temp']
temp
Out[20]:
Use high-level dates slice (like pandas)
In [21]:
temp.sel(ocean_time='2012-10-25')
Out[21]:
All pandas time-series goodies are present
temp.resample('1M', 'ocean_time')
You can even slice using integers and named the axis
In [22]:
t = temp.isel(s_rho=-1, ocean_time=0)
cs = t.plot()
http://dask.pydata.org/en/latest/
https://jakevdp.github.io/blog/2015/08/14/out-of-core-dataframes-in-python/
In [23]:
import dask
from dask.async import get_sync
dask.set_options(get=get_sync)
Out[23]:
The following lines do not download the data
In [24]:
from netCDF4 import Dataset
nc = Dataset(url)
t = nc['temp']
s = nc['salt']
t.shape, s.shape
Out[24]:
By making the netCDF4 variable a dask variable we can start computing stuff without loading the data! Note
In [25]:
import dask.array as da
chunks = (1, 16, 107, 345)
t = da.from_array(t, chunks=chunks)
s = da.from_array(s, chunks=chunks)
t
Out[25]:
Note that the result is a dask element wise object
In [26]:
(s + t) * 2
Out[26]:
You need to manually trigger the computations to load the data. Note that some operations, like plotting, will load the data automatically:
In [27]:
import gsw
import seawater as sw
import numpy.ma as ma
sigma0 = gsw.sigma0(s[-1, :, :, :], t[-1, :, :, :])
cs = plt.pcolormesh(ma.masked_equal(sigma0, -1000))
The first python library to actually solve the dateline problem
In [28]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
projection = ccrs.PlateCarree()
ax = plt.axes(projection=projection)
ax.coastlines(); ax.set_global()
kw = dict(lw=4, color='g', transform=ccrs.Geodetic())
l0 = plt.plot([-100, 50], [25, 25], label='GD1', **kw)
l1 = plt.plot([-38, 147], [-13, -42], label='GD2', **kw)
kw = dict(linewidth=4, color='b', transform=projection)
l2 = plt.plot([-100, 50], [25, 25], label='PC1', **kw)
l3 = plt.plot([-38, 147], [-13, -42], label='PC2', **kw)
leg = ax.legend(loc=(1.05, 0.5))
When creating a lot of maps it is convenient to wrap them in a plotting function
In [29]:
from oceans import cm
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import (LONGITUDE_FORMATTER,
LATITUDE_FORMATTER)
def make_map(projection=ccrs.PlateCarree()):
subplot_kw = dict(projection=projection)
fig, ax = plt.subplots(figsize=(9, 13),
subplot_kw=subplot_kw)
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
In [30]:
lon = nc['lon_rho'][:]
lat = nc['lat_rho'][:]
fig, ax = make_map()
ax.set_extent([lon.min(), lon.max(), lat.min(), lat.max()])
ax.coastlines('10m')
s0 = ma.masked_equal(sigma0, -1000)
cs = ax.pcolormesh(lon, lat, s0, cmap=cm.avhrr)
cbar = fig.colorbar(cs, shrink=0.45, extend='both')
Iris is an interpretation of the CF-Conventions.
The main object is the cube:
In [31]:
import iris
cubes = iris.load_raw(url)
print(cubes)
We can take advantage of the CF conventions to access variables by its standard_name
In [32]:
temp = cubes.extract_strict('sea_water_potential_temperature')
print(temp)
Like xray we have both high and low level slicing
(Although the high level slicing in iris is too sophisticated for my taste ;-)
In [33]:
t = temp[-1, -1, ...]
print(t) # The metadata is always propagated.
There are even some quick plotting routines to explore the data. Note the free colorbar, units label, and title:
In [34]:
import iris.quickplot as qplt
cs = qplt.pcolormesh(t)
Calling matlab
(See matlab_kernel for full matlab)
In [35]:
try:
%load_ext oct2py.ipython
octave = True
except AttributeError as e:
octave = False
print(e)
print("You cannot run the next 2 cells with the octave example :-(")
In [36]:
%%file fahr_to_kelvin.m
function ktemp = fahr_to_kelvin(ftemp)
ktemp = ((ftemp - 32) * (5/9)) + 273.15;
end
In [37]:
if octave:
%octave fahr_to_kelvin(32)
Calling Fortran
In [38]:
%load_ext fortranmagic
In [39]:
%%fortran
subroutine calc_sin(x, y)
real, intent(in) :: x
real, intent(out) :: y
y = sin(x)
end subroutine calc_sin
In [40]:
calc_sin(np.pi/2)
Out[40]:
Calling R
In [41]:
%load_ext rpy2.ipython
In [42]:
%%R
X=c(1,4,5,7)
Y = c(2,4,3,9)
print(summary(lm(Y~X)))
In [43]:
%%R
library(oce)
Sa <- 30; Ta <- 10; Sb <- 40;
rho0 <- swRho(Sa, Ta, 0)
Tb <- uniroot(function(T) rho0-swRho(Sb,T,0), lower=0, upper=100)$root
Sc <- (Sa + Sb) /2; Tc <- (Ta + Tb) /2
drho <- swRho(Sc, Tc, 0) - rho0
dT <- drho / rho0 / swAlpha(Sc, Tc, 0)
plotTS(as.ctd(c(Sa, Sb, Sc), c(Ta, Tb, Tc), 0), pch=20, cex=2)
drawIsopycnals(levels=rho0, col="red", cex=0)
segments(Sa, Ta, Sb, Tb, col="blue")
text(Sb, Tb, "b", pos=4); text(Sa, Ta, "a", pos=4); text(Sc, Tc, "c", pos=4)
legend("topleft",legend=sprintf("Sa=%.1f, Ta=%.1f, Sb=%.1f -> Tb=%.1f, drho=%.2f, dT=%.2f",Sa, Ta, Sb, Tb, drho, dT),bg="white")
In [44]:
HTML(css)
Out[44]: