If you want to know how to plot ocean data in Python, the short answer is this: load it with xarray or pandas, draw it with matplotlib, add geography with cartopy, and colour it with a cmocean colormap. The same four libraries handle a mooring CSV, a CTD cast, an Argo profile and a gridded sea surface temperature field.
That matters because a scripted plot can be regenerated. Change the year, the region or the colour scale, re-run one line, and the figure updates. A manual export from a mapping tool gives you a picture; a script gives you a result you can put in a report and reproduce six months later.
What follows is the order I set up every time: check the data, prepare the columns, draw the simplest figure that works, then add complexity only when the plot asks for it. Allow an hour for your first end-to-end run, mostly spent on setup rather than on plotting itself.
A few quick definitions before the code. A DataArray is a labelled array with named dimensions (time, depth, latitude, longitude). A Dataset is a bundle of those arrays sharing the same coordinates. NetCDF is the file format most satellite, model and mooring products ship in, and the CF conventions describe how those dimensions and units are labelled.
Table of Contents
- What You Need
- Step-by-Step: How to Plot Ocean Data in Python
- Common Mistakes
- Timestamps read as text sort alphabetically
- Celsius and Fahrenheit end up on the same figure
- Depth increases upward on a profile plot
- Straight vertical jumps in a time series
- pcolormesh raises a shape or dimension-order error
- The map is blank or the data never appears
- Data vanishes across the dateline
- The colormap hides the feature you are showing
- A lazy dask array plots as empty
- The chart is unreadable because everything is labelled
- Frequently Asked Questions
- Conclusion
What You Need
You need Python 3.10 or newer, a handful of open-source libraries, and one data file to start with. Nothing here costs money, and nothing needs a licence key.
Two data shapes cover almost everything oceanographic. The first is a row-per-reading file: a CSV log from a mooring, an ADCP, a temperature string or a field logger, with columns like timestamp, depth_m, temperature_c, salinity_psu and sometimes latitude and longitude. The second is a gridded file: NetCDF holding a four-dimensional field of a variable over time, depth, latitude and longitude.
Before plotting anything, print the column names and the first few rows and confirm which fields exist:
import pandas as pd
df = pd.read_csv("mooring_log.csv")
print(df.columns.tolist())
print(df.head())
print(df.isna().sum())
For a gridded file the equivalent check takes two lines:
import xarray as xr
ds = xr.open_dataset("sst_monthly.nc")
print(ds)
That printout lists every variable with its dimensions and units. If the depth axis is labelled in decibars rather than metres, or temperature is in Kelvin, note it now. Unit confusion is the single most expensive mistake in this workflow, and it is free to catch before the first figure.
A notebook works well for exploration because you keep state between cells. When the script runs clean from the top, save it as a .py file next to the data. That script is the reproducible version, and you want it before you write the report.
Step-by-Step: How to Plot Ocean Data in Python
The workflow has six stages: install the stack, load and inspect, prepare the columns, draw a time series, draw depth profiles, then add a map and export. Each stage has a visible output, so you know it worked before moving on. That check is what separates a working script from one that fails three steps later.
1. Install Python and Set Up the Plotting Libraries
Install the libraries into a dedicated environment so a future project cannot break them.
conda create -n ocean python=3.11
conda activate ocean
conda install -c conda-forge xarray numpy matplotlib cartopy cmocean netcdf4 pandas
On Windows and macOS the pip route works too:
python -m venv ocean-env
source ocean-env/bin/activate
pip install xarray numpy matplotlib cartopy cmocean netcdf4 pandas dask
Confirm the imports succeed before touching any data. On Windows, activate with ocean-envScriptsactivate instead.
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cmocean
print(xr.__version__, plt.matplotlib.__version__, ccrs.__name__)
If that line prints versions rather than a traceback, you are ready. On Colab, prefix the install with !pip install and run it in its own cell before any imports.
2. Load and Inspect Ocean Data
Reading a file is one line, but reading it well is where the quality control happens.
csv = pd.read_csv("mooring_log.csv", parse_dates=["timestamp"])
csv = csv.sort_values("timestamp")
print(csv.describe())
print(csv.isna().sum().sort_values(ascending=False).head())
For gridded data, open lazily and look at the coordinate ranges:
ds = xr.open_dataset("sst_monthly.nc")
print(ds.dims)
print(ds.time.values[[0, -1]])
print(float(ds.latitude.min()), float(ds.latitude.max()))
print(float(ds.longitude.min()), float(ds.longitude.max()))
Three things to check before plotting. Are there gaps in the time axis, or a regular interval? Do longitudes run from 0 to 360 or from -180 to 180? Is there a _FillValue or NaN over land? Knowing those three answers explains most of the odd-looking figures you will produce later.
Open datasets stay lazy, so nothing is read from disk until you ask for values. That is why ds.sst.isel(time=0).plot() works on a large file, and why closing the file when you are done matters.
3. Prepare Time, Depth, and Measurement Columns
Preparation means four things: parse the timestamps, sort by time, settle the units, and write down what you dropped.
csv["timestamp"] = pd.to_datetime(csv["timestamp"], utc=True)
csv = csv.sort_values("timestamp")
# Kelvin to Celsius, if the logger recorded Kelvin
if csv["temperature_c"].max() > 100:
csv["temperature_c"] = csv["temperature_c"] - 273.15
# Pressure in decibars to depth in metres, roughly 1 dbar per metre
csv["depth_m"] = csv["depth_dbar"] / 1.0
# Flagged readings come out before plotting, and the count goes in your notes
dropped = csv[csv["flag"] != 1]
csv = csv[csv["flag"] == 1]
print(f"dropped {len(dropped)} flagged rows")
Keep that exclusion count. When a reviewer asks why the figure ends in November, the number in your notes answers it.
For gridded data, sel picks by label and isel picks by position:
slice1 = ds.sst.sel(time="2024-07-01", method="nearest")
slice2 = ds.sst.isel(time=6)
region = ds.sst.sel(latitude=slice(30, 60), longitude=slice(-20, 10))
Use sel when you know the date or the coordinates and want the figure to keep working when the file updates. Use isel when you mean “the seventh frame”, such as the last month in the record. method="nearest" matters because satellite passes are not always on the first of the month.
4. Plot a Time Series of Ocean Measurements
A time series is the fastest way to see whether your data is worth plotting at all.
fig, ax = plt.subplots(figsize=(10, 4.5))
ax.plot(csv["timestamp"], csv["temperature_c"], linewidth=1.0, label="Temperature")
ax.set_xlabel("Date")
ax.set_ylabel("Temperature (degC)")
ax.set_title("Mooring surface temperature, 2024")
ax.grid(True, alpha=0.3)
ax.legend()
fig.autofmt_xdate()
plt.tight_layout()
plt.show()
Read the figure rather than admire it. A straight vertical jump between two points means a gap in the record, not a real warming event. A flat line over several hours often means the sensor stopped reporting. A daily sinusoid that matches the tide is a working instrument, not a broken one.
Two measurements on one axis need a second y axis, and you should only use it when both share units:
fig, ax1 = plt.subplots(figsize=(10, 4.5))
ax1.plot(csv["timestamp"], csv["temperature_c"], color="crimson", label="Temperature")
ax1.set_ylabel("Temperature (degC)", color="crimson")
ax2 = ax1.twinx()
ax2.plot(csv["timestamp"], csv["salinity_psu"], color="teal", label="Salinity")
ax2.set_ylabel("Salinity (PSU)", color="teal")
ax1.grid(True, alpha=0.3)
fig.autofmt_xdate()
plt.tight_layout()
Temperature in degC and salinity in practical salinity units do not belong on one axis. Two separate stacked panels are clearer.
5. Plot Temperature and Salinity by Depth
Depth profiles show the vertical structure: the mixed layer on top, the thermocline beneath it, and near-constant deep water below.
cast = csv[csv["cast_id"] == cast_number]
fig, ax = plt.subplots(figsize=(5.5, 7))
ax.plot(cast["temperature_c"], cast["depth_m"], marker="o", markersize=3)
ax.invert_yaxis()
ax.set_xlabel("Temperature (degC)")
ax.set_ylabel("Depth (m)")
ax.set_title(f"CTD cast {cast_number}")
ax.grid(True, alpha=0.3)
plt.tight_layout()
Call invert_yaxis() and depth now increases downward, matching every CTD plot you have seen in a paper. Skip it and the thermocline appears to sit above the surface.
Overlaying several casts in separate colours is how you show an upwelling event:
fig, ax = plt.subplots(figsize=(6, 7))
colours = plt.cm.viridis(np.linspace(0, 1, n_casts))
for colour, (cast_id, group) in zip(colours, csv.groupby("cast_id")):
ax.plot(group["temperature_c"], group["depth_m"], color=colour, linewidth=1.2, label=str(cast_id))
ax.invert_yaxis()
ax.set_xlabel("Temperature (degC)")
ax.set_ylabel("Depth (m)")
ax.legend(fontsize="small", title="Cast")
ax.grid(True, alpha=0.3)
plt.tight_layout()
A temperature-salinity diagram does the same job in one view and separates water masses by their shape:
fig, ax = plt.subplots(figsize=(6, 6))
sc = ax.scatter(csv["salinity_psu"], csv["temperature_c"], c=csv["depth_m"], s=12, cmap=cmocean.cm.deep)
fig.colorbar(sc, ax=ax, label="Depth (m)")
ax.set_xlabel("Salinity (PSU)")
ax.set_ylabel("Temperature (degC)")
ax.set_title("Temperature-salinity diagram")
ax.grid(True, alpha=0.3)
Colour by depth and you can read which water mass sits deepest without a legend of twelve lines.
6. Add a Map and Export the Final Figures
Maps are where most ocean plots end up, and cartopy handles the projection and coastlines.
sst = ds.sst.sel(time="2024-07-01", method="nearest")
fig, ax = plt.subplots(figsize=(11, 6), subplot_kw={"projection": ccrs.PlateCarree()})
mesh = ax.pcolormesh(
sst.longitude, sst.latitude, sst,
transform=ccrs.PlateCarree(),
cmap=cmocean.cm.thermal,
shading="auto",
)
ax.add_feature(cfeature.COASTLINE.with_scale("50m"), linewidth=0.6)
ax.add_feature(cfeature.BORDERS.with_scale("50m"), linewidth=0.3)
fig.colorbar(mesh, ax=ax, orientation="vertical", shrink=0.75, label="Sea surface temperature (degC)")
ax.set_title("Sea surface temperature, July 2024")
plt.tight_layout()
The transform=ccrs.PlateCarree() argument is not optional. Your data sits in longitude and latitude coordinates; the axes sit in whatever projection you chose. Without that line cartopy either draws nothing or throws.
Current vectors need speed as well as direction, and np.hypot computes the magnitude cleanly:
u = ds.uo.sel(time="2024-07-01", method="nearest")
v = ds.vo.sel(time="2024-07-01", method="nearest")
speed = np.hypot(u, v)
speed_kmh = speed * 3.6
ax.quiver(
u.longitude, u.latitude, u, v,
transform=ccrs.PlateCarree(),
color="white", scale=25, width=0.003,
)
qk = ax.quiverkey(u, 1.0, 1.02, 0.5, "0.5 m/s", labelpos="E", coordinates="axes")
Thin the vectors by taking every third point, otherwise a global grid turns into solid black. The quiverkey gives the reader a scale instead of leaving them to guess.
Export at full resolution and in a vector format for print:
fig.savefig("sst_july2024.png", dpi=300, bbox_inches="tight")
fig.savefig("sst_july2024.pdf", bbox_inches="tight")
plt.close()
For large files, chunk the dataset with ds.chunk({"time": 1}) so it stays lazy and only the slice you asked for is loaded.
Common Mistakes
Almost every bad ocean plot comes from one of the problems below, and each has a quick fix.
Timestamps read as text sort alphabetically
Pass parse_dates=["timestamp"] or call pd.to_datetime(), then check df.dtypes before plotting. A datetime axis cannot be formatted or resampled until the column is genuinely datetime.
Celsius and Fahrenheit end up on the same figure
Convert everything to degC in the preparation step and label the axis with units. Record the conversion in a comment, because you will forget which file was already in Celsius.
Depth increases upward on a profile plot
Call ax.invert_yaxis() once per axes. If the thermocline looks like it floats above the surface, this is why.
Straight vertical jumps in a time series
Those are gaps, not events. Use ax.plot rather than a bar or step chart, and mark the gap with csv["timestamp"].diff() so you know its length before you interpret anything.
pcolormesh raises a shape or dimension-order error
The array must be the last argument and its dimensions must run longitude by latitude. Fix it with field.transpose("longitude", "latitude") before plotting, and check field.dims first.
The map is blank or the data never appears
Add transform=ccrs.PlateCarree(). Also check that your longitude values fall inside the axes limits, since a 0-to-360 file plotted on axes set to -180 to 180 leaves an empty figure.
Data vanishes across the dateline
Normalise once: ds = ds.assign_coords(longitude=(((ds.longitude + 180) % 360) - 180)) and then sort the coordinate. Mixing the two conventions in one session is the usual cause.
The colormap hides the feature you are showing
Skip jet. It invents edges that are not in the data and most viewers cannot perceive the ordering reliably. Use cmocean.cm.thermal for temperature, haline for salinity, balance for anomalies and topo for bathymetry. Set the range explicitly with vmin and vmax so two figures are comparable.
A lazy dask array plots as empty
Chunked datasets stay lazy, and a slice with dask types sometimes needs computing before matplotlib sees values. Call .load() on the subset, or assign the result back. Only for the slice you are plotting, never for the whole file.
The chart is unreadable because everything is labelled
Tick labels on a dense grid do not scale. Show every nth label, thin the vectors, drop the minor ticks, and keep the units on the colorbar where readers look for them first.
Frequently Asked Questions
What is xarray and why is it used for ocean data?
xarray labels an array with its dimension names and coordinates, so you can pick a time or depth by name instead of counting rows. Ocean products arrive as gridded NetCDF files whose dimensions are time, depth, latitude and longitude, and xarray keeps those labels attached through every calculation. That is what lets you ask for one month of one depth across one bounding box in a single line.
Should I use basemap or cartopy?
Use cartopy. Basemap is deprecated and no longer receives updates, yet several older tutorials still rank well and send beginners into code that no longer installs cleanly. cartopy draws coastlines and projections inside matplotlib, so figure handling, colorbars and saving work the way you already expect. Translating a legacy Basemap script takes minutes: swap the draw calls for cartopy features and add the transform argument.
What is cmocean and why use it instead of jet?
cmocean is a set of colormaps built for oceanographic variables, with perceptually uniform steps so an equal visual difference means an equal data difference. jet creates bright artificial boundaries where no gradient exists, which makes a smooth thermocline look like a step change. Use thermal for temperature, haline for salinity, balance for anomalies and topo for bathymetry.
Is map lazy in Python?
Yes, in the sense that big ocean files stay lazy by default. Opening a NetCDF file with xarray reads only metadata, and a dask-chunked dataset computes values only when you index, plot or save. That keeps large files off your memory, but it means a slice can appear empty until something forces a compute. Call .load() on the small subset you are about to plot.
Where do I get ocean data to plot?
Satellite and reanalysis fields come as NetCDF from services such as Copernicus Marine, NOAA CoastWatch and ERDDAP. The World Ocean Atlas and Argo float profiles cover climatology and the upper 2000 metres, with the argopy package pulling Argo directly from a floating-point profile selector. For your own instruments, a mooring, ADCP or temperature string usually exports a CSV, which pandas reads without conversion.
How do I plot a very large NetCDF file without crashing?
Open it with xarray and chunk it before plotting, for example ds.chunk({“time”: 1}), which keeps the data lazy. Then slice to the region, time window and depth range you actually need before drawing, and pass method=“nearest” when your timestamps fall between coordinate values. Also check that the variable is a float rather than an unexpected integer type, and close the dataset when you finish.
Conclusion
Start with one file. Load it, print the column names and the first rows, confirm your timestamps parsed and your units are what the labels claim, then draw one time series with units on both axes. That single figure tells you whether the record is continuous and whether the sensor is alive.
From there the next steps are mechanical: a depth profile with an inverted axis for casts, a cmocean colormap and a cartopy map with coastlines for anything spatial. Save every figure at 300 dpi as both PNG and PDF, keep the script beside the data, and you have a result you can rerun and defend. That is the whole workflow, and once the first plot runs the rest is repetition.


