A Python tutorial

Mapping Sea Surface Chlorophyll in Python

Using MODIS-Aqua satellite data and xarray

An image of various phytoplankton species from Puget Sound, 2017. Photo taken through the viewer of a light microscope at 400x zoom. Image source: me!
The seasonal cycle of phytoplankton relative to variations in sunlight, nutrients, and zooplankton (copyright 2004 Pearson Prentice Hall, Inc.)
# Imports
import netCDF4 # pip install netCDF4
import xarray as xr # pip install xarray
import cmocean # pip install cmocean
import os
import numpy as np
import matplotlib.pyplot as plt
# Get file path and create list of data files
parent_dir = os.getcwd()
file_path = os.path.join(parent_dir, 'data')
files = [item for item in os.listdir(file_path) if not item.startswith('.')]
# Open the datasets with xarray
datasets = [xr.open_dataset('./data/' + file) for file in files]
# Use .data_vars to find variable names
datasets[0].data_vars
# Generate global snapshot 
datasets[0].chlor_a.plot(x='lon', y='lat', figsize=(26,12), vmin=0, vmax=5, cmap=cmocean.cm.algae);
# Add a title showing the year and month of data collection
plt.title(datasets[0].attrs['time_coverage_start'][:7])
# Enter coordinates of the Gulf of California
site_lat = 26.7
site_lon = -110.7
# Slice the data using the coordinates so that our computer doesn't have to process so much information
ds_slice = datasets[0].sel(lat=slice(site_lat+10, site_lat-10), lon=slice(site_lon-10, site_lon+10))
# Create a plot
ds_slice.chlor_a.plot(x='lon', y='lat', figsize=(12,12), vmin=0, vmax=3, cmap=cmocean.cm.algae);
# Add a title showing the year and month of data collection
plt.title('Gulf of California, ' + datasets[0].attrs['time_coverage_start'][:7])
# Set colorbar values. May take trial and error to get the level of detail you are aiming for.
vmin = 0.0
vmax = 3.0
# Set the lat/lon distance from site location to plot
box_lim = 7
# Set levels for resolution of colorbar. Change the 0.01 value for higher or lower resolution.
lvl = np.arange(vmin, vmax, 0.01).tolist()
# Create a grid of subplots, 3 rows x 4 columns
f, ((ax1, ax2, ax3, ax4), (ax5, ax6, ax7, ax8), (ax9, ax10, ax11, ax12)) = plt.subplots(3, 4, figsize=(24,16))
ax_list = [ax1, ax2, ax3, ax4, ax5, ax6, ax7, ax8, ax9, ax10, ax11, ax12]
# Loop through the subplot axes
for i in range(len(ax_list)):
ds = datasets[i]

# Slice the data using previous coordinates
ds_slice = ds.sel(lat=slice(site_lat+box_lim, site_lat-box_lim), lon=slice(site_lon-box_lim, site_lon+box_lim))

# Generate plot and title
ds_slice.chlor_a.plot.contourf(x='lon', y='lat', ax=ax_list[i], vmin=vmin, vmax=vmax, levels=lvl, cmap=cmocean.cm.algae)
ax_list[i].set_title('Gulf of California, ' + ds.attrs['time_coverage_start'][:7])
Chlorophyll concentrations of the Earth’s waters, collected by the NASA MODIS-Aqua satellite. This image was generated by the NEO website using measurements from February 2021, downloaded here. Lighter colors indicate higher chlorophyll concentrations (up to 60 mg/m³), dark blue indicates low concentration (as low as 0.01 mg/m³), and black represents land, ice, or cloud cover.

Oceanographer turned data scientist, doing my part in science communication! Bringing you my work, insights, and Python tutorials from Seattle, WA.

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store