In two previous posts (here and here) I showed how to download and plot World Ocean Atlas Data. This post is an update where we will be using the latest oceans library and the lasted Atlas database (WOA13).

First let's define a mapping function using cartopy.

In [2]:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import (LONGITUDE_FORMATTER,
                                   LATITUDE_FORMATTER)


LAND = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                    edgecolor='face',
                                    facecolor=cfeature.COLORS['land'])


def make_map(bbox, projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(8, 6),
                           subplot_kw=dict(projection=projection))
    ax.set_extent(bbox)
    ax.add_feature(LAND, facecolor='0.75')
    ax.coastlines(resolution='50m')
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax

Now we can import the new woa_subset from oceans. This version returns an iris cube instance instead of a Pandas Panel. The call below requests annual temperature climatology at 0.25 degrees resolution for the whole globe. The latest keyword argument requests the whole dataset and not only the objectively analyzed temperature (default is False).

In [3]:

import matplotlib.pyplot as plt
from oceans.colormaps import cm
from oceans.ff_tools import wrap_lon180

from oceans.datasets import woa_subset


bbox = [2.5, 357.5, -87.5, 87.5]
kw = dict(bbox=bbox, variable='temperature', clim_type='00',
          resolution='0.25', full=True)

cubes = woa_subset(**kw)
print(cubes)

0: sea_water_temperature / (degrees_celsius) (time: 1; depth: 102; latitude: 702; longitude: 1422)
1: sea_water_temperature / (1)         (time: 1; depth: 102; latitude: 702; longitude: 1422)
2: sea_water_temperature / (degrees_celsius) (time: 1; depth: 102; latitude: 702; longitude: 1422)
3: sea_water_temperature / (degrees_celsius) (time: 1; depth: 102; latitude: 702; longitude: 1422)
4: sea_water_temperature / (degrees_celsius) (time: 1; depth: 102; latitude: 702; longitude: 1422)
5: sea_water_temperature / (degrees_celsius) (time: 1; depth: 102; latitude: 702; longitude: 1422)
6: sea_water_temperature / (1)         (time: 1; depth: 102; latitude: 702; longitude: 1422)

/home/filipe/.virtualenvs/blog/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1619: UserWarning: NetCDF variable 't_sd' contains unknown cell method 'standard'
  warnings.warn('NetCDF variable {!r} contains unknown cell method {!r}'.format('{}'.format(cf_var_name), '{}'.format(d[CM_METHOD])))

The standard_name is confusing. All the variables are sea_water_temperature. To understand the dataset we need to inspect the long_name or the attribute description.

In [4]:

[(k, c.var_name, c.long_name) for k, c in enumerate(cubes)]

Out[4]:

[(0,
  u't_mn',
  u'Average of all unflagged interpolated values at each standard depth level for sea_water_temperature in each grid-square which contain at least one measurement.'),
 (1,
  u't_dd',
  u'The number of observations of sea_water_temperature in each grid-square at each standard depth level.'),
 (2,
  u't_sd',
  u'The standard deviation about the statistical mean of sea_water_temperature in each grid-square at each standard depth level.'),
 (3,
  u't_se',
  u'The standard error about the statistical mean of sea_water_temperature in each grid-square at each standard depth level.'),
 (4,
  u't_an',
  u'Objectively analyzed mean fields for sea_water_temperature at standard depth levels.'),
 (5,
  u't_oa',
  u'statistical mean value minus the objectively analyzed mean value for sea_water_temperature.'),
 (6,
  u't_gp',
  u'The number of grid-squares within the smallest radius of influence around each grid-square which contain a statistical mean for sea_water_temperature.')]

We want the "Objectively analyzed mean fields for sea_water_temperature at standard depth levels" or var_name==t_an. (Note that this is the only variable returned when calling with full=False.)

Let's plot the surface data.

In [5]:

import matplotlib.pyplot as plt
import iris.plot as iplt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from palettable import colorbrewer



def make_map(cube, projection=ccrs.PlateCarree(), figsize=(12, 10),
             cmap=cm.avhrr):
    fig, ax = plt.subplots(figsize=figsize,
                           subplot_kw=dict(projection=projection))
    ax.add_feature(cfeature.LAND, facecolor='0.75')
    cs = iplt.pcolormesh(cube, cmap=cmap)
    ax.coastlines()
    if isinstance(projection, ccrs.PlateCarree) or isinstance(projection, ccrs.Mercator):
        gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1.5,
                          color='gray', alpha=0.5, linestyle='--')
        gl.xlabels_top = gl.ylabels_right = False
        gl.xformatter = LONGITUDE_FORMATTER
        gl.yformatter = LATITUDE_FORMATTER
    cbar = dict(extend='both', shrink=0.5, pad=0.02,
                orientation='horizontal', fraction=0.1)
    cb = fig.colorbar(cs, **cbar)
    cb.ax.set_xlabel(r"$^{\circ}$ C")
    return fig, ax

Different from the previous woa_subset version, this one returns all depth levels and not the requested ones. The iris cube is "lazy", so there is no reason to slice it earlier because no data is downloaded until you actually need to.

In [6]:

cube = cubes[4]
c = cube[0, 0, ...]  # Slice singleton time and first level.

fig, ax = make_map(c, projection=ccrs.Robinson())

I also added a woa_profile version to quickly extract vertical profiles at a given position. Here is how to use it:

In [7]:

from oceans import woa_profile

kw = dict(variable='temperature', clim_type='00', resolution='1.00', full=False)
polar = woa_profile(-25.5, -70.5, **kw)
tempe = woa_profile(-25.5, -45.0, **kw)
equat = woa_profile(-25.5, 0, **kw)

In [8]:

from matplotlib import style
import matplotlib.pyplot as plt


style.use('ggplot')

def plot_profile(cube, label):
    z = cube.coord(axis='Z').points
    l = ax.plot(cube[0, :].data, z, label=label, linewidth=2)

fig, ax = plt.subplots(figsize=(2.25, 6))
plot_profile(polar, label='Polar')
plot_profile(tempe, label='Temperate')
plot_profile(equat, label='Equatorial')
ax.legend(bbox_to_anchor=(1.5, .95))
ax.invert_yaxis()
_ = ax.set_ylim(4000, 0)

Happy May the fourth!

In [9]:

from IPython.display import YouTubeVideo

YouTubeVideo("ngElkyQ6Rhs")

Out[9]:

In [10]:

HTML(html)

Out[10]:

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/.

python4oceanographers

Turning ripples into waves

Downloading World Ocean Atlas 2013 with oceans and iris

Comments