Commit 66c7c9dc authored by Clément Haëck's avatar Clément Haëck
Browse files

Clean temperature separation computation

parent bedf6483
"""Compute separation Temperature between South and North zones.
- Take distribution of SST values above 32°N.
The distribution presents a peak at the higher values, corresponding to
the jet and south of it.
- Fit the peak using a gaussian and define the threshold as the center of the
gaussian minus twice its standard deviation.
"""
import warnings
import numpy as np
from scipy.optimize import curve_fit
import xarray as xr
import lib
import lib.data.jet_zones
import lib.data.ostia
import lib.zones
def main():
args = lib.get_args(['region', 'days', 'year', 'fixes'])
args['fixes']['Y'] = args['year']
sst = lib.data.ostia.get_data(args)
sst['zone'] = lib.zones.get_data(args, grid=lib.data.ostia.grid)['total']
sst = sst.sel(lat=slice(32, None))
sst['analysed_sst'] = sst.analysed_sst.where(sst.zone)
sep = xr.apply_ufunc(
get_separation, sst.analysed_sst,
input_core_dims=[['lon', 'lat']],
output_core_dims=[[]],
dask='parallelized',
output_dtypes=[float])
sep.name = 'threshold'
ofile = lib.data.jet_zones.get_filename(args['fixes'])
lib.check_output_dir(ofile, file=True)
sep.to_dataset().to_netcdf(ofile)
return sep
def get_separation(sst):
n_time = sst.shape[0]
sep = np.zeros(n_time)
for i in range(n_time):
sep[i] = _get_separation(sst[i])
return sep
def _get_separation(sst):
hist, bins, n_bins = get_hist(sst, width=0.5)
# Clean up first values of histogram
# Lowest values can sometimes spike. Because of ice maybe ?
hist[:3] = hist[3]
# Find peak
imax = np.argmax(hist)
hmax = hist[imax]
# Select 5° of SST around it
n_step = int(2.5/np.diff(bins).mean())
slc = slice(max(imax-n_step, 0),
min(imax+n_step, n_bins))
x = bins[slc]
h = hist[slc]
# Baseline
B = min(hmax/4, np.amin(h))
# Guess
p0 = [B, hmax-B, bins[imax], 1.]
bounds = [[-1e-3, hmax/3., x[0], 0.1],
[hmax/4, 2*hmax, x[-1], 5]]
def gauss(x, B, A, x0, std):
return B + A*np.exp(-0.5*(x-x0)**2/std**2)
try:
p, _ = curve_fit(gauss, x, h, p0, bounds=bounds)
except RuntimeError:
# If could not fit, the guess will do
p = p0
warnings.warn('Could not fit.')
x0, std = p[2:]
threshold = x0 - 2*std
return threshold
def get_hist(sst, width=0.5):
vmin = np.nanmin(sst)
vmax = np.nanmax(sst)
n_bins = len(np.arange(vmin, vmax, width))
hist, bins = np.histogram(sst, bins=n_bins, density=False,
range=(vmin, vmax))
hist = hist / np.sum(hist*np.diff(bins))
return hist, bins, n_bins
if __name__ == '__main__':
sep = main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment