Accessing CEFI Cloud data with Python

Accessing CEFI Cloud data with Python#

Packages used#

To use python to access the ERDDAP server directly from your python script or jupyter-notebook, you will need

Xarray
cf_xarray
netcdf4
matplotlib
zarr
fsspec
s3sf
cartopy

Packages

The package netcdf4 developed by UNIDATA and the zarr engine are used at runtime. It is the essential package that support netCDF format output from Xarray and OPeNDAP access be available in the environment when using xarray. The package cartopy is also not needed to read the data, but is used later in the script for making the map plots.

On this page, we will explore the process of extracting a subset of model simualtion data produce from the regional MOM6 model. The model output is to support the Climate Ecosystem and Fishery Initiative. We will showcase how to utilize xarray and fsspec for accessing the data and visualize it on an detailed map. The currently available data can be viewed on

The contents of this folder encompass hindcast simulations derived from the regional MOM6 model spanning the years 1993 to 2019.

Import python packages#

import xarray as xr
import cf_xarray # This is a wrapper for xarray that allows access to the data through standard coordinate names
import fsspec

Info

Thanks to the xarray, fsspec, zarr and kerchunk, accessing data from the cloud directly from Xarray is seamless. For detailed usage guidance, the Xarray documentation offers excellent examples and explanations.

Quick view of the variable#

To check if the data makes sense and fits the chosen time and region, you can use Xarray’s .plot() method. This way, you see a graph, making it easier to understand.

ds_tob_subset.tob.plot(cmap='inferno')

<matplotlib.collections.QuadMesh at 0x7eff29d2b100>

../../../_images/f0678e0e067458d317e66f5878418c53eea762405f98c295176b4607b195ea48.png

Plotting the data on a detailed map#

In this section, we are going to use cartopy to plot the data with a detailed map backgournd.

import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs

Specify figure setup#

Use the cf wrapper to reference the latitude, longitude and time via their CF standard names.
Calculate the aspect ratio of the data to set the aspect ratio of the figure
Shows how to center the plot at 180, not applicable to this data set, but critical for others like the NE Pacific.

xmin = ds_tob_subset.cf['longitude'].min()
xmax = ds_tob_subset.cf['longitude'].max()
ymin = ds_tob_subset.cf['latitude'].min()
ymax = ds_tob_subset.cf['latitude'].max()
aspect = (xmax-xmin)/(ymax-ymin)

Add details to the map#

Use the default land feature to fill in the land mass
Draw the coasts with the 50m data (will trigger a download of the data the first time you use it).
Draw some graticules
Put the colorbar on top and size it to fit the plot
Create a title from the metadata in the file and the variable name
Draw the plot

plt.figure(figsize=(8*aspect,8))
proj = ccrs.PlateCarree(central_longitude=180)
proj180 = ccrs.PlateCarree()
ax = plt.axes(projection=proj)
ax.set_extent([xmin, xmax, ymin, ymax], crs=proj180)
# add some features to make the map a little more polished
ax.add_feature(cartopy.feature.LAND)
ax.coastlines('50m')
gl = ax.gridlines(draw_labels=True)
ct = ax.contourf(ds_tob_subset.cf['longitude'], ds_tob_subset.cf['latitude'], ds_tob_subset['tob'][0], levels=255, transform=proj180, cmap="inferno")
plt.colorbar(ct, orientation='horizontal',pad=0.08, aspect=35, fraction=.06, location='top')
plt.title(str(ds_tob_subset.attrs['title']) + '\ntob\n at t=' + str(ds_tob_subset['time'].values), y=1.25)
plt.show()

../../../_images/13f67d1e84f6a1dda55960e2347164c43d61ad82e461fa014f99e2161a07e62c.png

Quick view of the variable#

To check if the data makes sense and fits the chosen time and region, you can use Xarray’s .plot() method. This way, you see a graph instead of just numbers, making it easier to understand.

ds.tob.plot(x='longitude',y='latitude', cmap='inferno')

<matplotlib.collections.QuadMesh at 0x7eff1631e470>

../../../_images/f45182d71419fab1702b511e5b8b450ec59701ab72be8f2278ff3f981a5e08bd.png

Make a detailed plot showing the underlying grid#

long_name = ds.tob.attrs['long_name']
xmin = ds['longitude'].min()
xmax = ds['longitude'].max()
ymin = ds['latitude'].min()
ymax = ds['latitude'].max()
aspect = (xmax-xmin)/(ymax-ymin)
plt.figure(figsize=(8*aspect,8))
proj = ccrs.PlateCarree(central_longitude=180)
proj180 = ccrs.PlateCarree()
ax = plt.axes(projection=proj)
ax.set_extent([xmin, xmax, ymin, ymax], crs=proj180)
# add some features to make the map a little more polished
ax.add_feature(cartopy.feature.LAND)
ax.coastlines('50m')
gl = ax.gridlines(draw_labels=True)
ct = ax.contourf(ds['longitude'], ds['latitude'], ds['tob'][0], levels=255, transform=proj180, cmap="inferno")
plt.colorbar(ct, orientation='horizontal',pad=0.08, aspect=35, fraction=.06, location='top', cmap="inferno")
ax.autoscale(False) # keep the scalling info the same for the grid dots
ct = ax.scatter(ds['longitude'][::15], ds['latitude'][::15], .05, zorder=1, transform=proj180, color='grey')
plt.title(str(ds.attrs['title']) + '\n+long_name+\n at t=' + str(ds['time'].values)+'\nwith every 15th grid point overlaid', y=1.25)
plt.show()

../../../_images/34cfe6930d0d0d08057573d1c46b66db6ea82850fec7195994a810ec52e092f8.png

Accessing CEFI Cloud data with Python

Contents

Accessing CEFI Cloud data with Python#

Packages used#

Import python packages#

Access data (regular grid product)#

Quick view of the variable#

Plotting the data on a detailed map#

Specify figure setup#

Add details to the map#

Access data (raw grid product)#

Quick view of the variable#

Make a detailed plot showing the underlying grid#