Plotting Tiles¶

Objectives¶

Introduce several different methods for plotting ECCO v4 fields that are stored as tiles in Datasets or DataArrays. Emphasis is placed on fields stored on the ECCO v4 native llc90 grid and loaded from NetCDF tile files.

Introduction¶

“Over the years many different plotting modules and packages have been developed for Python. For most of that time there was no clear favorite package, but recently matplotlib has become the most widely used. Nevertheless, many of the others are still available and may suit your tastes or needs better. Some of these are interfaces to existing plotting libraries while others are Python-centered new implementations. – from : https://wiki.python.org/moin/NumericAndScientific/Plotting

The link above profiles a long list of Python tools for plotting. In this tutorial we use just two libraries, matplotlib and Cartopy.

Note: In this tutorial you will need to have downloaded monthly SSH, THETA, and SALT for the year 2000. The ShortNames of the datasets needed are ECCO_L4_SSH_LLC0090GRID_MONTHLY_V4R4 and ECCO_L4_TEMP_SALINITY_LLC0090GRID_MONTHLY_V4R4. You will also need the grid file.

matplotlib¶

“Matplotlib is a Python 2D plotting library which produces publication quality figures in a variety of hardcopy formats and interactive environments across platforms. Matplotlib can be used in Python scripts, the Python and IPython shell, the jupyter notebook, web application servers, and four graphical user interface toolkits.”

“For simple plotting the pyplot module provides a MATLAB-like interface, particularly when combined with [Juypter Notebooks]. For the power user, you have full control of line styles, font properties, axes properties, etc, via an object oriented interface or via a set of functions familiar to MATLAB users.” – from https://matplotlib.org/index.html

Matplotlib and pyplot even have a tutorial: https://matplotlib.org/users/pyplot_tutorial.html

Cartopy¶

“Cartopy is a Python package designed for geospatial data processing in order to produce maps and other geospatial data analyses.”

Cartopy makes use of the powerful PROJ.4, NumPy and Shapely libraries and includes a programmatic interface built on top of Matplotlib for the creation of publication quality maps.

Key features of cartopy are its object oriented projection definitions, and its ability to transform points, lines, vectors, polygons and images between those projections.

You will find cartopy especially useful for large area / small scale data, where Cartesian assumptions of spherical data traditionally break down. If you’ve ever experienced a singularity at the pole or a cut-off at the dateline, it is likely you will appreciate cartopy’s unique features!”*

– from https://scitools.org.uk/cartopy/docs/latest/

The default orientation of the lat-lon-cap tile fields¶

Before we begin plotting ECCOv4 fields on the native llc90 model grid we are reminded how how the 13 tiles are oriented with respect to their “local” x and y and with respect to each other.

llc90 tile layout

Tiles 7-12 are rotated 90 degrees counter-clockwise relative to tiles 0-5.

Note: The rotated orientation of tiles 7-12 presents some complications but don’t panic! The good news is that you don’t need to reorient tiles to plot them.

Plotting single tiles using `imshow`, `pcolormesh`, and `contourf`¶

First, let’s load the all 13 tiles for sea surface height and the model grid parameters.

[1]:

import numpy as np
import sys
import xarray as xr
import matplotlib.pyplot as plt
%matplotlib inline
import glob
import warnings
warnings.filterwarnings('ignore')

[2]:

# load some useful cartopy routines
from cartopy import config
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# and a new matplotlib routine
import matplotlib.path as mpath

[3]:

## Import the ecco_v4_py library into Python
## =========================================
##    If ecco_v4_py is not installed in your local Python library,
##    tell Python where to find it.  The example below adds
##    ecco_v4_py to the user's path if it is stored in the folder
##    ECCOv4-py under the user's home directory

from os.path import join,expanduser
user_home_dir = expanduser('~')

sys.path.append(join(user_home_dir,'ECCOv4-py'))

import ecco_v4_py as ecco

[4]:

## Set top-level file directory for the ECCO NetCDF files
## =================================================================

## currently set to ~/Downloads/ECCO_V4r4_PODAAC,
## the default if ecco_podaac_download was used to download dataset granules
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')

[5]:

## Load the model grid
ecco_grid = xr.open_dataset(glob.glob(join(ECCO_dir,'*GEOMETRY*','*.nc'))[0])

## Load one year of 2D monthly data: SSH, temperature and salinity
ds_SSH = xr.open_mfdataset(join(ECCO_dir,'*SSH*MONTHLY*','*_2000-*.nc'))
ds_temp_sal = xr.open_mfdataset(join(ECCO_dir,'*TEMP*SAL*MONTHLY*','*_2000-*.nc'))
## select only *surface* temperature and salinity (SST and SSS)
ds_SST_SSS = ds_temp_sal.isel(k=0)


## Copy ecco_ds from ecco_grid dataset
ecco_ds = ecco_grid.copy()
## Add SSH, SST, and SSS variables to ecco_ds
ecco_ds['SSH'] = ds_SSH['SSH']
ecco_ds['SST'] = ds_SST_SSS['THETA']
ecco_ds['SSS'] = ds_SST_SSS['SALT']

## Load ecco_ds into memory
ecco_ds = ecco_ds.load()

[6]:

ecco_ds

Plotting a single tile with `imshow`¶

First we’ll plot the average SSH for the first month (Jan 2000) on tiles 3, 7, and 8 using the basic imshow routine from pyplot. We are plotting these three different tiles to show that these lat-lon-cap tiles all have a different orientation in $x$ and $y$ .

Note: Theorigin=’lower’argument to ``imshow`` is required to make the :math:`y` origin at the bottom of the plot.

Tile 2 (Northeast Atlantic)¶

[7]:

plt.figure(figsize=(6,5), dpi= 90)

# Step 1, select the tile to plot using the **.isel( )** syntax.
tile_to_plot = ecco_ds.SSH.isel(tile=2, time=0)
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=2,k=0) !=0, np.nan)

# Step 2, use plt.imshow()
plt.imshow(tile_to_plot, origin='lower');

# Step 3, add colorbar, title, and x and y axis labels
plt.colorbar()
plt.title('SSH (m) ' + str(ecco_ds.time[0].values)[0:7] + '\n tile 2')
plt.xlabel('x -->')
plt.ylabel('y -->')

[7]:

Text(0, 0.5, 'y -->')

Tiles 0-5 are by default in a quasi-lat-lon orientation. + $x$ is to the east and + $y$ is to the north.

Tile 6 (the Arctic cap)¶

This time we’ll plot the Arctic cap tile 6. Notice the layout of the Arctic cap tile in $x$ and $y$ . We’ll follow the same procedure for plotting except we’ll use LaTeX to add arrows in the $x$ and $y$ axis labels (for fun).

[8]:

plt.figure(figsize=(6,5), dpi= 90)

# Step 1, select the tile to plot using the **.isel( )** syntax.
tile_to_plot = ecco_ds.SSH.isel(tile=6, time=0)
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=6,k=0) !=0, np.nan)

# Step 2, use plt.imshow()
plt.imshow(tile_to_plot, origin='lower');

# Step 3, add colorbar, title, and x and y axis labels
plt.colorbar()
plt.title('SSH (m) ' + str(ecco_ds.time[0].values)[0:7] + '\n tile 6')
plt.xlabel('x -->');
plt.ylabel('y -->');

Because tile 6 is the Arctic cap, $x$ and $y$ do not map to east and west throughout the domain.

Tile 7 (N. Pacific / Bering Sea / Chukchi Sea)¶

For tiles 7-12 , positive $x$ is southwards and positive $y$ is eastwards.

[9]:

plt.figure(figsize=(6,5), dpi= 90)

# pull out lats and lons
tile_num=7
tile_to_plot = ecco_ds.SSH.isel(tile=tile_num, time=1)
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=tile_num,k=0) !=0, np.nan)

plt.imshow(tile_to_plot)
plt.colorbar()
plt.title('SSH (m) ' + str(ecco_ds.time[1].values)[0:7] + '\n tile ' + str(tile_num))
plt.xlabel('longitude');
plt.ylabel('latitude');

Tiles 7-12 are are also in a quasi-lat-lon orientation except that + $x$ is roughly south and + $y$ is roughly east.

Plotting a single tile with `pcolor` and `contourf`¶

The pcolor and contourf routines allows us to add latitude and longitude to the figure. Because SSH is a ‘c’ point variable, its lat/lon coordinates are YC and XC

We can’t plot the Arctic cap tile with pcolor and contourf using latitude and longitude for the plot x and y axes because of the singularity at the pole and the 360 wrapping in longitude.

Instead, we will demonstrate pcolor and contourf for tile 2.

Tile 2 (Northeast N. Atlantic)¶

[10]:

fig=plt.figure(figsize=(10, 10))

tile_num=2
time_ind=1

# pull out lats and lons
lons = ecco_ds.XC.sel(tile=tile_num)
lats = ecco_ds.YC.sel(tile=tile_num)
tile_to_plot = ecco_ds.SSH.isel(tile=tile_num, time=time_ind)

# mask to NaN where hFacC is == 0
# syntax is actually "keep where hFacC is not equal to zero"
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=tile_num,k=0) !=0, np.nan)

# create subplot for pcolor
fig = plt.subplot(221)

# use pcolor with 'lons' and 'lats' for the plot x and y axes
plt.pcolor(lons, lats, tile_to_plot, vmin=-1, vmax=1, cmap='RdBu_r')
plt.colorbar()
plt.title('PCOLOR SSH (m) \n' + str(ecco_ds.time[time_ind].values)[0:7] + '\n tile ' + str(tile_num))
plt.xlabel('longitude')
plt.ylabel('latitude')

# create subplot for contourf
fig=plt.subplot(222)

# use contourf with 'lons' and 'lats' for the plot x and y axes
plt.contourf(lons, lats, tile_to_plot, np.linspace(-1,1, 20,endpoint=True), cmap='RdBu_r', vmin=-1, vmax=1)
plt.title('CONTOURF SSH (m) \n' + str(ecco_ds.time[time_ind].values)[0:7] + '\n tile ' + str(tile_num))
plt.xlabel('longitude')
plt.ylabel('latitude')
plt.colorbar()

# plot every 3rd model grid line to show how tile 3 is 'warped' above around 60N
plt.plot(ecco_ds.XG.isel(tile=tile_num)[::3,::3], ecco_ds.YG.isel(tile=tile_num)[::3,::3],'k-')
plt.plot(ecco_ds.XG.isel(tile=tile_num)[::3,::3].T, ecco_ds.YG.isel(tile=tile_num)[::3,::3].T,'k-')

# push the subplots away from each other a bit
plt.subplots_adjust(bottom=0, right=1.2, top=.9)

Tile 7 (N. Pacific / Bering Sea / Chukchi Sea)¶

If longitude and latitude are passed as the ‘x’ and ‘y’ arguments to pcolor and contourf then the fields will be oriented geographically.

[11]:

fig=plt.figure(figsize=(10, 10))

tile_num=7
time_ind=1

# pull out lats and lons
lons = np.copy(ecco_ds.XC.sel(tile=tile_num))

# we must convert the longitude coordinates from
# [-180 to 180] to [0 to 360]
# because of the crossing of the international date line.
lons[lons < 0] = lons[lons < 0]+360

lats = ecco_ds.YC.sel(tile=tile_num)
tile_to_plot = ecco_ds.SSH.isel(tile=tile_num, time=time_ind)

# mask to NaN where hFacC is == 0
# syntax is actually "keep where hFacC is not equal to zero"
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=tile_num,k=0) !=0, np.nan)

# create subplot for pcolor
fig = plt.subplot(221)

# use pcolor with 'lons' and 'lats' for the plot x and y axes
plt.pcolor(lons, lats, tile_to_plot, vmin=-1, vmax=1.1, cmap='RdBu_r')
plt.colorbar()
plt.title('PCOLOR SSH (m) \n' + str(ecco_ds.time[time_ind].values)[0:7] + '\n tile ' + str(tile_num))
plt.xlabel('longitude')
plt.ylabel('latitude')

# create subplot for contourf
fig=plt.subplot(222)

# use contourf with 'lons' and 'lats' for the plot x and y axes
plt.contourf(lons, lats, tile_to_plot, np.linspace(-1,1.1,22,endpoint=True), cmap='RdBu_r', vmin=-1, vmax=1.1)
plt.title('CONTOURF SSH (m) \n' + str(ecco_ds.time[time_ind].values)[0:7] + '\n tile ' + str(tile_num))
plt.xlabel('longitude')
plt.ylabel('latitude')

plt.colorbar()

# push the subplots away from each other a bit
plt.subplots_adjust(bottom=0, right=1.2, top=.9)

Plotting fields from one tile using `Cartopy`¶

The Cartopy package provides routines to make plots using different geographic projections. We’ll demonstrate plotting these three tiles again using Cartopy.

To see a list of Cartopy projections, see http://pelson.github.io/cartopy/crs/projections.html

Geographic Projections (AKA: plate carrée)¶

Cartopy works by transforming geographic coordintes (lat/lon) to new x,y coordinates associated with different projections. The most familiar projection is the so-called geographic projection (aka plate carree). When we plotted tiles using pcolor and contourf we were de-factor using the plate carree projection longitude and latitude were the ‘x’ and the ‘y’ of the plot.

With Cartopy we can make similar plots in the plate carree projection system and also apply some cool extra details, like land masks.

We’ll demonstrate on tiles 2 and 7 (again skipping tile 6 (Arctic cap) because we cannot use geographic coordinates as x and y when there is a polar singularity and 360 degrees of longitude.

Tile 2 with plate carree¶

[12]:

tile_num=2
time_ind=1

lons = ecco_ds.XC.isel(tile=tile_num)
lats = ecco_ds.YC.isel(tile=tile_num)

tile_to_plot = ecco_ds.SSH.isel(tile=tile_num, time=time_ind)
# mask to NaN where hFacC is == 0
# syntax is actually "keep where hFacC is not equal to zero"
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=tile_num,k=0) !=0, np.nan)

fig = plt.figure(figsize=(10,5), dpi= 90)

# here is where you specify what projection you want to use
ax = plt.axes(projection=ccrs.PlateCarree())

# here is here you tell Cartopy that the projection
# of your 'x' and 'y' are geographic (lons and lats)
# and that you want to transform those lats and lons
# into 'x' and 'y' in the projection
cf = plt.contourf(lons, lats, tile_to_plot, 60,
                  transform=ccrs.PlateCarree());

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--');
ax.coastlines()
ax.add_feature(cfeature.LAND)

# add separate axes for colorbar (to ensure it doesn't overlap main plot)
cbar_ax = fig.add_axes([0.88,0.1,0.02,0.8])
plt.colorbar(cf,ax=ax,cax=cbar_ax)

[12]:

<matplotlib.colorbar.Colorbar at 0x2926ca8f4f0>

Other features we could have added include:

cartopy.feature.BORDERS
    Country boundaries.

cartopy.feature.COASTLINE
    Coastline, including major islands.

cartopy.feature.LAKES
    Natural and artificial lakes.

cartopy.feature.LAND
    Land polygons, including major islands.

cartopy.feature.OCEAN
    Ocean polygons.

cartopy.feature.RIVERS
    Single-line drainages, including lake centerlines.

Let’s add geographic borders just to demonstrate how extra features can be added to a Cartopy map

[13]:

fig = plt.figure(figsize=(10,5), dpi= 90)

# here is where you specify what projection you want to use
ax = plt.axes(projection=ccrs.PlateCarree())

# here is here you tell Cartopy that the projection
# of your 'x' and 'y' are geographic (lons and lats)
# and that you want to transform those lats and lons
# into 'x' and 'y' in the projection
cf = plt.contourf(lons, lats, tile_to_plot, 60,
                  transform=ccrs.PlateCarree());

gl = ax.gridlines(crs=ccrs.PlateCarree(), \
                  draw_labels=True,
                  linewidth=2, color='gray', \
                  alpha=0.5, linestyle='--');
ax.coastlines()
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.BORDERS)

# add separate axes for colorbar (to ensure it doesn't overlap main plot)
cbar_ax = fig.add_axes([0.88,0.1,0.02,0.8])
plt.colorbar(cf,ax=ax,cax=cbar_ax)

[13]:

<matplotlib.colorbar.Colorbar at 0x2926cb8c940>

Tile 7 with plate carree¶

To use the plate carree projection across the international date line specify the central_longitude=-180 argument when defining the projection and for creating the gridlines. (see https://stackoverflow.com/questions/13856123/setting-up-a-map-which-crosses-the-dateline-in-cartopy)

[14]:

tile_num=7
time_ind=1

# pull out lats and lons
lons = np.copy(ecco_ds.XC.sel(tile=tile_num))

# we must convert the longitude coordinates from
# [-180 to 180] to [0 to 360]
# because of the crossing of the international date line.
lons[lons < 0] = lons[lons < 0]+360
lats = ecco_ds.YC.isel(tile=tile_num)

tile_to_plot = ecco_ds.SSH.isel(tile=tile_num, time=time_ind)
# mask to NaN where hFacC is == 0
# syntax is actually "keep where hFacC is not equal to zero"
tile_to_plot= tile_to_plot.where(ecco_ds.hFacC.isel(tile=tile_num,k=0) !=0, np.nan)

fig = plt.figure(figsize=(10,5), dpi= 90)

# here is where you specify what projection you want to use
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=-180))

# here is here you tell Cartopy that the projection of your 'x' and 'y' are geographic (lons and lats)
# and that you want to transform those lats and lons into 'x' and 'y' in the projection
cf = plt.contourf(lons, lats, tile_to_plot, 60,
                  transform=ccrs.PlateCarree())

gl = ax.gridlines(crs=ccrs.PlateCarree(central_longitude=-180), draw_labels=True,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--');
ax.coastlines()
ax.add_feature(cfeature.LAND)

# add separate axes for colorbar (to ensure it doesn't overlap main plot)
cbar_ax = fig.add_axes([0.88,0.1,0.02,0.8])
plt.colorbar(cf,ax=ax,cax=cbar_ax)

[14]:

<matplotlib.colorbar.Colorbar at 0x2925c6d9550>

Polar Stereographic Projection¶

[15]:

tile_num=6
time_ind=1

# use lower-left grid corner coordinates for flat pcolormesh shading --
# this avoids some funky effects when plot is transformed into
# the polar stereographic projection
lon_corners = ecco_ds.XG.isel(tile=tile_num)
lat_corners = ecco_ds.YG.isel(tile=tile_num)


tile_to_plot = ecco_ds.SSH.isel(tile=tile_num, time=time_ind)

# mask to NaN where hFacC is == 0
# syntax is actually "keep where hFacC is not equal to zero"
tile_to_plot = tile_to_plot.where(ecco_ds.hFacC.isel(tile=tile_num,k=0) !=0, \
                       np.nan)

plt.figure(figsize=(8,6), dpi= 90)


# Make a new projection, time of class "NorthPolarStereo"
ax = plt.axes(projection=ccrs.NorthPolarStereo(true_scale_latitude=70))

# here is here you tell Cartopy that the projection
# of your 'x' and 'y' are geographic (lons and lats)
# and that you want to transform those lats and lons
# into 'x' and 'y' in the projection
plt.pcolormesh(lon_corners, lat_corners, tile_to_plot[:-1,:-1],
               transform=ccrs.PlateCarree(),shading='flat');

# plot land
ax.add_feature(cfeature.LAND)
ax.gridlines()
ax.coastlines()
plt.colorbar()

# Limit the map to 60 degrees latitude and above.
ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())

# Compute a circle in axes coordinates, which we can use as a boundary
# for the map. We can pan/zoom as much as we like - the boundary will be
# permanently circular.
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)

#ax.set_boundary(circle, transform=ax.transAxes)

Plotting all 13 tiles simultaneously: No Projection¶

Plotting all 13 tiles with `plot_tiles`¶

The plot_tiles routine in the ecco_v4_py package makes plots of all 13 tiles of a field. By default the routine will plot all of the tiles in the lat-lon-cap layout shown earlier.

This routine will accept numpy arrays of dimension 13x90x90 or 2D slices of DataArrays with the same 13x90x90 dimension.

There are several additional arguments which we can access using the help command. Take a second to familiarize yourself with some of them.

[16]:

help(ecco.plot_tiles)

Help on function plot_tiles in module ecco_v4_py.tile_plot:

plot_tiles(tiles, cmap=None, layout='llc', rotate_to_latlon=False, Arctic_cap_tile_location=2, show_colorbar=False, show_cbar_label=False, show_tile_labels=True, cbar_label='', fig_size=9, less_output=True, **kwargs)
    Plots the 13 tiles of the lat-lon-cap (LLC) grid

    Parameters
    ----------
    tiles : numpy.ndarray or dask.array.core.Array or xarray.core.dataarray.DataArray
        an array of n=1..13 tiles of dimension n x llc x llc

            - If *xarray DataArray* or *dask Array* tiles are accessed via *tiles.sel(tile=n)*
            - If *numpy ndarray* tiles are acceed via [tile,:,:] and thus n must be 13.

    cmap : matplotlib.colors.Colormap, optional
        see plot_utils.assign_colormap for default
        a colormap for the figure

    layout : string, optional, default 'llc'
        a code indicating the layout of the tiles

        :llc:    situates tiles in a fan-like manner which conveys how the tiles
                 are oriented in the model in terms of x an y

        :latlon: situates tiles in a more geographically recognizable manner.
                 Note, this does not rotate tiles 7..12, it just places tiles
                 7..12 adjacent to tiles 0..5.  To rotate tiles 7..12
                 specifiy *rotate_to_latlon* as True

    rotate_to_latlon : boolean, default False
        rotate tiles 7..12 so that columns correspond with
        longitude and rows correspond to latitude.  Note, this rotates
        vector fields (vectors positive in x in tiles 7..12 will be -y
        after rotation).

    Arctic_cap_tile_location : int, default 2
        integer, which lat-lon tile to place the Arctic tile over. can be
        2, 5, 7 or 10.

    show_colorbar : boolean, optional, default False
        add a colorbar

    show_cbar_label : boolean, optional, default False
        add a label on the colorbar

    show_tile_labels : boolean, optional, default True
        show tiles numbers in subplot titles

    cbar_label : str, optional, default '' (empty string)
        the label to use for the colorbar

    less_output : boolean, optional, default True
        A debugging flag.  True = less debugging output

    cmin/cmax : floats, optional, default calculate using the min/max of the data
        the minimum and maximum values to use for the colormap

    fig_size : float, optional, default 9 inches
        size of the figure in inches

    fig_num : int, optional, default none
        the figure number to make the plot in.  By default make a new figure.

    Returns
    -------
    f : matplotlib figure object

    cur_arr : numpy ndarray
        numpy array of size:
            (llc*nrows, llc*ncols)
        where llc is the size of tile for llc geometry (e.g. 90)
        nrows, ncols refers to the subplot size
        For now, only implemented for llc90, otherwise None is returned

We’ve seen this routine used in a few earlier tutorials. We’ll provide some additional examples below:

Default ‘native grid’ layout¶

