index

Python Tutorial - How to work with OceanWatch data in Python

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 od chlorophyll-a concentration around the main Hawaiian islands

1. Downlading data from Python

Because ERDDAP includes RESTful services, you can download data listed on any ERDDAP platform from R using the URL structure. For example, the following page allows you to subset monthly Chlorophyll a data from the Aqua-MODIS sensor https://oceanwatch.pifsc.noaa.gov/erddap/griddap/OceanWatch_aqua_chla_monthly.html. Select your region and date range of interest, then select the '.nc' (NetCDF) file type and click on "Just Generate the URL".
In Python, run the following to download the data using the generated URL :
1
import urllib.request
2
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)]"
3
urllib.request.urlretrieve(url, "sst.nc")
Copied!
1
('sst.nc', <http.client.HTTPMessage at 0x2a195863be0>)
Copied!

2. Importing NetCDF4 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
1
import xarray as xr
2
import netCDF4 as nc
Copied!
  • Open the file and load it as an xarray dataset:
1
ds = xr.open_dataset('sst.nc',decode_cf=False)
Copied!
  • examine the data structure:
1
ds
Copied!
1
<xarray.Dataset>
2
Dimensions: (latitude: 261, longitude: 301, time: 12)
3
Coordinates:
4
* time (time) float64 1.515e+09 1.517e+09 ... 1.541e+09 1.544e+09
5
* latitude (latitude) float32 17.025 17.075 17.125 ... 29.975 30.025
6
* longitude (longitude) float32 195.025 195.075 ... 209.975 210.025
7
Data variables:
8
analysed_sst (time, latitude, longitude) float64 ...
9
Attributes:
10
acknowledgement: NOAA Coral Reef Watch Program
11
cdm_data_type: Grid
12
comment: This product is designed to improve on and re...
13
contributor_name: NOAA Coral Reef Watch program
14
contributor_role: Collecting source data and deriving products;...
15
Conventions: CF-1.6, ACDD-1.3, COARDS
16
creator_email: [email protected]
17
creator_institution: NOAA/NESDIS/STAR Coral Reef Watch program
18
creator_name: NOAA Coral Reef Watch program
19
creator_type: group
20
creator_url: https://coralreefwatch.noaa.gov/
21
data_source: NOAA Daily Global 5km Geo-Polar Blended Night...
22
date_created: 2018-01-01T00:00:00Z
23
date_issued: 2018-12-02T15:20:07Z
24
date_metadata_modified: 2018-09-01T00:00:00Z
25
date_modified: 2018-01-01T00:00:00Z
26
Easternmost_Easting: 210.025
27
geospatial_bounds: "POLYGON((-90.0 360.0, 90.0 360.0, 90.0 0.0, ...
28
geospatial_bounds_crs: EPSG:32663
29
geospatial_lat_max: 30.025
30
geospatial_lat_min: 17.025
31
geospatial_lat_resolution: 0.049999999999999996
32
geospatial_lat_units: degrees_north
33
geospatial_lon_max: 210.025
34
geospatial_lon_min: 195.025
35
geospatial_lon_resolution: 0.05000000000000001
36
geospatial_lon_units: degrees_east
37
history: Mon Mar 2 06:00:23 2020: ncatted -O -a geosp...
38
id: CoralTemp-v1.0
39
infoUrl: https://coralreefwatch.noaa.gov/satellite/ble...
40
institution: NOAA/NESDIS/STAR Coral Reef Watch program
41
instrument: ATSR-1, ATSR-2, AATSR, AVHRR, AVHRR-2, AVHRR-...
42
instrument_vocabulary: NOAA NODC Ocean Archive System Instruments
43
keywords: 5km, analysed, analysed_sst, analysis, blende...
44
keywords_vocabulary: GCMD Science Keywords
45
license: OSTIA Usage Statement (1985-2002): IMPORTANT ...
46
metadata_link: https://coralreefwatch.noaa.gov/satellite/ble...
47
naming_authority: gov.noaa.coralreefwatch
48
NCO: 4.3.7
49
nco_openmp_thread_number: 1
50
Northernmost_Northing: 30.025
51
platform: Ships, drifting buoys, moored buoys, TOGA-TAO...
52
platform_vocabulary: NOAA NODC Ocean Archive System Platforms
53
processing_level: L4
54
product_version: 1.0
55
program: NOAA Coral Reef Watch program
56
project: NOAA Coral Reef Watch program
57
publisher_email: [email protected]
58
publisher_institution: NOAA/NESDIS/STAR Coral Reef Watch program
59
publisher_name: NOAA Coral Reef Watch program
60
publisher_type: group
61
publisher_url: https://coralreefwatch.noaa.gov/
62
references: Donlon, et al., 2011. The Operational Sea Sur...
63
source: OSTIA Sea Surface Temperature Reanalysis (nig...
64
sourceUrl: (local files)
65
Southernmost_Northing: 17.025
66
standard_name_vocabulary: CF Standard Name Table v27
67
summary: CoralTemp 5km gap-free analysed blended sea s...
68
time_coverage_duration: P1D
69
time_coverage_end: 2018-12-01T12:00:00Z
70
time_coverage_resolution: P1D
71
time_coverage_start: 2018-01-01T12:00:00Z
72
title: Sea Surface Temperature, Coral Reef Watch, Co...
73
Westernmost_Easting: 195.025
Copied!
  • examine which coordinates and variables are included in the dataset:
1
ds.coords
Copied!
1
Coordinates:
2
* time (time) float64 1.515e+09 1.517e+09 ... 1.541e+09 1.544e+09
3
* latitude (latitude) float32 17.025 17.075 17.125 ... 29.925 29.975 30.025
4
* longitude (longitude) float32 195.025 195.075 195.125 ... 209.975 210.025
Copied!
1
ds.data_vars
Copied!
1
Data variables:
2
analysed_sst (time, latitude, longitude) float64 ...
Copied!
  • examine the structure of analysed_sst:
1
ds.analysed_sst.shape
Copied!
1
(12, 261, 301)
Copied!
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:
1
ds.time
Copied!
1
<xarray.DataArray 'time' (time: 12)>
2
array([1.514808e+09, 1.517486e+09, 1.519906e+09, 1.522584e+09, 1.525176e+09,
3
1.527854e+09, 1.530446e+09, 1.533125e+09, 1.535803e+09, 1.538395e+09,
4
1.541074e+09, 1.543666e+09])
5
Coordinates:
6
* time (time) float64 1.515e+09 1.517e+09 1.52e+09 ... 1.541e+09 1.544e+09
7
Attributes:
8
_CoordinateAxisType: Time
9
actual_range: [1.5148080e+09 1.5436656e+09]
10
axis: T
11
coverage_content_type: coordinate
12
ioos_category: Time
13
long_name: reference time of the sst field
14
standard_name: time
15
time_origin: 01-JAN-1970 00:00:00
16
units: seconds since 1970-01-01T00:00:00Z
Copied!
1
dates=nc.num2date(ds.time,ds.time.units)
2
dates
Copied!
1
array([datetime.datetime(2018, 1, 1, 12, 0),
2
datetime.datetime(2018, 2, 1, 12, 0),
3
datetime.datetime(2018, 3, 1, 12, 0),
4
datetime.datetime(2018, 4, 1, 12, 0),
5
datetime.datetime(2018, 5, 1, 12, 0),
6
datetime.datetime(2018, 6, 1, 12, 0),
7
datetime.datetime(2018, 7, 1, 12, 0),
8
datetime.datetime(2018, 8, 1, 12, 0),
9
datetime.datetime(2018, 9, 1, 12, 0),
10
datetime.datetime(2018, 10, 1, 12, 0),
11
datetime.datetime(2018, 11, 1, 12, 0),
12
datetime.datetime(2018, 12, 1, 12, 0)], dtype=object)
Copied!

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).
1
import pandas as pd
2
import numpy as np
3
from matplotlib import pyplot as plt
4
from matplotlib.colors import LinearSegmentedColormap
5
np.warnings.filterwarnings('ignore')
Copied!
  • set some color breaks
