Vector plotting methods

To plot vector quantities such as currents we can use a quiver plot in matplotlib. In the GBR4 tutorial dataset the u and v variables contain the x and y components of the vector for each cell in the dataset.

import matplotlib.pyplot as plt
import numpy
import pandas
import xarray
from shapely import box

import emsarray
from emsarray import plot
from emsarray.operations import point_extraction

dataset = emsarray.tutorial.open_dataset('gbr4')

surface_currents = dataset.ems.select_variables(["temp", "u", "v"]).isel(time=0, k=-1)
surface_currents.load()

<xarray.Dataset> Size: 4MB
Dimensions:    (j: 180, i: 600)
Coordinates:
    latitude   (j, i) float64 864kB -8.074 -8.107 -8.139 ... -28.56 -28.59
    longitude  (j, i) float64 864kB 142.2 142.2 142.2 ... 156.9 156.9 156.9
    time       datetime64[ns] 8B 2022-10-15T14:00:00
    zc         float64 8B 1.5
Dimensions without coordinates: j, i
Data variables:
    temp       (j, i) float64 864kB nan nan nan nan ... 20.9 20.92 20.93 20.93
    u          (j, i) float64 864kB nan nan nan ... -0.005159 0.03509 0.05466
    v          (j, i) float64 864kB nan nan nan nan ... 0.08804 0.08797 0.04464
