Example calculations with scalar quantities

Objectives

To demonstrate basic calculations using scalar fields (e.g., SSH, T, S) from the state estimate including: time series of mean quantities, spatial patterns of mean quantities, spatial patterns of linear trends, and spatial patterns of linear trends over different time periods.

Introduction

We will demonstrate global calculations with SSH (global mean sea level time series, mean dynamic topography, global mean sea level trend) and a regional calculation with THETA (The Nino 3.4 index).

For this tutorial we will need the grid file, monthly SSH and THETA (potential temperature) spanning the years 1993 through 2017, and daily SSH for the year 1994. The ShortNames of the datasets are:

  • ECCO_L4_GEOMETRY_LLC0090GRID_V4R4

  • ECCO_L4_SSH_LLC0090GRID_MONTHLY_V4R4 (1993-2017)

  • ECCO_L4_TEMP_SALINITY_LLC0090GRID_MONTHLY_V4R4 (1993-2017)

  • ECCO_L4_SSH_LLC0090GRID_DAILY_V4R4 (1994)

If you have set up the ecco_access library in your Python path, the notebook will access the data either by download or in-cloud direct access from S3, depending on the access mode you specify below.

Global calculations with SSH

First, load SSH and THETA fields and the model grid parameters.

[1]:

import numpy as np
import sys
import xarray as xr
from copy import deepcopy
import matplotlib.pyplot as plt
%matplotlib inline
import glob
import warnings
warnings.filterwarnings('ignore')

import ecco_access as ea

# indicate mode of access
# options are:
# 'download': direct download from internet to your local machine
# 'download_ifspace': like download, but only proceeds
#                     if your machine have sufficient storage
# 's3_open': access datasets in-cloud from an AWS instance
# 's3_open_fsspec': use jsons generated with fsspec and
#                   kerchunk libraries to speed up in-cloud access
# 's3_get': direct download from S3 in-cloud to an AWS instance
# 's3_get_ifspace': like s3_get, but only proceeds if your instance
#                   has sufficient storage
access_mode = 'download_ifspace'

[2]:

## Import the ecco_v4_py library into Python
## =========================================
##    If ecco_v4_py is not installed in your local Python library,
##    tell Python where to find it.  The example below adds
##    ecco_v4_py to the user's path if it is stored in the folder
##    ECCOv4-py under the user's home directory

from os.path import join,expanduser
user_home_dir = expanduser('~')

sys.path.append(join(user_home_dir,'ECCOv4-py'))

import ecco_v4_py as ecco

[3]:

## Set top-level file directory for the ECCO NetCDF files
## =================================================================

## currently set to ~/Downloads/ECCO_V4r4_PODAAC
ECCO_dir = join(user_home_dir,'Downloads','ECCO_V4r4_PODAAC')

# # for access_mode = 's3_open_fsspec', need to specify the root directory
# # containing the jsons
# jsons_root_dir = join('/efs_ecco','mzz-jsons')

[4]:

## Access datasets needed for this tutorial

ShortNames_list = ["ECCO_L4_GEOMETRY_LLC0090GRID_V4R4",\
                   "ECCO_L4_SSH_LLC0090GRID_MONTHLY_V4R4",\
                   "ECCO_L4_TEMP_SALINITY_LLC0090GRID_MONTHLY_V4R4"]
ShortNames_daily_list = ["ECCO_L4_SSH_LLC0090GRID_DAILY_V4R4"]

ds_dict = ea.ecco_podaac_to_xrdataset(ShortNames_list,\
                                              StartDate='1993-01',EndDate='2017-12',\
                                              mode=access_mode,\
                                              download_root_dir=ECCO_dir,\
                                              max_avail_frac=0.5)
ds_SSH_daily_1994 = ea.ecco_podaac_to_xrdataset(ShortNames_daily_list,\
                                                    StartDate='1994-01-01',EndDate='1994-12-31',\
                                                    mode=access_mode,\
                                                    download_root_dir=ECCO_dir,\
                                                    max_avail_frac=0.5)

Size of files to be downloaded to instance is 6.316 GB,
which is 4.4% of the 143.655 GB available storage.
Proceeding with file downloads via NASA Earthdata URLs

GRID_GEOMETRY_ECCO_V4r4_native_llc0090.nc already exists, and force=False, not re-downloading
DL Progress: 100%|########################| 1/1 [00:00<00:00, 12483.05it/s]

=====================================
total downloaded: 0.0 Mb
avg download speed: 0.0 Mb/s
Time spent = 0.0030794143676757812 seconds


DL Progress: 100%|#######################| 300/300 [01:38<00:00,  3.04it/s]

=====================================
total downloaded: 1775.35 Mb
avg download speed: 17.98 Mb/s
Time spent = 98.72319316864014 seconds


DL Progress:  100%|#######################| 300/300 [01:26<01:26,  2.59it/s]

=====================================
total downloaded: 1808.94 Mb
avg download speed: 8.23 Mb/s
Time spent = 219.87921047210693 seconds


Size of files to be downloaded to instance is 2.015 GB,
which is 1.47% of the 137.337 GB available storage.
Proceeding with file downloads via NASA Earthdata URLs
DL Progress: 100%|#######################| 365/365 [02:00<00:00,  3.03it/s]

=====================================
total downloaded: 2163.84 Mb
avg download speed: 17.94 Mb/s
Time spent = 120.64634680747986 seconds



[5]:

## Load the model grid
ecco_grid = ds_dict[ShortNames_list[0]].compute()

Now we’ll look at a dataset with grid parameters and daily SSH during 1994.

[6]:

%%time

ds_SSH_daily_1994 = ds_SSH_daily_1994.drop(['SSHNOIBC','SSHIBC','ETAN'])
ecco_daily_ds = xr.merge((ecco_grid,ds_SSH_daily_1994))
ecco_daily_ds

CPU times: user 87.7 ms, sys: 4.29 ms, total: 92 ms
Wall time: 315 ms

[6]:

<xarray.Dataset> Size: 242MB
Dimensions:    (tile: 13, j: 90, i: 90, k: 50, k_p1: 51, nb: 4, j_g: 90,
                i_g: 90, nv: 2, k_l: 50, k_u: 50, time: 365)
Coordinates: (12/22)
    XC         (tile, j, i) float32 421kB dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    XC_bnds    (tile, j, i, nb) float32 2MB dask.array<chunksize=(13, 90, 90, 4), meta=np.ndarray>
    XG         (tile, j_g, i_g) float32 421kB dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    YC         (tile, j, i) float32 421kB dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    YC_bnds    (tile, j, i, nb) float32 2MB dask.array<chunksize=(13, 90, 90, 4), meta=np.ndarray>
    YG         (tile, j_g, i_g) float32 421kB dask.array<chunksize=(13, 90, 90), meta=np.ndarray>
    ...         ...
  * k_l        (k_l) int32 200B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48 49
  * k_p1       (k_p1) int32 204B 0 1 2 3 4 5 6 7 8 ... 43 44 45 46 47 48 49 50
  * k_u        (k_u) int32 200B 0 1 2 3 4 5 6 7 8 ... 41 42 43 44 45 46 47 48 49
  * tile       (tile) int32 52B 0 1 2 3 4 5 6 7 8 9 10 11 12
  * time       (time) datetime64[ns] 3kB 1994-01-01T12:00:00 ... 1994-12-31T1...
    time_bnds  (time, nv) datetime64[ns] 6kB dask.array<chunksize=(365, 2), meta=np.ndarray>
Dimensions without coordinates: nb, nv
Data variables: (12/22)
    CS         (tile, j, i) float32 421kB 0.06158 0.06675 ... -0.9854 -0.9984
    Depth      (tile, j, i) float32 421kB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    PHrefC     (k) float32 200B 49.05 147.1 245.2 ... 5.357e+04 5.794e+04
    PHrefF     (k_p1) float32 204B 0.0 98.1 196.2 ... 5.57e+04 6.018e+04
    SN         (tile, j, i) float32 421kB -0.9981 -0.9978 ... -0.1705 -0.05718
    drC        (k_p1) float32 204B 5.0 10.0 10.0 10.0 ... 422.0 445.0 228.2
    ...         ...
    maskW      (k, tile, j, i_g) bool 5MB False False False ... False False
    rA         (tile, j, i) float32 421kB 3.623e+08 3.633e+08 ... 3.611e+08
    rAs        (tile, j_g, i) float32 421kB 1.802e+08 1.807e+08 ... 3.605e+08
    rAw        (tile, j, i_g) float32 421kB 3.617e+08 3.628e+08 ... 3.648e+08
    rAz        (tile, j_g, i_g) float32 421kB 1.799e+08 1.805e+08 ... 3.642e+08
    SSH        (time, tile, j, i) float32 154MB dask.array<chunksize=(223, 13, 90, 90), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                     CF-1.8, ACDD-1.3
    acknowledgement:                 This research was carried out by the Jet...
    author:                          Ian Fenty and Ou Wang
    cdm_data_type:                   Grid
    comment:                         Fields provided on the curvilinear lat-l...
    coordinates_comment:             Note: the global 'coordinates' attribute...
    ...                              ...
    references:                      ECCO Consortium, Fukumori, I., Wang, O.,...
    source:                          The ECCO V4r4 state estimate was produce...
    standard_name_vocabulary:        NetCDF Climate and Forecast (CF) Metadat...
    summary:                         This dataset provides geometric paramete...
    title:                           ECCO Geometry Parameters for the Lat-Lon...
    uuid:                            87ff7d24-86e5-11eb-9c5f-f8f21e2ee3e0

Now load monthly mean SSH, 1993-2017.

[7]:

%%time

## Load monthly SSH

ds_SSH_monthly = ds_dict[ShortNames_list[1]].drop(['SSHNOIBC','SSHIBC','ETAN'])

CPU times: user 152 µs, sys: 11 µs, total: 163 µs
Wall time: 166 µs

And monthly mean THETA, 1993-2017.

[8]:

%%time

## Load monthly THETA

# concatenate 2 datasets
ds_THETA_monthly = ds_dict[ShortNames_list[2]].drop('SALT')

## Merge the ecco_grid with the data variables to make ecco_monthly_ds
ecco_monthly_ds = xr.merge((ecco_grid,ds_SSH_monthly,ds_THETA_monthly))

CPU times: user 299 ms, sys: 56.5 ms, total: 355 ms
Wall time: 1.09 s

Display the first and last time entries in each dataset:

[9]:

print(ecco_daily_ds.time[0].values)
print(ecco_daily_ds.time[-1].values)

print(ecco_monthly_ds.time[0].values)
print(ecco_monthly_ds.time[-1].values)

1994-01-01T12:00:00.000000000
1994-12-31T12:00:00.000000000
1993-01-16T12:00:00.000000000
2017-12-16T06:00:00.000000000

Sea surface height

Global mean sea level

Global mean sea surface height at time t is defined as follows:

SSH_{\text{global mean}}(t) = \frac{\sum_{i} SSH(i,t) \,\, A(i)}{A_{\text{global ocean}}}

A_{\text{global ocean}} = \sum_{i} A(i)

Where SSH(i,t) is dynamic height at model grid cell i and time t, A(i) is the area (m^2) of model grid cell i

There are several ways of doing the above calculations. Since this is the first tutorial with actual calcuations, we’ll present a few different approaches for getting to the same answer.

Part 1: A_{\text{global ocean}}

Let’s start on the simplest quantity, the global ocean surface area A_{\text{global ocean}}. Our calculation uses SSH which is a ‘c’ point variable. The surface area of tracer grid cells is provided by the model grid parameter rA. rA is a two-dimensional field that is defined over all model grid points, including land.

To calculate the total ocean surface area we need to ignore the area contributions from land.

We will first construct a 3D mask that is True for model grid cells that are wet and False for model grid cells that are dry cells.

[10]:

# ocean_mask is ceiling of hFacC which is 0 for 100% dry cells and
# 0 > hFacC >= 1 for grid cells that are at least partially wet

# hFacC is the fraction of the thickness (h) of the grid cell which
# is wet.  we'll consider all hFacC > 0 as being a wet grid cell
# and so we use the 'ceiling' function to make all values > 0 equal to 1.

ocean_mask = np.ceil(ecco_monthly_ds.hFacC)
ocean_mask = ocean_mask.where(ocean_mask==1, np.nan)

[11]:

# print the dimensions and size of ocean_mask
print(type(ocean_mask))
print((ocean_mask.dims))

<class 'xarray.core.dataarray.DataArray'>
('k', 'tile', 'j', 'i')

[12]:

plt.figure(figsize=(12,6), dpi= 90)

ecco.plot_tiles(ocean_mask.isel(k=0),layout='latlon', rotate_to_latlon=True)

# select out the model depth at k=1, round the number and convert to string.
z = str((np.round(ecco_monthly_ds.Z.values[0])))
plt.suptitle('Wet (1) /dry (0) mask for k=' + str(0) + ',   z=' + z + 'm');

<Figure size 1080x540 with 0 Axes>
_images/ECCO_v4_Example_calculations_with_scalar_quantities_17_1.png

To calculate A_{\text{global ocean}} we must apply the surface wet/dry mask to rA.

[13]:

# Method 1: the array index method, []
#           select land_c at k index 0
total_ocean_area = np.sum(ecco_monthly_ds.rA*ocean_mask[0,:])

# these three methods give the same numerical result.  Here are
# three alternative ways of printing the result
print ('total ocean surface area ( m^2) %d  ' % total_ocean_area.values)
print ('total ocean surface area (km^2) %d ' % (total_ocean_area.values/1.0e6))

# or in scientific notation with 2 decimal points
print ('total ocean surface area (km^2) %.2E' % (total_ocean_area.values/1.0e6))

total ocean surface area ( m^2) 358013844062208
total ocean surface area (km^2) 358013844
total ocean surface area (km^2) 3.58E+08

This compares favorably with approx 3.60 x 10^8 km^2 from https://hypertextbook.com/facts/1997/EricCheng.shtml

Multiplication of DataArrays

You probably noticed that the multiplication of grid cell area with the land mask was done element by element. One useful feature of DataArrays is that their dimensions are automatically lined up when doing binary operations. Also, because rA and ocean_mask are both DataArrays, their inner product and their sums are also DataArrays.

Note:: ocean_mask has a depth (k) dimension while rA does not (horizontal model grid cell area does not change as a function of depth in ECCOv4). As a result, when rA is multiplied with ocean_mask xarray broadcasts rA to all k levels. The resulting matrix inherits the k dimension from ocean_mask.

Another way of summing over numpy arrays

As rA and ocean both store numpy arrays, you can also calculate the sum of their product by invoking the .sum() command inherited in all numpy arrays:

[14]:

total_ocean_area = (ecco_monthly_ds.rA*ocean_mask).isel(k=0).sum()
print ('total ocean surface area (km^2) ' + '%.2E' % (total_ocean_area.values/1e6))

total ocean surface area (km^2) 3.58E+08

Part2 : SSH_{\text{global mean}}(t)

The global mean SSH at each t is given by,

SSH_{\text{global mean}}(t) = \frac{\sum_{i} SSH(i,t) \,\, A(i)}{A_{\text{global ocean}}}

One way of calculating this is to take advantage of DataArray coordinate labels and use its .sum() functionality to explicitly specify which dimensions to sum over:

[15]:

# note no need to multiple RAC by land_c because SSH is nan over land
SSH_global_mean_mon = (ecco_monthly_ds.SSH*ecco_monthly_ds.rA).sum(dim=['i','j','tile'])/total_ocean_area

[16]:

# remove time mean from time series
SSH_global_mean_mon = SSH_global_mean_mon-SSH_global_mean_mon.mean(dim='time')

[17]:

# add helpful unit label
SSH_global_mean_mon.attrs['units']='m'

[18]:

# and plot for fun
SSH_global_mean_mon.plot(color='k');plt.grid()
_images/ECCO_v4_Example_calculations_with_scalar_quantities_26_0.png

Alternatively we can do the summation over the three non-time dimensions. The time dimension of SSH is along the first dimension (axis) of the array, axis 0.

[19]:

# note no need to multiple RAC by land_c because SSH is nan over land
SSH_global_mean = np.sum(ecco_monthly_ds.SSH*ecco_monthly_ds.rA,axis=(1,2,3))/total_ocean_area
SSH_global_mean = SSH_global_mean.compute()

Even though SSH has 3 dimensions (time, tile, j, i) and rA and ocean_mask.isel(k=0) have 2 (j,i), we can mulitply them. With xarray the element-by-element multiplication occurs over their common dimension.

The resulting SSH_{global-mean} DataArray has a single dimension, time.

Part 3 : Plotting the global mean sea level time series:

Before we plot the global mean sea level curve let’s remove its time-mean to make it global mean sea level anomaly (the absolute value has no meaning here anyway).

[20]:

plt.figure(figsize=(8,4), dpi= 90)

# Method 1: .mean() method of `DataArrays`
SSH_global_mean_anomaly = SSH_global_mean - SSH_global_mean.mean()

# Method 2: numpy's `mean` method
SSH_global_mean_anomaly = SSH_global_mean - np.mean(SSH_global_mean)

SSH_global_mean_anomaly.plot()
plt.grid()
plt.title('ECCO v4r4 Global Mean Sea Level Anomaly');
plt.ylabel('m');
plt.xlabel('year');
_images/ECCO_v4_Example_calculations_with_scalar_quantities_30_0.png

Mean Dynamic Topography

Mean dynamic topography is calculated as follows,

MDT(i) = \frac{\sum_{t} SSH(i,t) - SSH_{\text{global mean}}(t)}{nt}

Where nt is the number of time records.

For MDT we preserve the spatial dimensions. Summation and averaging are over the time dimensions (axis 0).

[21]:

## Two equivalent methods

# Method 1, specify the axis over which to average
MDT = np.mean(ecco_monthly_ds.SSH - SSH_global_mean,axis=0)

# Method 2, specify the coordinate label over which to average
MDT_B = (ecco_monthly_ds.SSH - SSH_global_mean).mean(dim=['time'])

# which can be verified using the '.equals()' method to compare Datasets and DataArrays
print(MDT.equals(MDT_B))

True

As expected, MDT has preserved its spatial dimensions:

[22]:

MDT.dims

[22]:

('tile', 'j', 'i')

Before plotting the MDT field remove its spatial mean since its spatial mean conveys no dynamically useful information.

[23]:

MDT_no_spatial_mean = MDT - MDT*ecco_monthly_ds.rA/total_ocean_area

[24]:

MDT_no_spatial_mean.shape

[24]:

(13, 90, 90)

[25]:

plt.figure(figsize=(12,6), dpi= 90)

# mask land points to Nan
MDT_no_spatial_mean = MDT_no_spatial_mean.where(ocean_mask[0,:] !=0)

ecco.plot_proj_to_latlon_grid(ecco_monthly_ds.XC, \
                              ecco_monthly_ds.YC, \
                              MDT_no_spatial_mean*ocean_mask.isel(k=0), \
                              user_lon_0=0,\
                              plot_type = 'pcolormesh', \
                              cmap='seismic',show_colorbar=True,\
                              dx=2,dy=2);

plt.title('ECCO v4r4 Mean Dynamic Topography [m]');
_images/ECCO_v4_Example_calculations_with_scalar_quantities_38_0.png

Constructing Monthly means from Daily means