1
np.nanmin(ds.analysed_sst)
Copied!
1
17.922142857142862
Copied!
1
np.nanmax(ds.analysed_sst)
Copied!
1
28.390645161290323
Copied!
1
levs = np.arange(17.5, 28.5, 0.05)
Copied!
  • define a color palette
1
jet=["blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"]
Copied!
  • set color scale using the jet palette
1
cm = LinearSegmentedColormap.from_list('my_jet', jet, N=len(levs))
Copied!
  • plot the SST map
1
plt.contourf(ds.longitude, ds.latitude, ds.analysed_sst[0,:,:], levs,cmap=cm)
2
#plot the color scale
3
plt.colorbar()
4
#example of how to add points to the map
5
plt.scatter(range(202,206),np.repeat(26,4),c='black')
6
#example of how to add a contour line
7
plt.contour(ds.longitude, ds.latitude, ds.analysed_sst[0,:,:],levels=20,linewidths=1)
8
#plot title
9
plt.title("Monthly Sea Surface Temperature " + dates[0].strftime('%b %Y'))
10
plt.show()
Copied!

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:
1
lat_bnds, lon_bnds = [18, 23], [200, 206]
2
da=ds.sel(latitude=slice(*lat_bnds), longitude=slice(*lon_bnds))
Copied!
  • let's plot the subset:
1
plt.contourf(da.longitude, da.latitude, da.analysed_sst[0,:,:], levs,cmap=cm)
2
plt.colorbar()
3
plt.title("Monthly Sea Surface Temperature " + dates[0].strftime('%b %Y'))
4
plt.show()
Copied!
  • let's compute the monthly mean over the bounding region:
1
res=np.mean(da.analysed_sst,axis=(1,2))
Copied!
  • let's plot the time-series:
1
plt.figure(figsize=(8,4))
2
plt.scatter(dates,res)
3
plt.ylabel('SST (ºC)')
Copied!
1
Text(0,0.5,'SST (ºC)')
Copied!

Creating a map of average SST over a year

  • let's compute the yearly mean for the region:
1
mean_sst=np.mean(ds.analysed_sst,axis=0)
Copied!
1
mean_sst.shape
Copied!
1
(261, 301)
Copied!
  • let's plot the map of the 2018 average SST in the region:
1
plt.contourf(ds.longitude, ds.latitude, mean_sst, levs,cmap=cm)
2
plt.colorbar()
3
plt.title("Mean SST " + dates[0].strftime('%Y-%m')+' - '+dates[11].strftime('%Y-%m'))
4
plt.show()
Copied!
Last modified 4mo ago