[17]:

# optional arguments:
#   cbar       - show the colorbar
#   cmin, cmax - color range min and max
#   fsize      - figure size in inches

# pull out surface salinity
tmp_plt = ecco_ds.SSS.isel(time=3)

# mask to NaN where hFacC is == 0
# syntax is actually "keep where hFacC is not equal to zero"
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) != 0, np.nan)

[18]:

ecco.plot_tiles(tmp_plt, \
                cmin=32, \
                cmax=35.5, \
                cmap='RdYlBu_r', \
                show_colorbar=True);

# use `suptitle` (super title) to make a title over subplots.
plt.suptitle('SSS ' + str(ecco_ds.time[3].values)[0:7]);

lat-lon layout¶

Another option of plot_tiles is to show tiles 7-12 rotated and lined up tiles 0-5

Note: Rotation of tiles 7-13 is only forplotting. These arrays are not rotated using this routine. We’ll show to how actually rotate these tiles in a later tutorial.

[19]:

# optional arguments:
#   cbar       - show the colorbar
#   cmin, cmax - color range min and max
#   fsize      - figure size in inches

tmp_plt = ecco_ds.SSS.isel(time=3)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) != 0, np.nan)
ecco.plot_tiles(tmp_plt, \
                cmin=32, cmax=35.5, cmap='RdYlBu_r', \
                show_colorbar=True, fig_size=8,\
                layout='latlon', \
                rotate_to_latlon=True);

# use `suptitle` (super title) to make a title over subplots.
plt.suptitle('SSS ' + str(ecco_ds.time[3].values)[0:7]);

The version of plot_tiles is to remove the tile labels and put the titles together in a tight formation and sticks the Arctic tile over tile 10

[20]:

# optional arguments:
#   cbar       - show the colorbar
#   cmin, cmax - color range min and max
#   fsize      - figure size in inches

tmp_plt = ecco_ds.SSS.isel(time=3)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) != 0, np.nan)
ecco.plot_tiles(tmp_plt, cmin=32, cmax=35.5, cmap='RdYlBu_r', \
                show_colorbar=True, fig_size=8,\
               layout='latlon',rotate_to_latlon=True,\
               show_tile_labels=False, \
               Arctic_cap_tile_location=10)

# use `suptitle` (super title) to make a title over subplots.
plt.suptitle('SSS ' + str(ecco_ds.time[3].values)[0:7]);

Almost ready for the hyperwall!

You can plot a subset of tiles using plot_tiles, but you need to pass it a full xarray DataArray or Numpy array with 13 tiles, and the undesired tiles masked out.

[21]:

# optional arguments:
#   cbar       - show the colorbar
#   cmin, cmax - color range min and max
#   fsize      - figure size in inches

tmp_plt = ecco_ds.SSS.isel(time=3)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) != 0, np.nan)
tiles_to_subset = [1,2,4,5]
# add dimensions to vector of tile numbers, so it will readily broadcast across grid
tile_to_broadcast = np.expand_dims(ecco_ds.tile.values,axis=(-2,-1))
# mask tiles not in tiles_to_subset with NaNs
tmp_plt_subset = tmp_plt.where(np.isin(tile_to_broadcast,tiles_to_subset),np.nan)

# select a subset of tiles
ecco.plot_tiles(tmp_plt_subset, cmin=32, cmax=35.5, \
                cmap='RdYlBu_r', show_colorbar=True, fig_size=8,\
                layout='latlon',rotate_to_latlon=True,\
                show_tile_labels=False)

# use `suptitle` (super title) to make a title over subplots.
plt.suptitle('SSS ' + str(ecco_ds.time[3].values)[0:7]);

Plotting all 13 tiles with `plot_proj_to_latlon_grid`¶

Our routine plot_proj_to_latlon_grid takes numpy arrays or DataArrays with 13 tiles and creates global plots with one of three types of projections (passed as arguments to the function): ~~~ projection_type : string, optional denote the type of projection, options include ‘robin’ - Robinson ‘PlateCaree’ - flat 2D projection ‘Mercator’ ‘cyl’ - Lambert Cylindrical ‘ortho’ - Orthographic ‘stereo’ - polar stereographic projection, see lat_lim for choosing ‘InterruptedGoodeHomolosi ~~~

Before plotting this routine interpolates the the filed onto a lat-lon grid (default resoution 0.25 degree) to conform with Cartopy's requirement that the fields to be transformed be on regular square grid.

There are only three argument required of plot_proj_to_latlon_grid, an array of longitudes, an array of latitudes, and an array of the field you wish to plot. The arrays can be either numpy arrays or DataArrays.

Let’s again spend a second to look at the optional arguments available to us in this routine:

[22]:

help(ecco.plot_proj_to_latlon_grid)

