Prerequisites¶
To run this notebook your python environment needs the following packages:
- numpy
- xarray
- cftime
- matplotlib
- cartopy
- netcdf4
A suitable environment can be installed using anaconda or mamba via:
mamba create -n py_env_xarray numpy cftime=1.2.1 jupyter xarray matplotlib cartopy netcdf4
source activate py_env_xarray
In [1]:
# importing packages
import numpy as np
import xarray as xr
import cftime
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
In [7]:
# change matplotlib standard parameters
plt.rcParams['figure.figsize'] = 12,8
Plot Data¶
Loading data from disk. One PISM output file for the Antarctic and one for the Greenland ice sheet.
In [3]:
ds_ais = xr.open_dataset('data/ais_extra.nc')
ds_gris = xr.open_dataset('data/gris_extra.nc')
xarray's standard plotting functionality¶
xarray automatically uses available metadata to label axes, adds a colorbar, etc.
In [4]:
fig, axis = plt.subplots(1,1, figsize=(8,6))
ds_ais['velsurf_mag'].plot(ax=axis)
axis.set_aspect('equal')
Plot data on a projection with cartopy¶
The cartopy package allows to plot data on maps and projections, which is useful when visualising spheric coordinates for example. For the Antarctic dataset, we choose a stereographic southpole projection and add a progressive colorbar.
In [5]:
map_proj = ccrs.SouthPolarStereo()
levels = [0, 1, 10, 100, 1000, 2000, 4000]
ax = ds_ais['velsurf_mag'].plot(
transform=ccrs.PlateCarree(),
x='lon', y='lat',
cmap='cividis',
levels=levels,
cbar_kwargs={'ticks': levels},
subplot_kws={'projection': map_proj})
ax.axes.gridlines(draw_labels=True)
Out[5]:
<cartopy.mpl.gridliner.Gridliner at 0x7f968466d730>
Adding coastlines with cartopy¶
Cartopy also comes with a predefined set of coastlines which can sometimes be useful to refer the data geographically. We plot here the ice sheet thickness of the Greenland data set.
In [6]:
map_proj = ccrs.Orthographic(central_longitude=-45.0, central_latitude=80.0)
levels = [0, 10, 100, 1000, 1500, 2000, 2500, 3000]
ax = ds_gris['thk'].where(ds_gris['thk']>0,np.nan).plot(
transform=ccrs.PlateCarree(),
x='lon', y='lat',
levels=levels,
cbar_kwargs={'ticks': levels},
subplot_kws={'projection': map_proj}
)
ax.axes.coastlines()
gl = ax.axes.gridlines(draw_labels=True)
gl.bottom_labels = False
gl.right_labels = False
In [ ]: