Zooplankton biomass

A copepod, a type of zooplankton. Art credit: Kristen Krumhardt

Overview

Zooplankton are tiny oceanic animals, making up the next step up after phytoplankton in the food web. Here we evaluate modeled zooplankton biomass and compare it to observational data.

  1. General setup

  2. Subsetting

  3. Processing - long-term mean

  4. Mapping zooplankton biomass at the surface

  5. Comparing mesozooplankton biomass to observations

  6. Making monthly climatology maps to compare to observations

Prerequisites

Concepts

Importance

Notes

Matplotlib

Necessary

Intro to Cartopy

Necessary

Dask Cookbook

Helpful

Intro to Xarray

Helpful

  • Time to learn: 30 min

Imports

import xarray as xr
import glob
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import pop_tools
from dask.distributed import LocalCluster
import s3fs
import netCDF4

from module import adjust_pop_grid

General setup (see intro notebooks for explanations)

Connect to cluster

cluster = LocalCluster()
client = cluster.get_client()

Bring in POP grid utilities

ds_grid = pop_tools.get_grid('POP_gx1v7')
lons = ds_grid.TLONG
lats = ds_grid.TLAT
depths = ds_grid.z_t * 0.01

Load the data

jetstream_url = 'https://js2.jetstream-cloud.org:8001/'

s3 = s3fs.S3FileSystem(anon=True, client_kwargs=dict(endpoint_url=jetstream_url))

# Generate a list of all files in CESM folder
s3path = 's3://pythia/ocean-bgc/cesm/g.e22.GOMIPECOIAF_JRA-1p4-2018.TL319_g17.4p2z.002branch/ocn/proc/tseries/month_1/*'
remote_files = s3.glob(s3path)

# Open all files from folder
fileset = [s3.open(file) for file in remote_files]

# Open with xarray
ds = xr.open_mfdataset(fileset, data_vars="minimal", coords='minimal', compat="override", parallel=True,
                       drop_variables=["transport_components", "transport_regions", 'moc_components'], decode_times=True)

ds

<xarray.Dataset> Size: 28GB
Dimensions:                         (nlat: 384, nlon: 320, time: 120, z_t: 60,
                                     z_t_150m: 15)