Help on function plot_proj_to_latlon_grid in module ecco_v4_py.tile_plot_proj:

plot_proj_to_latlon_grid(lons, lats, data, projection_type='robin', dx=0.25, dy=0.25, radius_of_influence=112000, plot_type='pcolormesh', cmap=None, cmin=None, cmax=None, user_lon_0=0, user_lat_0=None, lat_lim=50, parallels=None, show_coastline=True, show_colorbar=False, show_land=True, show_grid_lines=True, show_grid_labels=False, show_coastline_over_data=True, show_land_over_data=True, grid_linewidth=1, grid_linestyle='--', colorbar_label=None, subplot_grid=None, less_output=True, **kwargs)
    Plot a field of data from an arbitrary grid with lat/lon coordinates
    on a geographic projection after resampling it to a regular lat/lon grid.


    Parameters
    ----------
    lons, lats : numpy ndarray or xarray DataArrays, required
        the longitudes and latitudes of the data to plot

    data : numpy ndarray or xarray DataArray, required
        the field to be plotted

    dx, dy : float, optional, default 0.25 degrees
        latitude, longitude spacing of the new lat/lon grid onto which the
        field 'data' will be resampled.

    radius_of_influence : float, optional, default 112000 m
        to map values from 'data' to the new lat/lon grid, we use use a
        nearest neighbor approach with the constraint that we only use values
        from 'data' that fall within a circle with radius='radius_of_influence'
        from the center of each new lat/lon grid cell.
        for the llc90, with 1 degree resolution,
            radius_of_influence = 1/2 x sqrt(2) x 112e3 km
        would suffice.

    projection_type : string, optional
        denote the type of projection, see Cartopy docs.
        options include
            'robin' - Robinson
            'PlateCarree' - flat 2D projection
            'LambertConformal'
            'Mercator'
            'EqualEarth'
            'Mollweide'
            'AlbersEqualArea'
            'cyl' - Lambert Cylindrical
            'ortho' - Orthographic
            'stereo' - polar stereographic projection, see lat_lim for choosing
            'InterruptedGoodeHomolosine'
                North or South

    plot_type : string, optional
        denotes type of plot ot make with the data
        options include
            'pcolormesh' - pcolormesh
            'contourf' - filled contour
            'points' - plot points at lat/lon locations

    cmap : matplotlib.colors.Colormap, optional, default None
        a colormap for the figure.

    cmin/cmax : floats, optional, default None
        the minimum and maximum values to use for the colormap
        if not specified, use the full range of the data

    user_lon_0 : float, optional, default 0 degrees
        denote central longitude

    user_lat_0 : float, optional, default None
        denote central latitude (for relevant projections only, see Cartopy)

    lat_lim : int, optional, default 50 degrees
        for stereographic projection, denote the Southern (Northern) bounds for
        North (South) polar projection or cutoff for LambertConformal projection

    parallels : float, optional,
        standard_parallels, one or two latitudes of correct scale
        (for relevant projections only, see Cartopy docs)

    show_coastline : logical, optional, default True
        show coastline or not

    show_colorbar : logical, optional, default False
        show a colorbar or not,

    show_land : logical, optional, default True
        show land or not

    show_grid_lines : logical, optional, default True
        True only possible for some cartopy projections

    show_grid_labels: logical, optional, default False
        True only possible for some cartopy projections

    show_coastline_over_data : logical, optional, default True
        draw coastline over the data or under the data

    show_land_over_data: logical, optional, default True
        draw land over the data or under the data

    grid_linewidth : float, optional, default 1.0
        width of grid lines

    grid_linestyle : string, optional, default = '--'
        pattern of grid lines,

    subplot_grid : dict or list, optional
        specifying placement on subplot as
            dict:
                {'nrows': rows_val, 'ncols': cols_val, 'index': index_val}

            list:
                [nrows_val, ncols_val, index_val]

            equates to

                matplotlib.pyplot.subplot(
                    row=nrows_val, col=ncols_val,index=index_val)

    less_output : string, optional
        debugging flag, don't print if True

plot_proj_to_latlon_grid(lons, lats, data, projection_type=’robin’, plot_type=’pcolormesh’, user_lon_0=-66, lat_lim=50, levels=20, cmap=’jet’, dx=0.25, dy=0.25, show_colorbar=False, show_grid_lines=True, show_grid_labels=True, subplot_grid=None, less_output=True, **kwargs)

Robinson projection¶

First we’ll demonstrate the Robinson projection interpolated to a 2x2 degree grid

[23]:

plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, \
                              ecco_ds.YC, \
                              tmp_plt, \
                              plot_type = 'pcolormesh', \
                              dx=2,\
                              dy=2, \
                              projection_type = 'robin',\
                              less_output = False);

