This tutorial will show the steps to grab data in ERDDAP from Python, how to work with NetCDF files in Python and how to make some maps and time-series of sea surface temperature (SST) around the main Hawaiian islands.

This tutorial is also available as a Jupyter notebook.

You will need to install the urllib, xarray and NetCDF4 packages if you don't already have them.

Downloading data in Python

Because ERDDAP includes RESTful services, you can download data listed on any ERDDAP platform from Python using the URL structure.

For example, the following page allows you to subset monthly SST data:

Select your region and date range of interest, then select the '.nc' (NetCDF) file type and click on "Just Generate the URL".

In this specific example, the URL we generated is :

https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v1_0_monthly.nc?analysed_sst[(2018-01-01T12:00:00Z):1:(2018-12-01T12:00:00Z)][(17):1:(30)][(195):1:(210)]

You can also edit this URL manually.

In Python, run the following to download the data using the generated URL :

import urllib.request

url="https://oceanwatch.pifsc.noaa.gov/erddap/griddap/CRW_sst_v1_0_monthly.nc?analysed_sst[(2018-01-01T12:00:00Z):1:(2018-12-01T12:00:00Z)][(17):1:(30)][(195):1:(210)]"

urllib.request.urlretrieve(url, "sst.nc")

Importing the downloaded data in Python

Now that we've downloaded the data locally, we can import it and extract our variables of interest:

The xarray package makes it very convenient to work with NetCDF files. Documentation is available here:
http://xarray.pydata.org/en/stable/why-xarray.html

import xarray as xr import netCDF4 as nc

Open the file and load it as an xarray dataset:

ds = xr.open_dataset('sst.nc',decode_cf=False)

examine the data structure:

ds

Dimensions: (latitude: 261, longitude: 301, time: 12) 
Coordinates: 
    *time (time) float64 1.515e+09 1.517e+09 ... 1.541e+09 1.544e+09 
    *latitude (latitude) float32 17.025 17.075 17.125 ... 29.975 30.025 
    *longitude (longitude) float32 195.025 195.075 ... 209.975 210.025 