Attributes: (12/16)
    title:                           gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_Dhnd
    description:                     eReefs GBR 4k grid with 64 active river ...
    paramhead:                       eReefs GBR 4k grid, with Atmospheric for...
    paramfile:                       /g/data/et4/projects/eReefs/gbr4_H4p0_AB...
    ems_version:                     v1.4.0 rev(7365)
    date_created:                    Mon Oct 21 15:47:23 2024
    ...                              ...
    bald__isPrefixedBy:              prefix_list
    metadata_link:                   https://marlin.csiro.au/geonetwork/srv/e...
    prefix_list_puv__:               https://w3id.org/env/puv#
    DODS_EXTRA.Unlimited_Dimension:  time
    history:                         Wed Jan 14 14:15:20 2026: ncks -d time,0...
    NCO:                             netCDF Operators version 5.3.6 (Homepage...

xarray.Dataset

• Dimensions:
  • j: 180
  • i: 600
• Coordinates: (4)
  • latitude
    (j, i)
    float64
    -8.074 -8.107 ... -28.56 -28.59

    puv__uom :

    http://vocab.nerc.ac.uk/collection/P06/current/DEGN/

    units :

    degrees_north

    long_name :

    Latitude

    standard_name :

    latitude

    coordinate_type :

    latitude

    projection :

    geographic

    array([[ -8.0741344,  -8.1069128,  -8.1390488, ...,         nan,
                    nan,         nan],
           [ -8.063217 ,  -8.0957244,  -8.1278305, ...,         nan,
                    nan,         nan],
           [ -8.0526946,  -8.0849348,  -8.1169279, ...,         nan,
                    nan,         nan],
           ...,
           [ -7.0133871,  -7.0358495,  -7.0584407, ..., -28.5233981,
            -28.5594013, -28.5953144],
           [ -7.0127399,  -7.0348327,  -7.0571806, ..., -28.5227067,
            -28.5586112, -28.5943624],
           [ -7.0119082,  -7.0335482,  -7.0556973, ..., -28.5221398,
            -28.5579669, -28.5935179]], shape=(180, 600))

  • longitude
    (j, i)
    float64
    142.2 142.2 142.2 ... 156.9 156.9

    puv__uom :

    http://vocab.nerc.ac.uk/collection/P06/current/DEGE/

    units :

    degrees_east

    long_name :

    Longitude

    standard_name :

    longitude

    coordinate_type :

    longitude

    projection :

    geographic

    array([[142.168788 , 142.1804468, 142.19211  , ...,         nan,
                    nan,         nan],
           [142.2013204, 142.2123752, 142.2236158, ...,         nan,
                    nan,         nan],
           [142.2333969, 142.2441535, 142.25509  , ...,         nan,
                    nan,         nan],
           ...,
           [147.8135068, 147.8142988, 147.8153794, ..., 156.8042355,
            156.8051407, 156.8061819],
           [147.8353694, 147.8363158, 147.8375886, ..., 156.84433  ,
            156.8451114, 156.8460513],
           [147.8566896, 147.8580363, 147.8596273, ..., 156.8843511,
            156.8849406, 156.8856378]], shape=(180, 600))

  • time
    ()
    datetime64[ns]
    2022-10-15T14:00:00

    long_name :

    Time

    standard_name :

    time

    coordinate_type :

    time

    puv__uom :

    http://vocab.nerc.ac.uk/collection/P06/current/UTAA/

    _ChunkSizes :

    512

    array('2022-10-15T14:00:00.000000000', dtype='datetime64[ns]')

  • zc
    ()
    float64
    1.5

    puv__uom :

    http://vocab.nerc.ac.uk/collection/P06/current/ULAA/

    units :

    m

    positive :

    up

    long_name :

    Z coordinate

    axis :

    Z

    coordinate_type :

    Z

    array(1.5)

• Data variables: (3)
  • temp
    (j, i)
    float64
    nan nan nan ... 20.92 20.93 20.93

    _ChunkSizes :

    [ 1 24 90 300]

    long_name :

    Temperature

    puv__parameter :

    http://vocab.nerc.ac.uk/collection/P01/current/TEMPMP01/

    units :

    degrees C

    valid_range :

    [ 0. 40.]

    array([[        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           ...,
           [        nan,         nan,         nan, ..., 20.92445135,
            20.93568182, 20.93665838],
           [        nan,         nan,         nan, ..., 20.92445135,
            20.93421698, 20.93714666],
           [        nan,         nan,         nan, ..., 20.91517401,
            20.93177557, 20.93324041]], shape=(180, 600))

  • u
    (j, i)
    float64
    nan nan nan ... 0.03509 0.05466

    _ChunkSizes :

    [ 1 24 90 300]

    long_name :

    Eastward current

    standard_name :

    eastward_sea_water_velocity

    vector_name :

    Currents

    puv__parameter :

    http://vocab.nerc.ac.uk/collection/P01/current/LCEWMP01/

    units :

    ms-1

    vector_components :

    u v

    valid_range :

    [-100. 100.]

    array([[        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           ...,
           [        nan,         nan,         nan, ..., -0.00702071,
             0.06094193,  0.10937333],
           [        nan,         nan,         nan, ..., -0.00833297,
             0.06732011,  0.1100142 ],
           [        nan,         nan,         nan, ..., -0.00515914,
             0.03509355,  0.05465531]], shape=(180, 600))

  • v
    (j, i)
    float64
    nan nan nan ... 0.08797 0.04464

    _ChunkSizes :

    [ 1 24 90 300]

    long_name :

    Northward current

    standard_name :

    northward_sea_water_velocity

    vector_name :

    Currents

    puv__parameter :

    http://vocab.nerc.ac.uk/collection/P01/current/LCNSMP01/

    units :

    ms-1

    vector_components :

    u v

    valid_range :

    [-100. 100.]

    array([[       nan,        nan,        nan, ...,        nan,        nan,
                   nan],
           [       nan,        nan,        nan, ...,        nan,        nan,
                   nan],
           [       nan,        nan,        nan, ...,        nan,        nan,
                   nan],
           ...,
           [       nan,        nan,        nan, ..., 0.03984872, 0.04570809,
            0.0269703 ],
           [       nan,        nan,        nan, ..., 0.04747811, 0.0544056 ,
            0.03197518],
           [       nan,        nan,        nan, ..., 0.08803597, 0.08797494,
            0.04463997]], shape=(180, 600))

• Attributes: (16)

  title :

  gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_Dhnd

  description :

  eReefs GBR 4k grid with 64 active river flows. Hydro 4.0.

  paramhead :

  eReefs GBR 4k grid, with Atmospheric forcing using BARRAr2, open Ocean forcing using BRAN2020, river Flows forcing using G2G and QDES SOURCE; Deployed as a hindcast. Hydro GBR4 parameter file 4.0.

  paramfile :

  /g/data/et4/projects/eReefs/gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_Dhnd/tran/gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3.tran

  ems_version :

  v1.4.0 rev(7365)

  date_created :

  Mon Oct 21 15:47:23 2024

  Conventions :

  CF-1.0

  Run_ID :

  4.0

  Run_code :

  SHOC grid|G3.00|H4.0|S-3.00|B0.00

  Parameter_File_Revision :

  $Revision$

  bald__isPrefixedBy :

  prefix_list

  metadata_link :

  https://marlin.csiro.au/geonetwork/srv/eng/catalog.search#/metadata/dcc8462c-8dc5-4d78-b9d8-6f1c31ae0001

  prefix_list_puv__ :

  https://w3id.org/env/puv#

  DODS_EXTRA.Unlimited_Dimension :

  time

  history :

  Wed Jan 14 14:15:20 2026: ncks -d time,0,0 https://thredds.nci.org.au/thredds/dodsC/fx3/gbr4_H4p0_ABARRAr2_OBRAN2020_FG2Gv3_Dhnd/gbr4_simple_2022-10-16.nc gbr4.temp.nc

  NCO :

  netCDF Operators version 5.3.6 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco, Citation = 10.1016/j.envsoft.2008.03.004)

Plotting every cell is straight forward for small areas:

figure = plt.figure(figsize=(7, 8), layout="constrained")
axes = figure.add_subplot(projection=dataset.ems.data_crs)
figure.suptitle("Surface water temperature and currents near the Whitsundays")

temp = surface_currents["temp"]
temp_artist = dataset.ems.make_artist(
    axes, temp,
    cmap="coolwarm", clim=(24, 27))

current_artist = dataset.ems.make_artist(
    axes, (surface_currents["u"], surface_currents["v"]),
    scale=20, scale_units="width")

plot.add_coast(axes)
plot.add_gridlines(axes)

# Just show a small area over the Whitsundays
view_box = box(148.7, -20.4, 149.6, -19.8)
axes.set_aspect("equal", adjustable="datalim")
axes.set_extent(plot.bounds_to_extent(view_box.bounds))

Plotting the entire dataset this way leads to the current vectors becoming a confusing mess however:

figure = plt.figure(figsize=(7, 10), layout="constrained")
axes = figure.add_subplot(projection=dataset.ems.data_crs)
figure.suptitle("A bad plot of surface water temperature and currents over the entire reef")

temp = surface_currents["temp"]
temp_artist = dataset.ems.make_artist(axes, temp, cmap="coolwarm")

current_artist = dataset.ems.make_artist(
    axes, (surface_currents["u"], surface_currents["v"]),
    scale=20, scale_units="width")

plot.add_coast(axes)
plot.add_gridlines(axes)

# Show the entire model domain
axes.set_aspect("equal", adjustable="datalim")
axes.autoscale()

Sampling gridded datasets

For datasets on a two dimensional grid like this we can sample the current data at regular intervals to display only a subset of the vectors. The sampled points will follow the geometry of the dataset. For curvilinear datasets such as GBR4 the vectors will follow the curves of the dataset shape:

# Make an empty array of the same shape as the data, then select every 10th cell in there.
face_grid = dataset.ems.get_grid(surface_currents["u"])
samples = xarray.DataArray(numpy.full(face_grid.size, False))
samples = face_grid.wind(samples)
samples[::10, ::10] = True
samples = dataset.ems.ravel(samples)

# Select the (x, y) coordinates and the (u, v) components of the sampled cells
x, y = face_grid.centroid_coordinates[samples].T
u = dataset.ems.ravel(surface_currents["u"]).values[samples]
v = dataset.ems.ravel(surface_currents["v"]).values[samples]

# Plot the points
figure = plt.figure(figsize=(7, 10), layout="constrained")
axes = figure.add_subplot(projection=dataset.ems.data_crs)
figure.suptitle("Surface water temperature and currents across the entire dataset domain")

temp = surface_currents["temp"]
temp_artist = dataset.ems.make_artist(axes, temp, cmap="coolwarm")

quiver = plt.quiver(
    x, y, u, v,
    scale=40, scale_units="width")
axes.add_collection(quiver)

plot.add_coast(axes)
plot.add_gridlines(axes)

# Show the entire model domain
axes.set_aspect("equal", adjustable="datalim")
axes.autoscale()

Sampling the dataset domain

Another approach is to plot vectors at regular points across the dataset domain by sampling at regular intervals. The vector locations are not tied to the grid geometry. This approach will work with unstructured grids unlike the previous method. point_extraction.extract_dataframe() can be used for this:

# Generate a mesh of points across the dataset domain
domain = box(*dataset.ems.bounds)
x = numpy.arange(domain.bounds[0], domain.bounds[2], 0.4)
y = numpy.arange(domain.bounds[1], domain.bounds[3], 0.4)
xx, yy = numpy.meshgrid(x, y)
points = pandas.DataFrame({
    'x': xx.flatten(),
    'y': yy.flatten(),
})

# Extract the surface current components at these locations
components = point_extraction.extract_dataframe(
    surface_currents, points, ('x', 'y'), missing_points='drop')

# Plot the points
figure = plt.figure(figsize=(7, 10), layout="constrained")
axes = figure.add_subplot(projection=dataset.ems.data_crs)
figure.suptitle("Surface water temperature and currents across the entire dataset domain")

temp = surface_currents["temp"]
temp_artist = dataset.ems.make_artist(axes, temp, cmap="coolwarm")

quiver = plt.quiver(
    components['x'], components['y'], components['u'], components['v'],
    scale=30, scale_units="width")
axes.add_collection(quiver)

plot.add_coast(axes)
plot.add_gridlines(axes)

axes.set_aspect("equal", adjustable="datalim")
axes.autoscale()