It has been a while that I wanted to play with Jake's JSAnimation module. (You can see his example here.)
My animation will be a 1-year sequence of Sea Surface Height Anomaly using
AVISO data. For that I'll use iris to
download the data via OpenDAP, and cartopy to plot it.
import iris
url = 'http://%s:%s@opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-dt-upd-global-merged-msla-h' % (user, passwd)
longitude = lambda cell: 300 <= cell <= 360
latitude = lambda cell: -60 <= cell <= 0
constraint = iris.Constraint(coord_values=dict(longitude=longitude,
latitude=latitude), name='H')
cube = iris.load_cube(url, constraint)
print(cube)
To save download time I stored a local copy of the data (just the last 52 weeks, that is 1 year of data). Note also that in this block we use 3 cool iris features:
- 1-line data conversion (from cm $\rightarrow$ meters).
- 1-line CF-compliant netcdf file!
- The data was not downloaded in the block above. The download only happens when we try to access the data like in the unit conversion step.
import os
if os.path.isfile('./data/altimetry.nc'):
cube = iris.load_cube('./data/altimetry.nc')
else:
cube = cube[-52:, ...]
cube.convert_units('meters')
iris.save(cube, 'altimetry.nc')
time = cube.dim_coords[0]
lon = cube.dim_coords[1].points
lat = cube.dim_coords[2].points
time = time.units.num2date(time.points)
print(cube)
Now let's plot just the first time slice to check the data.
import iris.plot as iplt
import cartopy.crs as ccrs
import iris.quickplot as qplt
import matplotlib.pyplot as plt
from cartopy.feature import NaturalEarthFeature, LAND
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
k = 0
fig, ax = plt.subplots(figsize=(4, 4), subplot_kw=dict(projection=ccrs.PlateCarree()))
qplt.contourf(cube[k, ...], cmap=plt.cm.GnBu, extend='both')
ax.add_feature(LAND)
ax.coastlines('50m', color='k')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
_ = iplt.citation(time[k].strftime('%Y-%m-%d'))
Before jumping into the animation we need to create a colormap with the white color centered at zero values to visualize better the anomalies. Luckily Phil Elson (the guy behind cartopy and iris) already did that. We will just change the usual Red/White/Blue to Green/White/Yellow (Yep, World Cup fever!!!).
from matplotlib.colors import LinearSegmentedColormap
levels = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 0.9] # Must be even so zero is white.
colors = [(1, 1, 0), (1., 1., 1.), (0, 1, 0)]
cmap = LinearSegmentedColormap.from_list(colors=colors,
N=len(levels)-1,
name='brasil_cmap',)
This block imports the JSAnimation and FuncAnimation modules, and defines the
update figure function in the matplotlib animation way.
from JSAnimation import IPython_display
from matplotlib.animation import FuncAnimation
def draw_coastline(ax):
ax.coastlines('50m', color='k')
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1.5, color='gray', alpha=0.5, linestyle='--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
ax.add_feature(LAND)
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
name='admin_1_states_provinces_shp')
ax.add_feature(states, edgecolor='gray')
def update_contourf(k):
global ax, fig
ax.cla()
cs = ax.contourf(lon, lat, cube[k, ...].data.T, **kw)
text = ax.text(-7.5, -7.5, time[k].strftime('%Y-%m-%d'),
ha='center', va='center')
draw_coastline(ax)
return cs
And now the final result!
kw = dict(cmap=cmap, clim=(-1, 1), levels=levels)
fig, ax = plt.subplots(figsize=(6, 6), subplot_kw=dict(projection=ccrs.PlateCarree()))
fig.patch.set_visible(False)
ax.axis([-59, 0, -59, 0])
cs = ax.contourf(lon, lat, cube[0, ...].data.T, **kw)
fig.colorbar(cs, shrink=0.7)
text = ax.text(-7.5, -7.5, time[0].strftime('%Y-%m-%d'), ha='center', va='center')
draw_coastline(ax)
FuncAnimation(fig, update_contourf, interval=100, frames=52)
This last part is just to show how slow these anomalies propagate in our domain.
dist_degree = 15
degree2km = 111.
time_weeks = 43
week2hour = 7 * 24
total_time = (dist_degree * degree2km) / (time_weeks * week2hour)
print("Velocity of the anomaly while crossing 15 longitude degrees is %0.2f km/h" % total_time)
time_years = (time_weeks * 7 * 3) / 365.
print(u"Time it takes an anomaly from 0 degrees of longitude to reach\n"
u"São Paulo at 45 degrees of longitude is %1.2f years" % time_years)
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.