Data variables: analysed_sst (time, latitude, longitude) float64 ... 
Attributes: 
acknowledgement: NOAA Coral Reef Watch Program 
cdm_data_type: Grid 
comment: This product is designed to improve on and re... 
contributor_name: NOAA Coral Reef Watch program 
contributor_role: Collecting source data and deriving products;... 
Conventions: CF-1.6, ACDD-1.3, COARDS 
creator_email: coralreefwatch@noaa.gov 
creator_institution: NOAA/NESDIS/STAR Coral Reef Watch program 
creator_name: NOAA Coral Reef Watch program 
creator_type: group 
creator_url: https://coralreefwatch.noaa.gov/ 
data_source: NOAA Daily Global 5km Geo-Polar Blended Night... 
date_created: 2018-01-01T00:00:00Z 
date_issued: 2018-12-02T15:20:07Z 
date_metadata_modified: 2018-09-01T00:00:00Z 
date_modified: 2018-01-01T00:00:00Z 
Easternmost_Easting: 210.025 
geospatial_bounds: "POLYGON((-90.0 360.0, 90.0 360.0, 90.0 0.0, ... 
geospatial_bounds_crs: EPSG:32663 
geospatial_lat_max: 30.025 
geospatial_lat_min: 17.025 
geospatial_lat_resolution: 0.049999999999999996 
geospatial_lat_units: degrees_north 
geospatial_lon_max: 210.025 
geospatial_lon_min: 195.025 
geospatial_lon_resolution: 0.05000000000000001 
geospatial_lon_units: degrees_east 
history: Mon Mar 2 06:00:23 2020: ncatted -O -a geosp... 
id: CoralTemp-v1.0 infoUrl: https://coralreefwatch.noaa.gov/satellite/ble... 
institution: NOAA/NESDIS/STAR Coral Reef Watch program 
instrument: ATSR-1, ATSR-2, AATSR, AVHRR, AVHRR-2, AVHRR-... 
instrument_vocabulary: NOAA NODC Ocean Archive System Instruments 
keywords: 5km, analysed, analysed_sst, analysis, blende... 
keywords_vocabulary: GCMD Science Keywords 
license: OSTIA Usage Statement (1985-2002): IMPORTANT ... 
metadata_link: https://coralreefwatch.noaa.gov/satellite/ble... 
naming_authority: gov.noaa.coralreefwatch 
NCO: 4.3.7 
nco_openmp_thread_number: 1 
Northernmost_Northing: 30.025 
platform: Ships, drifting buoys, moored buoys, TOGA-TAO... 
platform_vocabulary: NOAA NODC Ocean Archive System Platforms 
processing_level: L4 product_version: 1.0 
program: NOAA Coral Reef Watch program 
project: NOAA Coral Reef Watch program 
publisher_email: coralreefwatch@noaa.gov 
publisher_institution: NOAA/NESDIS/STAR Coral Reef Watch program 
publisher_name: NOAA Coral Reef Watch program 
publisher_type: group publisher_
url: https://coralreefwatch.noaa.gov/ 
references: Donlon, et al., 2011. The Operational Sea Sur... 
source: OSTIA Sea Surface Temperature Reanalysis (nig... 
sourceUrl: (local files) 
Southernmost_Northing: 17.025 
standard_name_vocabulary: CF Standard Name Table v27 
summary: CoralTemp 5km gap-free analysed blended sea s... 
time_coverage_duration: P1D 
time_coverage_end: 2018-12-01T12:00:00Z 
time_coverage_resolution: P1D 
time_coverage_start: 2018-01-01T12:00:00Z 
title: Sea Surface Temperature, Coral Reef Watch, Co... 
Westernmost_Easting: 195.025

examine which coordinates and variables are included in the dataset:

ds.coords

Coordinates: 
*time (time) float64 1.515e+09 1.517e+09 ... 1.541e+09 1.544e+09 
*latitude (latitude) float32 17.025 17.075 17.125 ... 29.925 29.975 30.025 
*longitude (longitude) float32 195.025 195.075 195.125 ... 209.975 210.025

ds.data_vars

Data variables: analysed_sst (time, latitude, longitude) float64 ...

examine the structure of analysed_sst:

ds.analysed_sst.shape

(12, 261, 301)

Our dataset is a 3-D array with 261 rows corresponding to latitudes and 301 columns corresponding to longitudes, for each of the 12 time steps.

get the dates for each time step:

dates=nc.num2date(ds.time,ds.time.units,only_use_cftime_datetimes=False) dates

array([datetime.datetime(2018, 1, 1, 12, 0), datetime.datetime(2018, 2, 1, 12, 0), datetime.datetime(2018, 3, 1, 12, 0), datetime.datetime(2018, 4, 1, 12, 0), datetime.datetime(2018, 5, 1, 12, 0), datetime.datetime(2018, 6, 1, 12, 0), datetime.datetime(2018, 7, 1, 12, 0), datetime.datetime(2018, 8, 1, 12, 0), datetime.datetime(2018, 9, 1, 12, 0), datetime.datetime(2018, 10, 1, 12, 0), datetime.datetime(2018, 11, 1, 12, 0), datetime.datetime(2018, 12, 1, 12, 0)], dtype=object)

Working with the extracted data

Creating a map for one time step

Let's create a map of SST for January 2018 (our first time step).

import pandas as pd 
import numpy as np 
from matplotlib import pyplot as plt 
from matplotlib.colors import LinearSegmentedColormap 
np.warnings.filterwarnings('ignore')

set some color breaks

np.nanmin(ds.analysed_sst)

17.922142857142862

np.nanmax(ds.analysed_sst)

28.390645161290323

levs = np.arange(17.5, 28.5, 0.05)

define a color palette

jet=["blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"]

set color scale using the jet palette

cm = LinearSegmentedColormap.from_list('my_jet', jet, N=len(levs))

plot the SST map

plt.contourf(ds.longitude, ds.latitude, ds.analysed_sst[0,:,:], levs,cmap=cm)

plot color scale

plt.colorbar()

example of how to add points to the map

plt.scatter(range(202,206),np.repeat(26,4),c='black')

example of how to add a contour line

plt.contour(ds.longitude, ds.latitude, ds.analysed_sst[0,:,:],levels=[20],linewidths=1)

plot title

plt.title("Monthly Sea Surface Temperature " + dates[0].strftime('%b %Y')) 

plt.show()

Plotting a time series

Let's pick the following box : 18-23N, 200-206E. We are going to generate a time series of mean SST within that box.

first, let subset our data:

da=ds.sel(latitude=slice(18,23), longitude=slice(200,206))

let's plot the subset:

plt.contourf(da.longitude, da.latitude, da.analysed_sst[0,:,:], levs,cmap=cm)
plt.colorbar()
plt.title("Monthly Sea Surface Temperature " + dates[0].strftime('%b %Y'))
plt.show()

let's compute the monthly mean over the bounding region:

res=np.mean(da.analysed_sst,axis=(1,2))

let's plot the time-series:

plt.figure(figsize=(8,4)) 
plt.scatter(dates,res) 
plt.ylabel('SST (ºC)')

Creating a map of average SST over a year

let's compute the yearly mean for the region:

mean_sst=np.mean(ds.analysed_sst,axis=0)

let's plot the map of the 2018 average SST in the region:

plt.contourf(ds.longitude, ds.latitude, mean_sst, levs,cmap=cm) 
plt.colorbar() 
plt.title("Mean SST " + dates[0].strftime('%Y-%m')+' - '+dates[11].strftime('%Y-%m')) 
plt.show()