# Copyright 2013-2021 The Salish Sea NEMO Project and
# The University of British Columbia
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# A module to interpolate Live Ocean results onto Salish Sea NEMO grid and
# save boundary forcing files.
import datetime
import logging
import os
import gsw
import netCDF4 as nc
import numpy as np
import xarray as xr
from salishsea_tools import LiveOcean_grid as grid
from scipy import interpolate
logger = logging.getLogger(__name__)
logger.addHandler(logging.NullHandler())
# ---------------------- Interpolation functions ------------------------
[docs]
def load_SalishSea_boundary_grid(imin, imax, rim, meshfilename):
"""Load the Salish Sea NEMO model boundary depth, latitudes and longitudes.
:arg imin int: first grid point of western boundary
:arg imax int: one past last grid point of western boundary
:arg rim int: rim width, note rim starts at second grid point
:arg fname str: name of boundary file
:returns: numpy arrays depth, lon, lat and a tuple shape
"""
with nc.Dataset(meshfilename) as meshfile:
lonBC = meshfile.variables['nav_lon'][imin:imax, 1:rim + 1]
latBC = meshfile.variables['nav_lat'][imin:imax, 1:rim + 1]
depBC = meshfile.variables['gdept_1d'][0]
shape = lonBC.shape
return depBC, lonBC, latBC, shape
[docs]
def load_LiveOcean(
date,
LO_dir='/results/forcing/LiveOcean/downloaded/',
LO_file='low_passed_UBC.nc'
):
"""Load a time series of Live Ocean results represented by a date,
location and filename
:arg date: date as yyyy-mm-dd
:type date: str
:arg LO_dir: directory of Live Ocean file
:type LO_dir: str
:arg LO_file: Live Ocean filename
:type LO_file: str
:returns: xarray dataset of Live Ocean results
"""
# Choose file and load
sdt = datetime.datetime.strptime(date, '%Y-%m-%d')
file = os.path.join(LO_dir, sdt.strftime('%Y%m%d'), LO_file)
T = grid.get_basic_info(file, only_T=True) # note: grid.py is from Parker
d = xr.open_dataset(file)
return d
[docs]
def interpolate_to_NEMO_depths(
dataset, depBC, var_names
):
""" Interpolate variables in var_names from a Live Ocean dataset to NEMO
depths. LiveOcean land points (including points lower than bathymetry) are
set to np.nan and then masked.
:arg dataset: Live Ocean dataset
:type dataset: xarray Dataset
:arg depBC: NEMO model depths
:type depBC: 1D numpy array
:arg var_names: list of Live Ocean variable names to be interpolated,
e.g ['salt', 'temp']
:type var_names: list of str
:returns: dictionary containing interpolated numpy arrays for each variable
"""
interps = {}
for var_name in var_names:
var_interp = np.zeros((depBC.shape[0], dataset[var_name][0, 0].shape[0],
dataset[var_name][0, 0].shape[1]))
for j in range(var_interp.shape[1]):
for i in range(var_interp.shape[2]):
LO_depths = dataset.z_rho.values[0, :, j, i]
var = dataset[var_name].values[0, :, j, i]
var_interp[:, j, i] = np.interp(
-depBC, LO_depths, var, left=np.nan
)
# NEMO depths are positive, LiveOcean are negative
interps[var_name] = np.ma.masked_invalid(var_interp)
return interps
[docs]
def remove_south_of_Tatoosh(interps, imask=6, jmask=17):
"""Removes points south of Tatoosh point because the water masses there
are different"
:arg interps: dictionary of 3D numpy arrays.
Key represents the variable name.
:type var_arrrays: dict
:arg imask: longitude points to be removed
:type imask: int
:arg jmask: latitude points to be removed
:type jmask: int
:returns: interps with values south-west of Tatoosh set to Nan.
"""
for var in interps.keys():
for i in range(imask):
for j in range(jmask):
interps[var][:, i, j] = np.nan
interps[var] = np.ma.masked_invalid(interps[var][:])
return interps
[docs]
def fill_box(interps, maxk=35):
"""Fill all NaN in the LiveOcean with nearest points (constant depth),
go as far down as maxk
:arg interps: dictionary of 3D numpy arrays.
Key represents the variable name.
:type interps: dictionary
:arg maxk: maximum Live Ocean depth with data
:type maxk: int
:returns interps with filled values
"""
for var in interps.keys():
for k in range(maxk):
array = np.ma.masked_invalid(interps[var][k])
xx, yy = np.meshgrid(
range(interps[var].shape[2]), range(interps[var].shape[1])
)
x1 = xx[~array.mask]
y1 = yy[~array.mask]
newarr = array[~array.mask]
interps[var][k] = interpolate.griddata((x1, y1),
newarr.ravel(), (xx, yy),
method='nearest')
return interps
[docs]
def convect(sigma, interps):
"""Convect interps based on density (sigma).
Ignores variations in cell depths and convects vertically
:arg interps: dictionary of 3D numpy arrays.
Key represents the variable name.
:type interps: dictionary
:arg sigma: sigma-t, density, 3D array
:type sigma: numpy array
:returns sigma, interps stabilized
"""
small = 0.01
var_names = interps.keys()
kmax, imax, jmax = sigma.shape
good = False
while not good:
good = True
for k in range(kmax - 1):
for i in range(imax):
for j in range(jmax):
if sigma[k, i, j] > sigma[k + 1, i, j]:
good = False
for var in var_names:
interps[var][k, i, j], interps[var][
k + 1, i, j
] = interps[var][k + 1, i, j], interps[var][k, i, j
]
sigma[k, i, j], sigma[k + 1, i, j] = sigma[
k + 1, i, j
], sigma[k, i, j]
return sigma, interps
[docs]
def stabilize(sigma, interps):
"""Add a little salt to stabilize marginally
stable cells
:arg interps: dictionary of 3D numpy arrays.
Key represents the variable name.
:type interps: dictionary
:arg sigma: sigma-t, density, 3D array
:type sigma: numpy array
:returns interps stabilized
"""
small = 0.01 # stabilize for delta sigma less than this
kl = 25 # stabilize for low delta sigma higher than this
add_salt = 0.01 # add this much salt
kmax, imax, jmax = sigma.shape
for k in range(kl - 1):
for i in range(imax):
for j in range(jmax):
if sigma[k+1, i, j] - sigma[k, i, j] < small:
interps['salt'][:k+1, i, j] += -add_salt / (k+1)
interps['salt'][k+1:, i, j] += add_salt / (kmax - k+1)
return interps
[docs]
def extend_to_depth(interps, maxk=35):
"""Fill all values below level with LiveOcean data with data from above
start at maxk
:arg interps: dictionary of 3D numpy arrays.
Key represents the variable name.
:type interps: dictionary
:arg maxk: maximum Live Ocean depth with data
:type maxk: int
:returns interps extended to depth
"""
for var in interps.keys():
interps[var][maxk:] = interps[var][maxk - 1]
return interps
[docs]
def interpolate_to_NEMO_lateral(interps, dataset, NEMOlon, NEMOlat, shape):
"""Interpolates arrays in interps laterally to NEMO grid.
Assumes these arrays have already been interpolated vertically.
Note that by this point interps should be a full array
:arg interps: dictionary of 4D numpy arrays.
Key represents the variable name.
:type interps: dictionary
:arg dataset: LiveOcean results. Used to look up lateral grid.
:type dataset: xarray Dataset
:arg NEMOlon: array of NEMO boundary longitudes
:type NEMOlon: 1D numpy array
:arg NEMOlat: array of NEMO boundary longitudes
:type NEMOlat: 1D numpy array
:arg shape: the lateral shape of NEMO boundary area.
:type shape: 2-tuple
:returns: a dictionary, like var_arrays, but with arrays replaced with
interpolated values
"""
# LiveOcean grid
lonsLO = dataset.lon_rho.values.ravel()
latsLO = dataset.lat_rho.values.ravel()
# interpolate each variable
interpl = {}
for var in interps.keys():
var_new = np.zeros((interps[var].shape[0], shape[0], shape[1]))
for k in range(var_new.shape[0]):
var_grid = interps[var][k, :, :].ravel()
var_new[k, ...] = interpolate.griddata(
(lonsLO, latsLO), var_grid,
(NEMOlon, NEMOlat), method='linear'
)
interpl[var] = var_new
return interpl
[docs]
def calculate_Si_from_NO3(NO3, SA, a=6.46, b=1.35, c=0, sigma=1, tsa=29):
"""Use a simple fit to calculate Si from NO3
:arg NO3: 3-D array of nitate values
:type NO3: array
:arg SA: 3-D array of absolute salinities
:type SA: array
:arg float a: constant in Si from NO3 fit units uM
:arg float b: linear term in Si form NO3 fit units none
:arg float c: magnitude for salinity additional impact units uM
:arg float sigma: 1/width of tanh for salinity impact units /(g/kg)
:arg float tsa: centre of salnity correction units g/kg
:returns: a 3-D array of silicon values
"""
Si = a + b * NO3 + c * np.tanh(sigma * (SA - tsa))
Si[Si < 0] = 0
return Si
[docs]
def correct_high_NO3(NO3, smax=100, nmax=120):
"""Correct LiveOcean nitrates that are higher than smax, so that
the largest nitrate is nmax. Defaults cause no correction.
:arg NO3: 3-D array of nitrate values
:type NO3: array
:arg smax: highest nitrate value corrected
:type smax: float
:arg nmax: maximum nitrate value allowed
:type nmax: float
:returns: a 3-D array of corrected nitrate values"""
#correction = np.array([(nitrate - smax) if nitrate > smax else 0 for
# nitrate in NO3])
correction = NO3 - smax
correction[NO3 < smax] = 0.
newnitrate = NO3 - correction * correction / (correction + nmax - smax)
return newnitrate
[docs]
def prepare_dataset(interpl, var_meta, LO_to_NEMO_var_map, depBC, time):
"""Prepare xarray dataset for Live Ocean file
:arg interpl: dictionary of 4D numpy arrays.
Key represents the variable name.
:type interpl: dictionary SSS
:arg var_meta: metadata for each variable in var_arrays.
Keys are NEMO variable names.
:type var_meta: a dictionary of dictionaries with key-value pairs of
metadata
:arg LO_to_NEMO_var_map: a dictionary mapping between LO variable names
(keys) and NEMO variable names (values)
:type LO_to_NEMO_var_map: a dictionary with string key-value pairs
:arg depBC: NEMO model depths
:type depBC: 1D numpy array
:arg time: time from Live Ocean dataset
:type time: float
returns dataset
"""
# Add some global attributes
ds_attrs = {
'acknowledgements':
'Live Ocean http://faculty.washington.edu/pmacc/LO/LiveOcean.html',
'creator_email':
'sallen@eoas.ubc.ca',
'creator_name':
'Salish Sea MEOPAR Project Contributors',
'creator_url':
'https://salishsea-meopar-docs.readthedocs.org/',
'institution':
'UBC EOAS',
'institution_fullname': (
'Earth, Ocean & Atmospheric Sciences,'
' University of British Columbia'
),
'summary': (
'Temperature, Salinity, Nitrate, Oxygen, DIC and TALK'
'from the Live Ocean model'
' interpolated in space onto the Salish Sea NEMO Model'
' western open boundary. Silicon from Nitrate.'
),
'source': (
'http://nbviewer.jupyter.org/urls/bitbucket.org/'
'salishsea/.../LiveOceanNew'
),
'history': (
'[{}] File creation.'
.format(datetime.datetime.today().strftime('%Y-%m-%d'))
)
}
da = {}
var_names = (var for var in interpl.keys() if var != 'NH4')
for var in var_names:
da[var] = xr.DataArray(
data=interpl[var],
name=LO_to_NEMO_var_map[var],
dims=('time_counter', 'deptht', 'yb', 'xbT'),
coords={
'time_counter': time,
'deptht': depBC,
'yb': [1],
'xbT': np.arange(interpl[var].shape[3])
},
attrs=var_meta[LO_to_NEMO_var_map[var]]
)
ds = xr.Dataset(
data_vars={
'vosaline': da['salt'],
'votemper': da['temp'],
'NO3': da['NO3'],
'Si': da['Si'],
'OXY': da['oxygen'],
'DIC': da['TIC'],
'TA': da['alkalinity']
},
coords={
'time_counter': time,
'deptht': depBC,
'yb': [1],
'xbT': np.arange(interpl['salt'].shape[3])
},
attrs=ds_attrs
)
return ds
[docs]
def write_out_file(ds, date, file_template, bc_dir):
"""Write out the Live Ocean Data File
:arg ds: xarray dataset containing data
:type ds: xarray dataset
:arg date: date for file
:type date: str
:arg file_template: filename template for the saved files;
it will be formatted with a datetime.
:type file_template: str
:arg bc_dir: directory for boundary condition file
:type bc_dir: str
"""
sdt = datetime.datetime.strptime(date, '%Y-%m-%d')
filename = file_template.format(sdt)
filepath = os.path.join(bc_dir, filename)
encoding = {var: {'zlib': True} for var in ds.data_vars}
encoding['time_counter'] = {'units': 'minutes since 1970-01-01 00:00'}
ds.to_netcdf(
path=filepath,
unlimited_dims=('time_counter'),
encoding=encoding,
)
logger.debug('Saved {}'.format(filename))
return filepath
# ------------------ Creation of files ------------------------------
[docs]
def create_LiveOcean_TS_BCs(
date,
file_template='LiveOcean_v201712_{:y%Ym%md%d}.nc',
meshfilename='/results/nowcast-sys/grid/mesh_mask201702.nc',
bc_dir='/results/forcing/LiveOcean/boundary_conditions/',
LO_dir='/results/forcing/LiveOcean/downloaded/',
LO_to_SSC_parameters = {'NO3': {'smax' : 100.,
'nmax' : 120.,},
'Si' : {'a' : 6.46,
'b' : 1.35,
'c' : 0.,
'sigma' : 1.,
'tsa' : 29}}
):
"""Create a Live Ocean boundary condition file for date
for use in the NEMO model.
:arg str date: date in format 'yyyy-mm-dd'
:arg str file_template: filename template for the saved files;
it will be formatted with a datetime.
:arg str bc_dir: the directory in which to save the results.
:arg str LO_dir: the directory in which Live Ocean results are stored.
:arg dict LO_to_SSC_parameters: a dictionary of parameters to convert
Live Ocean values to Salish Sea Cast
:returns: Boundary conditions files that were created.
:rtype: list
"""
# Create metadeta for temperature and salinity
# (Live Ocean variables, NEMO grid)
var_meta = {
'vosaline': {
'grid': 'SalishSea2',
'long_name': 'Practical Salinity',
'units': 'psu'
},
'votemper': {
'grid': 'SalishSea2',
'long_name': 'Potential Temperature',
'units': 'deg C'
},
'NO3': {
'grid': 'SalishSea2',
'long_name': 'Nitrate',
'units': 'muM'
},
'Si': {
'grid': 'SalishSea2',
'long_name': 'Dissolved Silicon',
'units': 'muM'
},
'OXY': {
'grid': 'SalishSea2',
'long_name': 'Oxygen',
'units': 'muM'
},
'DIC': {
'grid': 'SalishSea2',
'long_name': 'Dissolved Inorganic Carbon',
'units': 'muM'
},
'TA': {
'grid': 'SalishSea2',
'long_name': 'Total Alkalinity',
'units': 'muM'
},
}
# Mapping from LiveOcean TS names to NEMO TS names
LO_to_NEMO_var_map = {
'salt': 'vosaline',
'temp': 'votemper',
'NO3': 'NO3',
'NH4': 'NH4',
'Si': 'Si',
'oxygen': 'OXY',
'TIC': 'DIC',
'alkalinity': 'TA',
}
# Load BC information
depBC, lonBC, latBC, shape = load_SalishSea_boundary_grid(
imin=376 - 1, imax=470, rim=10, meshfilename=meshfilename
)
# Load the Live Ocean File
d = load_LiveOcean(date, LO_dir)
# Depth interpolation
interps = interpolate_to_NEMO_depths(d, depBC, var_names=(var for var in LO_to_NEMO_var_map if var != 'Si'))
# Change to TEOS-10
var_meta, interps['salt'], interps['temp'] = _convert_TS_to_TEOS10(
var_meta, interps['salt'], interps['temp']
)
# Remove South of Tatoosh
interps = remove_south_of_Tatoosh(interps)
# Fill whole LiveOcean data box horizontally
interps = fill_box(interps)
# Calculate the density (sigma) and convect
sigma = gsw.sigma0(interps['salt'][:], interps['temp'][:])
sigma, interps = convect(sigma, interps)
# Fill Live Ocean Vertically
interps = extend_to_depth(interps)
# Interpolate Laterally onto NEMO grid
interpl = interpolate_to_NEMO_lateral(interps, d, lonBC, latBC, shape)
# Convect Again
sigmal = gsw.sigma0(interpl['salt'][:], interpl['temp'][:])
sigmal, interpl = convect(sigmal, interpl)
interpl = stabilize(sigmal, interpl)
# Rework indexes for NEMO
for var in interpl.keys():
interpl[var] = np.swapaxes(interpl[var], 1, 2)
interpl[var] = interpl[var].reshape(
1, interpl[var].shape[0], 1,
interpl[var].shape[2] * interpl[var].shape[1]
)
# Due to change in LiveOcean (May 22, 2023) add NH4 to NO3 to
# preserve previous behaviour (note LiveOcean NH4+NO3 evaluates
# better than NO3 against NO3 obs)
interpl['NO3'] = interpl['NO3'] + interpl['NH4']
# Calculate Si from NO3 using LiveOcean nitrate
interpl['Si'] = calculate_Si_from_NO3(
interpl['NO3'], interpl['salt'],
a=LO_to_SSC_parameters['Si']['a'],
b=LO_to_SSC_parameters['Si']['b'],
c=LO_to_SSC_parameters['Si']['c'],
sigma=LO_to_SSC_parameters['Si']['sigma'],
tsa=LO_to_SSC_parameters['Si']['tsa']
)
# Correct NO3 values
interpl['NO3'] = correct_high_NO3(
interpl['NO3'],
smax=LO_to_SSC_parameters['NO3']['smax'],
nmax=LO_to_SSC_parameters['NO3']['nmax']
)
# Prepare dataset
ts = d.ocean_time.data
ds = prepare_dataset(interpl, var_meta, LO_to_NEMO_var_map, depBC, ts)
# Write out file
filepath = write_out_file(ds, date, file_template, bc_dir)
return filepath
def _convert_TS_to_TEOS10(var_meta, sal, temp):
"""Convert Practical Salinity and potential temperature to Reference
Salinity and Conservative Temperature using gsw functions.
:arg var_meta: dictionary of metadata for salinity and temperature.
Must have keys vosaline and votemper, each with a sub
dictionary with keys long_name and units
:type var_meta: dictionary of dictionaries
:arg sal: salinity data
:type sal: numpy array
:arg temp: temperature daya
:type temp: numpy array
:returns: updated meta data, salinity and temperature"""
# modify metadata
new_meta = var_meta.copy()
new_meta['vosaline']['long_name'] = 'Reference Salinity'
new_meta['vosaline']['units'] = 'g/kg'
new_meta['votemper']['long_name'] = 'Conservative Temperature'
# Convert salinity from practical to reference salinity
sal_ref = gsw.SR_from_SP(sal[:])
# Convert temperature from potential to conservative
temp_cons = gsw.CT_from_pt(sal_ref[:], temp[:])
return new_meta, sal_ref, temp_cons