ACINN workshop, Tue 07.02.2017
Fabien Maussion
Slides: http://fabienmaussion.info/acinn_xarray_workshop
Notebook: On GitHub
# Ignore numpy warnings
import warnings
warnings.filterwarnings('ignore')
import matplotlib.pyplot as plt
%matplotlib inline
import xarray as xr
# Some defaults:
plt.rcParams['figure.figsize'] = (12, 6) # Default plot size|
xr.set_options(display_width=64); # same here
Documentation: http://xarray.pydata.org
Repository: https://github.com/pydata/xarray
Initial release: 03.05.2014
Latest release: v0.9.1 (20.01.2017)
53 contributors (latest release: 24)
Umbrellas: Python for data & NumFOCUS (but no funding...)
import numpy as np
a = np.array([[1, 3, 9], [2, 8, 4]])
a
a[1, 2]
a.mean(axis=0)
import xarray as xr
da = xr.DataArray(a, dims=['lat', 'lon'],
coords={'lon':[11, 12, 13], 'lat':[1, 2]})
da
da.sel(lon=13, lat=2).values
da.mean(dim='lat')
f = 'ERA-Int-MonthlyAvg-4D-TUVWZ.nc'
ds = xr.open_dataset(f)
ds
ds.t.sel(month=8, level=850)
ds.t.isel(month=7, level=11)
ds.t.sel(level=1001, latitude=47.26, longitude=11.38, method='nearest')
ds.t[7, 11, :, :]
ds.u.mean(dim=['month', 'longitude']).plot.contourf(levels=13)
plt.ylim([1000, 100]);
u_avg = ds.u.mean(dim=['month', 'longitude'])
u_avg_masked = u_avg.where(u_avg > 12)
u_avg_masked.plot.contourf(levels=13)
plt.ylim([1000, 100]);
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
a = xr.DataArray(np.arange(3), dims='time',
coords={'time':np.arange(3)})
b = xr.DataArray(np.arange(5), dims='time',
coords={'time':np.arange(5)+1})
a + b
ts = ds.t.sel(level=1001, latitude=47.26, longitude=11.38, method='nearest')
ts.plot();
import cartopy.crs as ccrs
ax = plt.axes(projection=ccrs.Robinson())
ds.z.sel(level=1000, month=8).plot(ax=ax, transform=ccrs.PlateCarree());
ax.coastlines();
Opening all files in a directory...
mfs = '/home/mowglie/disk/Data/Gridded/GPM/3BDAY_sorted/*.nc'
dsmf = xr.open_mfdataset(mfs)
... results in a consolidated dataset ...
dsmf
... on which all usual operations can be applied:
dsmf = dsmf.sel(time='2015')
dsmf
Yes, even computations!
ts = dsmf.precipitationCal.mean(dim=['lon', 'lat'])
ts
Computations are done "lazily"
No actual computation has happened yet:
ts.data
But they can be triggered:
ts = ts.load()
ts
For more information: http://xarray.pydata.org/en/stable/dask.html
ts.plot();
ts.rolling(time=31, center=True).mean().plot();
Taken from: http://ajdawson.github.io/eofs/examples/nao_xarray.html
from eofs.xarray import Eof
from eofs.examples import example_data_path
# Read geopotential height data using the xarray module
filename = example_data_path('hgt_djf.nc')
z_djf = xr.open_dataset(filename)['z']
# Compute anomalies by removing the time-mean.
z_djf = z_djf - z_djf.mean(dim='time')
# Create an EOF solver to do the EOF analysis.
coslat = np.cos(np.deg2rad(z_djf.coords['latitude'].values)).clip(0., 1.)
solver = Eof(z_djf, weights=np.sqrt(coslat)[..., np.newaxis])
# Get the leading EOF
eof1 = solver.eofsAsCovariance(neofs=1)
# Leading EOF expressed as covariance in the European/Atlantic domain
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-20, central_latitude=60))
ax.coastlines() ; ax.set_global()
eof1[0, 0].plot.contourf(ax=ax, levels=np.linspace(-75, 75, 11),
cmap=plt.cm.RdBu_r, add_colorbar=False,
transform=ccrs.PlateCarree())
ax.set_title('EOF1 expressed as covariance', fontsize=16);
http://salem.readthedocs.io/en/latest/
Try it out:
pip install salem
# importing salem adds a new "toolbox" to xarray objects
import salem
pday = dsmf.precipitationCal.sel(time='2015-02-01')
cm = pday.salem.quick_map(cmap='Blues', vmax=100);
shdf = salem.read_shapefile(salem.get_demo_file('world_borders.shp'))
shdf = shdf.loc[shdf['CNTRY_NAME'].isin(['Peru'])]
dsmfperu = dsmf.salem.subset(shape=shdf, margin=10)
pday = dsmfperu.precipitationCal.sel(time='2015-02-01')
cm = pday.salem.quick_map(cmap='Blues', vmax=100);
dsmfperu = dsmfperu.salem.roi(shape=shdf)
pday = dsmfperu.precipitationCal.sel(time='2015-02-01')
cm = pday.salem.quick_map(cmap='Blues', vmax=100);
prpc_a = dsmfperu.precipitationCal.sum(dim=['time']).load()
prpc_a.salem.quick_map(cmap='Blues', vmax=5000);
Problems:
f = 'wrfpost_d01_2005-09-21_00-00-00_25h.nc'
ds = xr.open_dataset(f)
ds
wrf = salem.open_wrf_dataset(f)
wrf
wrf.T2C.mean(dim='time', keep_attrs=True).salem.quick_map();
ws_h = wrf.isel(time=5).salem.wrf_zlevel('WS', levels=10000.)
ws_h.salem.quick_map(cmap='Reds');