Coordinates:
    TLAT                            (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
    TLONG                           (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
  * time                            (time) object 960B 2010-01-16 12:00:00 .....
  * z_t                             (z_t) float32 240B 500.0 ... 5.375e+05
  * z_t_150m                        (z_t_150m) float32 60B 500.0 ... 1.45e+04
Dimensions without coordinates: nlat, nlon
Data variables: (12/45)
    FG_CO2                          (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    Fe                              (time, z_t, nlat, nlon) float32 4GB dask.array<chunksize=(30, 15, 96, 80), meta=np.ndarray>
    NO3                             (time, z_t, nlat, nlon) float32 4GB dask.array<chunksize=(30, 15, 96, 80), meta=np.ndarray>
    PO4                             (time, z_t, nlat, nlon) float32 4GB dask.array<chunksize=(30, 15, 96, 80), meta=np.ndarray>
    POC_FLUX_100m                   (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    SALT                            (time, z_t, nlat, nlon) float32 4GB dask.array<chunksize=(30, 15, 96, 80), meta=np.ndarray>
    ...                              ...
    sp_Fe_lim_Cweight_avg_100m      (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    sp_Fe_lim_surf                  (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    sp_N_lim_Cweight_avg_100m       (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    sp_N_lim_surf                   (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    sp_P_lim_Cweight_avg_100m       (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>
    sp_P_lim_surf                   (time, nlat, nlon) float32 59MB dask.array<chunksize=(60, 192, 160), meta=np.ndarray>

Subsetting

variables =['mesozooC', 'microzooC']

keep_vars=['z_t','z_t_150m','dz','time_bound','time','TAREA','TLAT','TLONG'] + variables
ds = ds.drop_vars([v for v in ds.variables if v not in keep_vars])

Processing - long-term mean

Pull in the function we defined in the nutrients notebook…

def year_mean(ds):
    """
    Source: https://ncar.github.io/esds/posts/2021/yearly-averages-xarray/
    """
    
    # Make a DataArray with the number of days in each month, size = len(time)
    month_length = ds.time.dt.days_in_month

    # Calculate the weights by grouping by 'time.year'
    weights = (
        month_length.groupby("time.year") / month_length.groupby("time.year").sum()
    )

    # Test that the sum of the weights for each year is 1.0
    np.testing.assert_allclose(weights.groupby("time.year").sum().values, np.ones((len(ds.groupby("time.year")), )))

    # Calculate the weighted average
    return (ds * weights).groupby("time.year").sum(dim="time")

# Take the long-term mean of our data set, processing years and months separately

ds_annual = year_mean(ds).mean("year")

Plot mesozooplankton and microzooplankton biomass at the surface

fig = plt.figure(figsize=(8,5))

ax = fig.add_subplot(2,1,1, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('microzooC at surface', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  ds_annual.microzooC.isel(z_t_150m=0))
pc=ax.pcolormesh(lon, lat, field, cmap='Blues',vmin=0,vmax=2,transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(pc, ax=ax,extend='max',label='microzooC (mmol m$^{-3}$)')
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)


ax = fig.add_subplot(2,1,2, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('mesozooC at surface', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  ds_annual.mesozooC.isel(z_t_150m=0))
pc=ax.pcolormesh(lon, lat, field, cmap='Oranges',vmin=0,vmax=4,transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(pc, ax=ax,extend='max',label='mesozooC (mmol m$^{-3}$)')
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

<cartopy.mpl.feature_artist.FeatureArtist at 0x7ff6acf29c40>
../_images/1a2d76bb70818d056a87caa8e4bf5f181d15ab747fc5f03b5b647ccbbe6c63ec.png

Compare mesozooplankton biomass to COPEPOD database

We use data compiled through the COPEPOD project (Moriarty & O’Brien, 2013). This data has been pre-processed, but the raw data is available on the COPEPOD website.

Read in COPEPOD data

copepod_obs_path = 's3://pythia/ocean-bgc/obs/copepod-2012__cmass-m00-qtr.zarr'

copepod_obs = s3fs.S3Map(root=copepod_obs_path, s3=s3)

ds_copepod = xr.open_dataset(copepod_obs, engine="zarr")

### converting grams to moles
ds_copepod['copepod_C']=ds_copepod.copepod_C/12.011

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[10], line 5
      1 copepod_obs_path = 's3://pythia/ocean-bgc/obs/copepod-2012__cmass-m00-qtr.zarr'
      3 copepod_obs = s3fs.S3Map(root=copepod_obs_path, s3=s3)
----> 5 ds_copepod = xr.open_dataset(copepod_obs, engine="zarr")
      7 ### converting grams to moles
      8 ds_copepod['copepod_C']=ds_copepod.copepod_C/12.011

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/xarray/backends/api.py:679, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
    667 decoders = _resolve_decoders_kwargs(
    668     decode_cf,
    669     open_backend_dataset_parameters=backend.open_dataset_parameters,
   (...)
    675     decode_coords=decode_coords,
    676 )
    678 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 679 backend_ds = backend.open_dataset(
    680     filename_or_obj,
    681     drop_variables=drop_variables,
    682     **decoders,
    683     **kwargs,
    684 )
    685 ds = _dataset_from_backend_dataset(
    686     backend_ds,
    687     filename_or_obj,
   (...)
    697     **kwargs,
    698 )
    699 return ds

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/xarray/backends/zarr.py:1564, in ZarrBackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, synchronizer, consolidated, chunk_store, storage_options, zarr_version, zarr_format, store, engine, use_zarr_fill_value_as_mask, cache_members)
   1562 filename_or_obj = _normalize_path(filename_or_obj)
   1563 if not store:
-> 1564     store = ZarrStore.open_group(
   1565         filename_or_obj,
   1566         group=group,
   1567         mode=mode,
   1568         synchronizer=synchronizer,
   1569         consolidated=consolidated,
   1570         consolidate_on_close=False,
   1571         chunk_store=chunk_store,
   1572         storage_options=storage_options,
   1573         zarr_version=zarr_version,
   1574         use_zarr_fill_value_as_mask=None,
   1575         zarr_format=zarr_format,
   1576         cache_members=cache_members,
   1577     )
   1579 store_entrypoint = StoreBackendEntrypoint()
   1580 with close_on_error(store):

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/xarray/backends/zarr.py:703, in ZarrStore.open_group(cls, store, mode, synchronizer, group, consolidated, consolidate_on_close, chunk_store, storage_options, append_dim, write_region, safe_chunks, zarr_version, zarr_format, use_zarr_fill_value_as_mask, write_empty, cache_members)
    678 @classmethod
    679 def open_group(
    680     cls,
   (...)
    696     cache_members: bool = True,
    697 ):
    698     (
    699         zarr_group,
    700         consolidate_on_close,
    701         close_store_on_close,
    702         use_zarr_fill_value_as_mask,
--> 703     ) = _get_open_params(
    704         store=store,
    705         mode=mode,
    706         synchronizer=synchronizer,
    707         group=group,
    708         consolidated=consolidated,
    709         consolidate_on_close=consolidate_on_close,
    710         chunk_store=chunk_store,
    711         storage_options=storage_options,
    712         zarr_version=zarr_version,
    713         use_zarr_fill_value_as_mask=use_zarr_fill_value_as_mask,
    714         zarr_format=zarr_format,
    715     )
    717     return cls(
    718         zarr_group,
    719         mode,
   (...)
    727         cache_members,
    728     )

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/xarray/backends/zarr.py:1761, in _get_open_params(store, mode, synchronizer, group, consolidated, consolidate_on_close, chunk_store, storage_options, zarr_version, use_zarr_fill_value_as_mask, zarr_format)
   1759 if consolidated is None:
   1760     try:
-> 1761         zarr_group = zarr.open_consolidated(store, **open_kwargs)
   1762     except (ValueError, KeyError):
   1763         # ValueError in zarr-python 3.x, KeyError in 2.x.
   1764         try:

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/zarr/api/synchronous.py:212, in open_consolidated(use_consolidated, *args, **kwargs)
    207 def open_consolidated(*args: Any, use_consolidated: Literal[True] = True, **kwargs: Any) -> Group:
    208     """
    209     Alias for :func:`open_group` with ``use_consolidated=True``.
    210     """
    211     return Group(
--> 212         sync(async_api.open_consolidated(*args, use_consolidated=use_consolidated, **kwargs))
    213     )

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/zarr/core/sync.py:142, in sync(coro, loop, timeout)
    139 return_result = next(iter(finished)).result()
    141 if isinstance(return_result, BaseException):
--> 142     raise return_result
    143 else:
    144     return return_result

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/zarr/core/sync.py:98, in _runner(coro)
     93 """
     94 Await a coroutine and return the result of running it. If awaiting the coroutine raises an
     95 exception, the exception will be returned.
     96 """
     97 try:
---> 98     return await coro
     99 except Exception as ex:
    100     return ex

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/zarr/api/asynchronous.py:346, in open_consolidated(use_consolidated, *args, **kwargs)
    341 if use_consolidated is not True:
    342     raise TypeError(
    343         "'use_consolidated' must be 'True' in 'open_consolidated'. Use 'open' with "
    344         "'use_consolidated=False' to bypass consolidated metadata."
    345     )
--> 346 return await open_group(*args, use_consolidated=use_consolidated, **kwargs)

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/zarr/api/asynchronous.py:800, in open_group(store, mode, cache_attrs, synchronizer, path, chunk_store, storage_options, zarr_version, zarr_format, meta_array, attributes, use_consolidated)
    797 if chunk_store is not None:
    798     warnings.warn("chunk_store is not yet implemented", RuntimeWarning, stacklevel=2)
--> 800 store_path = await make_store_path(store, mode=mode, storage_options=storage_options, path=path)
    802 if attributes is None:
    803     attributes = {}

File ~/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/zarr/storage/_common.py:316, in make_store_path(store_like, path, mode, storage_options)
    314     else:
    315         msg = f"Unsupported type for store_like: '{type(store_like).__name__}'"  # type: ignore[unreachable]
--> 316         raise TypeError(msg)
    318     result = await StorePath.open(store, path=path_normalized, mode=mode)
    320 if storage_options and not used_storage_options:

TypeError: Unsupported type for store_like: 'FSMap'

ds_copepod

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 ds_copepod

NameError: name 'ds_copepod' is not defined

Plot

fig = plt.figure(figsize=(12,3))

ax = fig.add_subplot(1,2,1, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('COPEPOD dataset', fontsize=12)
pc=ax.pcolormesh(ds_copepod.lon, ds_copepod.lat, ds_copepod.copepod_C, cmap='Reds',vmin=0,vmax=2,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(1,2,2, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('CESM ${\it Mesozooplankton}$ biomass', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats, ds_annual.mesozooC.mean(dim='z_t_150m'))
pc=ax.pcolormesh(lon, lat, field, cmap='Reds',vmin=0,vmax=2,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.7])
fig.colorbar(pc, cax=cbar_ax,extend='max', label='top 150m/200m mean (mmol m$^{-3}$)');

<>:10: SyntaxWarning: invalid escape sequence '\i'
<>:10: SyntaxWarning: invalid escape sequence '\i'
/tmp/ipykernel_3705/3097243711.py:10: SyntaxWarning: invalid escape sequence '\i'
  ax.set_title('CESM ${\it Mesozooplankton}$ biomass', fontsize=12)
/tmp/ipykernel_3705/3097243711.py:10: SyntaxWarning: invalid escape sequence '\i'
  ax.set_title('CESM ${\it Mesozooplankton}$ biomass', fontsize=12)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[12], line 5
      3 ax = fig.add_subplot(1,2,1, projection=ccrs.Robinson(central_longitude=305.0))
      4 ax.set_title('COPEPOD dataset', fontsize=12)
----> 5 pc=ax.pcolormesh(ds_copepod.lon, ds_copepod.lat, ds_copepod.copepod_C, cmap='Reds',vmin=0,vmax=2,transform=ccrs.PlateCarree())
      6 land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
      7 ax.add_feature(land)

NameError: name 'ds_copepod' is not defined
../_images/9738aabfdc8cedf87e07d9c820b47c072e8a68ecb6bbc32a3f2b1bbd1536e852.png

Making monthly climatology maps to compare to observations

Compare to observation-based GLMM (Generalized Linear Mixed Model) of global mesozooplankton biomass climatology

This data is from Heneghan et al., 2020, which includes the COPEPOD dataset we used previously as well as additional observations, with some pre-processing.

mesozoo_obs_path = 'data/obsglmm_zmeso_vint_200m_monthly_climatology.nc'

ds_copepod_clim = xr.open_dataset(mesozoo_obs_path)
ds_copepod_clim.zmeso200.attrs['units'] = 'mgC m-2'

months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']

Make our CESM data into a monthly climatology

mon_ds = ds.copy()
mon_ds = ds.groupby('time.month').mean('time')

### depth integrate and convert model to mg C/m2
mon_ds['mesozooC_zint'] = ((mon_ds.mesozooC) * 10.).sum(dim='z_t_150m') #in mmol/m2
mon_ds['mesozooC_zint'] = mon_ds['mesozooC_zint'] * 12.011 #convert to mgC/m2
mon_ds['mesozooC_zint'].attrs['units'] = 'mgC m-2'

Plot

fig = plt.figure(figsize=(5,18))

for row in np.arange(1,13):
    
    ts=row-1
    
    plot = row*2 - 1
    ax = fig.add_subplot(12,2,plot, projection=ccrs.Robinson(central_longitude=305.0))
    ax.set_title(months[ts]+' obs', fontsize=12)
    pc=ax.pcolormesh(ds_copepod_clim.Lon, ds_copepod_clim.Lat, ds_copepod_clim.zmeso200.isel(month=ts), 
                     cmap='Reds',vmin=0,vmax=4000,transform=ccrs.PlateCarree())
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
    ax.add_feature(land)
    
    plot = row*2
    ax = fig.add_subplot(12,2,plot, projection=ccrs.Robinson(central_longitude=305.0))
    ax.set_title(months[ts]+' CESM', fontsize=12)
    tmp = mon_ds.mesozooC_zint.isel(month=ts)
    lon, lat, field = adjust_pop_grid(lons, lats,  tmp)
    pc=ax.pcolormesh(lon, lat, field, cmap='Reds',vmin=0,vmax=4000,transform=ccrs.PlateCarree())
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
    ax.add_feature(land)

cbar_ax = fig.add_axes([0.92, 0.15, 0.03, 0.7])
fig.colorbar(pc, cax=cbar_ax,extend='max', label='Depth-integrated copepod biomass (mg m$^{-2}$)');
../_images/8718bb924c3a05a9f096e68b5fbab2b9346ca5c096dd5a54aca027c29af9210a.png

And close the Dask cluster we spun up at the beginning.

cluster.close()

Summary

You’ve learned how to evaluate zooplankton biomass modeled by CESM-MARBL and compare it to observations.

Resources and references

On this page