_create_projection_axis: projection_type robin
_create_projection_axis: user_lon_0, user_lat_0 0 None
_create_projection_axis: parallels None
_create_projection_axis: lat_lim 50
Projection type:  robin
-179.0 179.0
-180.0 180.0
-89.0 89.0
-90.0 90.0

Setting lon_0 = 110 or -66 yield a global centering that is is more usesful for plotting ocean basins.

[24]:

plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC,
                 ecco_ds.YC,
                 tmp_plt,user_lon_0=-66,
                 plot_type = 'pcolormesh', dx=2,dy=2);

-179.0 113.0
-180.0 114.0
-89.0 89.0
-90.0 90.0
115.03125984375001 178.96875015625
114.00001 180.0
-89.0 89.0
-90.0 90.0

Cylindrical projection¶

Try the Cylindrical Projection with an interpolated lat-lon resolution of 0.25 degrees and pcolormesh.

[25]:

plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, ecco_ds.YC, \
                              tmp_plt, \
                              user_lon_0=-66,\
                              projection_type='cyl',\
                              plot_type = 'pcolormesh', \
                              dx=.25,dy=.25);

-179.875 113.875
-180.0 114.0
-89.875 89.875
-90.0 90.0
114.1254852661597 179.8745247338403
114.00001 180.0
-89.875 89.875
-90.0 90.0

Polar stereographic projection¶

Another projection built into plot_proj_to_latlon_grid is polar stereographic. The argument lat_lim determines the limit of this type of projection. If lat_lim is postive, the projection is centered around the north pole and vice versa.

Northern Hemisphere¶

[26]:

plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, ecco_ds.YC, \
                              tmp_plt, \
                              projection_type='stereo',\
                              plot_type = 'contourf', \
                              show_colorbar=True,
                              dx=1, dy=1,cmin=-1, cmax=1,\
                              lat_lim=40);

-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0

Southern Hemisphere¶

The final example is a south-pole centered plot. Note that lat_lim is now negative.

[27]:

plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, ecco_ds.YC, \
                              tmp_plt, \
                              projection_type='stereo',\
                              plot_type = 'contourf', \
                              show_colorbar=True,
                              dx=1, dy=1,\
                              lat_lim=-40,cmin=-2,cmax=2);

-179.5 179.5
-180.0 180.0
-89.5 89.5
-90.0 90.0

Conclusion¶

You now know several ways of plotting ECCO state estimate fields. There is a lot more to explore with Cartopy - dive in and start making your own cool plots!

[28]:

plt.figure(figsize=(16,6), dpi=90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, ecco_ds.YC, \
                              tmp_plt, \
                              user_lon_0=-66,\
                              projection_type='InterruptedGoodeHomolosine',\
                              plot_type = 'pcolormesh', \
                              show_colorbar=True,
                              dx=1, dy=1);

plt.title('ECCO SSH [m] -- like never before! :)');

-179.5 113.5
-180.0 114.0
-89.5 89.5
-90.0 90.0
114.50770223076924 179.49230776923076
114.00001 180.0
-89.5 89.5
-90.0 90.0

Plotting Tiles¶

Objectives¶

Introduction¶

matplotlib¶

Cartopy¶

The default orientation of the lat-lon-cap tile fields¶

Plotting single tiles using imshow, pcolormesh, and contourf¶

Plotting a single tile with imshow¶

Tile 2 (Northeast Atlantic)¶

Tile 6 (the Arctic cap)¶

Tile 7 (N. Pacific / Bering Sea / Chukchi Sea)¶

Plotting a single tile with pcolor and contourf¶

Tile 2 (Northeast N. Atlantic)¶

Tile 7 (N. Pacific / Bering Sea / Chukchi Sea)¶

Plotting fields from one tile using Cartopy¶

Geographic Projections (AKA: plate carrée)¶

Tile 2 with plate carree¶

Tile 7 with plate carree¶

Polar Stereographic Projection¶

Plotting all 13 tiles simultaneously: No Projection¶

Plotting all 13 tiles with plot_tiles¶

Default ‘native grid’ layout¶

lat-lon layout¶

Plotting all 13 tiles with plot_proj_to_latlon_grid¶

Robinson projection¶

Cylindrical projection¶

Polar stereographic projection¶

Northern Hemisphere¶

Southern Hemisphere¶

Conclusion¶

Plotting single tiles using `imshow`, `pcolormesh`, and `contourf`¶

Plotting a single tile with `imshow`¶

Plotting a single tile with `pcolor` and `contourf`¶

Plotting fields from one tile using `Cartopy`¶

Plotting all 13 tiles with `plot_tiles`¶

Plotting all 13 tiles with `plot_proj_to_latlon_grid`¶