We can also construct our own monthly means from the daily means using this command: (See http://xarray.pydata.org/en/stable/generated/xarray.Dataset.resample.html for more information)

[32]:

# note no need to multiple RAC by land_c because SSH is nan over land
SSH_global_mean_day = (ecco_daily_ds.SSH*ecco_daily_ds.rA).sum(dim=['i','j','tile'])/total_ocean_area

[33]:

# remove time mean from time series and replace with 1994 time mean from monthly SSH
SSH_global_mean_day = SSH_global_mean_day-SSH_global_mean_day.mean(dim='time')\
                        +SSH_global_mean_mon.sel(time='1994').mean(dim='time')

[34]:

# add helpful unit label
SSH_global_mean_day.attrs['units']='m'

[35]:

# and plot for fun
SSH_global_mean_day.plot(color='k');plt.grid()
_images/ECCO_v4_Example_calculations_with_scalar_quantities_55_0.png 
[36]:

SSH_global_mean_mon_alt = SSH_global_mean_day.resample(time='1M', loffset='-15D').mean()

Plot to compare.

[37]:

SSH_global_mean_mon.sel(time='1994').plot(color='r', marker='.');
SSH_global_mean_mon_alt.sel(time='1994').plot(color='g', marker='o');
plt.grid()
_images/ECCO_v4_Example_calculations_with_scalar_quantities_58_0.png

These small differences are simply an artifact of the time indexing. We used loffset=’15D’ to shift the time of the monthly mean SSH back 15 days, close to the center of the month. The SSH_global_mean_mon field is centered exactly at the middle of the month, and since months aren’t exactly 30 days long, this results in a small discrepancy when plotting with a time x-axis. If we plot without a time object x axis we find the values to be the same. That’s because ECCO monthly means are calculated over calendar months.

[38]:

print ('date in middle of month')
print(SSH_global_mean_mon.time.values[0:2])
print ('\ndate with 15 day offset from the end of the month')
print(SSH_global_mean_mon_alt.time.values[0:2])

date in middle of month
['1993-01-16T12:00:00.000000000' '1993-02-15T00:00:00.000000000']

date with 15 day offset from the end of the month
['1994-01-16T00:00:00.000000000' '1994-02-13T00:00:00.000000000']

[39]:

plt.plot(SSH_global_mean_mon.sel(time='1994').values, color='r', marker='.');
plt.plot(SSH_global_mean_mon_alt.sel(time='1994').values, color='g', marker='o');
plt.xlabel('months since 1994-01');
plt.ylabel('m')
plt.grid()
_images/ECCO_v4_Example_calculations_with_scalar_quantities_61_0.png

Regional calculations with THETA

[40]:

lat_bounds = np.logical_and(ecco_monthly_ds.YC >= -5, ecco_monthly_ds.YC <= 5)
lon_bounds = np.logical_and(ecco_monthly_ds.XC >= -170, ecco_monthly_ds.XC <= -120)

SST = ecco_monthly_ds.THETA.isel(k=0)
SST_masked=SST.where(np.logical_and(lat_bounds, lon_bounds))

[41]:

plt.figure(figsize=(12,5), dpi= 90)

ecco.plot_proj_to_latlon_grid(ecco_monthly_ds.XC, \
                              ecco_monthly_ds.YC, \
                              SST_masked.isel(time=0),\
                              user_lon_0 = -66,\
                              show_colorbar=True);

plt.title('SST in Niño 3.4 box: \n %s ' % str(ecco_monthly_ds.time[0].values));
_images/ECCO_v4_Example_calculations_with_scalar_quantities_64_0.png 
[42]:

# Create the same mask for the grid cell area
rA_masked=ecco_monthly_ds.rA.where(np.logical_and(lat_bounds, lon_bounds));

# Calculate the area-weighted mean in the box
SST_masked_mean=(SST_masked*rA_masked).sum(dim=['tile','j','i'])/np.sum(rA_masked)

# Substract the temporal mean from the area-weighted mean to get a time series, the Nino 3.4 index
SST_nino_34_anom_ECCO_monthly_mean = SST_masked_mean - np.mean(SST_masked_mean)

Load up the Niño 3.4 index values from ESRL

[43]:

# https://psl.noaa.gov/gcos_wgsp/Timeseries/Data/nino34.long.anom.data
# NINO34
# 5N-5S 170W-120W
# HadISST
#  Anomaly from 1981-2010
#  units=degC

import urllib.request
data = urllib.request.urlopen('https://psl.noaa.gov/gcos_wgsp/Timeseries/Data/nino34.long.anom.data')

# the following code parses the ESRL text file and puts monthly-mean nino 3.4 values into an array
start_year = 1993
end_year = 2017
num_years = end_year-start_year+1
nino34_noaa = np.zeros((num_years, 12))
for i,l in enumerate(data):
    line_str = str(l, "utf-8")
    x=line_str.split()
    try:
        year = int(x[0])
        row_i = year-start_year
        if row_i >= 0 and year <= end_year:

            print('loading Niño 3.4 for year %s  row %s' % (year, row_i))

            for m in range(0,12):
                nino34_noaa[row_i, m] = float(x[m+1])
    except:
        continue

loading Niño 3.4 for year 1993  row 0
loading Niño 3.4 for year 1994  row 1
loading Niño 3.4 for year 1995  row 2
loading Niño 3.4 for year 1996  row 3
loading Niño 3.4 for year 1997  row 4
loading Niño 3.4 for year 1998  row 5
loading Niño 3.4 for year 1999  row 6
loading Niño 3.4 for year 2000  row 7
loading Niño 3.4 for year 2001  row 8
loading Niño 3.4 for year 2002  row 9
loading Niño 3.4 for year 2003  row 10
loading Niño 3.4 for year 2004  row 11
loading Niño 3.4 for year 2005  row 12
loading Niño 3.4 for year 2006  row 13
loading Niño 3.4 for year 2007  row 14
loading Niño 3.4 for year 2008  row 15
loading Niño 3.4 for year 2009  row 16
loading Niño 3.4 for year 2010  row 17
loading Niño 3.4 for year 2011  row 18
loading Niño 3.4 for year 2012  row 19
loading Niño 3.4 for year 2013  row 20
loading Niño 3.4 for year 2014  row 21
loading Niño 3.4 for year 2015  row 22
loading Niño 3.4 for year 2016  row 23
loading Niño 3.4 for year 2017  row 24

[44]:

SST_nino_34_anom_ECCO_monthly_mean.plot();plt.grid()
plt.title('Niño3.4 index from ECCO')

[44]:

Text(0.5, 1.0, 'Niño3.4 index from ECCO')
_images/ECCO_v4_Example_calculations_with_scalar_quantities_68_1.png

we’ll make a new DataArray for the NOAA SST nino_34 data by copying the DataArray for the ECCO SST data and replacing the values

[45]:

SST_nino_34_anom_NOAA_monthly_mean = xr.DataArray(data=nino34_noaa.ravel(),\
                                                  dims=SST_nino_34_anom_ECCO_monthly_mean.dims,\
                                                  coords=SST_nino_34_anom_ECCO_monthly_mean.coords)

[46]:

SST_nino_34_anom_NOAA_monthly_mean.plot();plt.grid()
plt.title('Niño3.4 index from NOAA ESRL')

[46]:

Text(0.5, 1.0, 'Niño3.4 index from NOAA ESRL')
_images/ECCO_v4_Example_calculations_with_scalar_quantities_71_1.png

Plot the ECCOv4r4 and ESRL Niño 3.4 index

[47]:

# calculate correlation between time series
nino_corr = np.corrcoef(SST_nino_34_anom_ECCO_monthly_mean, SST_nino_34_anom_NOAA_monthly_mean)[1]
nino_ev   = 1 - np.var(SST_nino_34_anom_ECCO_monthly_mean-SST_nino_34_anom_NOAA_monthly_mean)/np.var(SST_nino_34_anom_NOAA_monthly_mean)
plt.figure(figsize=(8,5), dpi= 90)
plt.plot(SST_nino_34_anom_ECCO_monthly_mean.time, \
             SST_nino_34_anom_ECCO_monthly_mean - SST_nino_34_anom_ECCO_monthly_mean.mean(),'b.-')
plt.plot(SST_nino_34_anom_NOAA_monthly_mean.time, \
             SST_nino_34_anom_NOAA_monthly_mean - SST_nino_34_anom_NOAA_monthly_mean.mean(),'r.-')
plt.title('Niño 3.4 SST Anomaly \n correlation: %s \n explained variance: %s' % (np.round(nino_corr[0],3), \
                                                                                 np.round(nino_ev.values,3)));
plt.legend(('ECCO','NOAA'))
plt.ylabel('deg C');
plt.xlabel('year');
plt.grid()
_images/ECCO_v4_Example_calculations_with_scalar_quantities_73_0.png

ECCO is able to match the NOAA Niño 3.4 index fairly well.

Conclusion

You should now be familiar with doing some doing calculations using scalar quantities.

Suggested exercises

  1. Create the global mean SSH time series from 1992-2015

  2. Create the global mean sea level trend (map) from 1992-2015

  3. Create the global mean sea level trend (map) for two epochs 1992-2003, 2003-2015

  4. Compare other Niño indices