python4oceanographers

Turning ripples into waves

Using JSAnimation with Sea Surface Height data

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.

In [2]:
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)
H / (cm)                            (time: 1087; longitude: 179; latitude: 227)
     Dimension coordinates:
          time                           x                -              -
          longitude                      -                x              -
          latitude                       -                -              x
     Attributes:
          CreatedBy: SSALTO/DUACS
          CreatedOn: 03-DEC-2013 09:23:10:000000
          Date_CNES_JD: 23222.0
          FileType: GRID_DOTS_MERCATOR
          OriginalName: dt_upd_global_merged_msla_h_20130731_20130731_20131203.nc
          date: 2013-07-31 00:00:00.000000 UTC
          title: SSALTO/DUACS - DT MSLA - Merged Product - Homogeneous Global Processin...

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.
In [3]:
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)
H / (meters)                        (time: 52; longitude: 179; latitude: 227)
     Dimension coordinates:
          time                           x              -              -
          longitude                      -              x              -
          latitude                       -              -              x
     Attributes:
          Conventions: CF-1.5
          CreatedBy: SSALTO/DUACS
          CreatedOn: 02-OCT-2013 08:30:38:000000
          Date_CNES_JD: 23110.0
          FileType: GRID_DOTS_MERCATOR
          OriginalName: dt_upd_global_merged_msla_h_20130410_20130410_20131002.nc
          date: 2013-04-10 00:00:00.000000 UTC
          title: SSALTO/DUACS - DT MSLA - Merged Product - Homogeneous Global Processin...

Now let's plot just the first time slice to check the data.

In [4]:
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!!!).

In [5]:
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.

In [6]:
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!

In [7]:
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)
Out[7]:


Once Loop Reflect