Spatial selection with tripolar ACCESS-OM2 grid

The ACCESS-OM2 model collection utilises a tripolar grid north of 65N to avoid a point of convergence over the ocean at the north pole.

This means any plotting or analysis in this region needs to use 2D curvilinear latitude and longitude coordinates.

This notebook will cover how to use the curvilinear coordinates to select data within latitude and longitude limits, and how to plot data in the tripole region correctly. There is also a full tutorial on plotting with projections.

Below the tripole area it is sufficient to use 1D latitude and longitude coordinates, which are much more convenient to use as they can use the xarray’s sel method with slice notation.

This notebook will also demonstrate how to do this.

The data output from the CICE model component do not have a convenient to use 1D latitude and longitude coordinates. As the ice and ocean models share a grid the coordinates from each model can be used interchangeably with the other. A method to add 1D latitude and longitude coordinates to an ice `xarray.DataArray <http://xarray.pydata.org/en/stable/generated/xarray.DataArray.html>`__ is also demonstrated.

[1]:

import intake
import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import cartopy.crs as ccrs

from dask.distributed import Client

Start dask client

[2]:

client = Client(threads_per_worker=1)
client

/g/data/hh5/public/apps/miniconda3/envs/analysis3-24.04/lib/python3.10/site-packages/distributed/node.py:182: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 38929 instead
  warnings.warn(

[2]:

Client

Client-f0c02835-4588-11ef-a74a-000007aafe80

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /proxy/38929/status

Cluster Info

Load data

[3]:

catalog = intake.cat.access_nri

The specific experiment does not matter a great deal. This is the main IAF control run with JRA55 v1.4 forcing

[4]:

experiment = '01deg_jra55v140_iaf'

Load the 1D latitude/longitude variables from some ocean data to use them with ice data. The dimensions are renamed to match the required dimension in the ice data

[5]:

var_search = catalog[experiment].search(variable='area_t') #area_t is arbitrary, we are only interested in the lat/lon dimension
var_search = var_search.search(path=var_search.df['path'][0]) #workaround to only load from one file
ds = var_search.to_dask()

[6]:

yt_ocean = ds['yt_ocean'].rename({'yt_ocean': 'latitude'})
xt_ocean = ds['xt_ocean'].rename({'xt_ocean': 'longitude'})

Both the ocean and the ice output variables have masked curvilinear coordinates, which means there are missing values for the coordinates over land. These variables are not suitable for plotting with cartopy. A full curvilinear grid is available from the CICE model input grid:

[7]:

ice_grid = xr.open_dataset('/g/data/ik11/inputs/access-om2/input_eee21b65/cice_01deg/grid.nc')

Extract the full curvilinear coordinates, from the ice grid. Note that these grid coordinates are in radians we need to convert them to degrees. Again, the dimensions are renamed to match the dataset to which they will be added:

[8]:

geolon_t = np.degrees(ice_grid.tlon.rename({'nx': 'longitude', 'ny': 'latitude'}))
geolat_t = np.degrees(ice_grid.tlat.rename({'nx': 'longitude', 'ny': 'latitude'}))

Load an ice variable to use as an example dataset. In this case only the first 12 months for purposes of illustration.

[9]:

var_search = catalog[experiment].search(variable='aice_m')
ds = var_search.to_dask(
    xarray_combine_by_coords_kwargs=dict(
        compat="override",
        data_vars="minimal",
        coords="minimal"
    )
)
aice_m = ds['aice_m'].isel(time=slice(0, 12))

Printing the variable shows it has non-informative dimension names and no CF coordinate variables. As a result when the first time slice is plotted there is no useful units automatically added to the plot.

[10]:

aice_m

[10]:

<xarray.DataArray 'aice_m' (time: 12, nj: 2700, ni: 3600)> Size: 467MB
dask.array<getitem, shape=(12, 2700, 3600), dtype=float32, chunksize=(1, 675, 900), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 96B 1958-02-01 1958-03-01 ... 1959-01-01
    TLON     (nj, ni) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
    TLAT     (nj, ni) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
    ULON     (nj, ni) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
    ULAT     (nj, ni) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
Dimensions without coordinates: nj, ni
Attributes:
    units:          1
    long_name:      ice area  (aggregate)
    cell_measures:  area: tarea
    cell_methods:   time: mean
    time_rep:       averaged

Define colormap for plotting. We define the ice_cmap method which accepts two options:

  1. 'ice': Standard ice colormap from cmocean

  2. 'blues': Colormap from dark blue (open ocean) to white (100% sea ice cover)

[11]:

def ice_cmap(colormap):
    """Returns a colormap appropriate for sea ice. Accepts either 'ice' or 'blues'."""

    from palettable.colorbrewer.sequential import Blues_9_r
    import cmocean

    if colormap == 'ice':
        cmap = cmocean.cm.ice
    elif colormap == 'blues':
        cmap = Blues_9_r.mpl_colormap
    else:
        raise ValueError("colormap must be 'ice' or 'blues'")

    cmap.set_bad('lightgrey') # set color for land

    return cmap

[12]:

cmap = ice_cmap('blues')
cmap

[12]:
Blues_r
Blues_r colormap
under
bad
over 
[13]:

aice_m.isel(time=0).plot(aspect=2, size=10, cmap=cmap);
../_images/Tutorials_Spatial_selection_22_0.png

Add spatial coordinates

To improve usability CF compatible spatial coordinates from the ocean data loaded above can be added to the ice data.

The unintuitive dimensions are renamed, and importantly match the dimension names chosen for the coordinates. The coordinates loaded above are added to the DataArray, and saved in a new variable, aice_m_coords.

[14]:

aice_m_coords = aice_m.rename({'nj': 'latitude',
                               'ni': 'longitude'}).assign_coords({'latitude': yt_ocean,
                                                                  'longitude': xt_ocean,
                                                                  'geolat_t': geolat_t,
                                                                  'geolon_t': geolon_t})
aice_m_coords

[14]:

<xarray.DataArray 'aice_m' (time: 12, latitude: 2700, longitude: 3600)> Size: 467MB
dask.array<getitem, shape=(12, 2700, 3600), dtype=float32, chunksize=(1, 675, 900), chunktype=numpy.ndarray>
Coordinates:
  * time       (time) datetime64[ns] 96B 1958-02-01 1958-03-01 ... 1959-01-01
    TLON       (latitude, longitude) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
    TLAT       (latitude, longitude) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
    ULON       (latitude, longitude) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
    ULAT       (latitude, longitude) float32 39MB dask.array<chunksize=(675, 900), meta=np.ndarray>
  * latitude   (latitude) float64 22kB -81.11 -81.07 -81.02 ... 89.94 89.98
  * longitude  (longitude) float64 29kB -279.9 -279.8 -279.7 ... 79.85 79.95
    geolat_t   (latitude, longitude) float64 78MB -81.11 -81.11 ... 65.06 65.02
    geolon_t   (latitude, longitude) float64 78MB -279.9 -279.8 ... 80.0 80.0
Attributes:
    units:          1
    long_name:      ice area  (aggregate)
    cell_measures:  area: tarea
    cell_methods:   time: mean
    time_rep:       averaged

Now when the first time slice is plotted there are informative axes as there are now CF compliant dimensions latitude and longitude.

[15]:

aice_m_coords.isel(time=0).plot(aspect=2, size=10, cmap=cmap);
../_images/Tutorials_Spatial_selection_27_0.png

Selection using sel and slice

Also selection using slices of latitude and longitude are easy using these new dimensions. e.g. selecting the Weddell Sea in Antarctica

[16]:

aice_m_coords.isel(time=0).sel(longitude=slice(-80, 0), latitude=slice(None, -50)).plot(size=10, cmap=cmap);
../_images/Tutorials_Spatial_selection_30_0.png

The same data can also be plotted and projected using cartopy which works fine using the 1D latitude and longitude coordinates as these are correct for all values south of 65N, so all of the Southern Hemisphere.

[17]:

fig = plt.figure(figsize=(14, 8))
ax = plt.axes(projection=ccrs.Orthographic(-40, -90))
ax.coastlines()
ax.gridlines()
aice_m_coords.isel(time=0).sel(longitude=slice(-80, 0),
                               latitude=slice(None, -50)).plot(ax=ax,
                                                               transform=ccrs.PlateCarree(),
                                                               cmap=cmap);
../_images/Tutorials_Spatial_selection_32_0.png

Simple sel with slice can also be used to select areas in the northern hemisphere, but they appear distorted as this area is in the tripole, so the 1D latitude and longitude variables are not correct

[18]:

aice_m_coords.isel(time=0).sel(longitude=slice(-80, 80),
                               latitude=slice(50, 85)).plot(size=8, aspect=2, cmap=cmap);
../_images/Tutorials_Spatial_selection_34_0.png

If the same selection is plotted in an Orthographic projection with cartopy it is obvious that land/sea boundaries are not correct, and artefacts from the tripole are visible

[19]:

fig = plt.figure(figsize=(18, 8))
ax = plt.axes(projection=ccrs.Orthographic(0, 90))
ax.coastlines()
ax.gridlines()
aice_m_coords.isel(time=0).sel(longitude=slice(-80, 80),
                               latitude=slice(50, 85)).plot(ax=ax,
                                                            cmap=cmap,
                                                            transform=ccrs.PlateCarree());
../_images/Tutorials_Spatial_selection_36_0.png

Using the true curvilinear coordinates to plot, the artefacts are removed, and the land/sea registration is correct, but the simple sel using the 1D coordinates means there are not constant lines of latitude or longitude in the selected data.

[20]:

fig = plt.figure(figsize=(18, 8))
ax = plt.axes(projection=ccrs.Orthographic(0,90))
ax.coastlines()
ax.gridlines()
aice_m_coords.isel(time=0).sel(longitude=slice(-80, 80),
                               latitude=slice(50, 85)).plot(ax=ax, cmap=cmap,
                                                            transform=ccrs.PlateCarree(),
                                                            x='geolon_t', y='geolat_t');
../_images/Tutorials_Spatial_selection_38_0.png

Selection using where

To do proper spatial selection with 2D coordinates in the tripole area it is necessary to use the xarray `where <http://xarray.pydata.org/en/stable/generated/xarray.DataArray.where.html>`__ method, which will mask out the data where the stipulated condition is not met. With drop=True coordinate labels with only False values are removed.

[21]:

aice_m_coords.isel(time=0).where((geolon_t >= -80) &
                                 (geolon_t <= 80) &
                                 (geolat_t >= 50) &
                                 (geolat_t <= 85), drop=True).plot(size=8, aspect=2,
                                                                   x='geolon_t',
                                                                   y='geolat_t',
                                                                   cmap=cmap)
plt.xlim([-80, 80])
plt.ylim([50, 85]);
../_images/Tutorials_Spatial_selection_41_0.png

Now when plotted in Orthographic projection with cartopy the selection is clearly along lines of constant latitude, and the left hand side selection is a constant longitude.

[22]:

fig = plt.figure(figsize=(18, 8))
ax = plt.axes(projection=ccrs.Orthographic(0, 90))
ax.coastlines()
ax.gridlines()
aice_m_coords.isel(time=0).where((aice_m_coords.geolon_t >= -80) &
                                 (aice_m_coords.geolon_t <= 80) &
                                 (aice_m_coords.geolat_t >= 50) &
                                 (aice_m_coords.geolat_t <= 85), drop=True).plot(ax=ax,
                                                                                 transform=ccrs.PlateCarree(),
                                                                                 x='geolon_t',
                                                                                 y='geolat_t',
                                                                                 cmap = cmap);
../_images/Tutorials_Spatial_selection_43_0.png 
[23]:

client.close()