Commit 1624029a authored by Clément Haëck's avatar Clément Haëck
Browse files

Add script for frt proba total (for all data)

parent c4e99d90
......@@ -23,7 +23,7 @@ def main():
args = lib.get_args(['region', 'days', 'year',
'scale', 'number', 'fixes'],
add_args=add_args)
add_args)
if args['period'] != 'monthly':
raise ValueError("Script only supports monthly probability")
......
"""Compute probability for a front to be present.
For all data available.
For both S/N zones.
"""
from os import path
import numpy as np
import xarray as xr
import lib
import lib.data.hi
import lib.data.ostia
import lib.data.SN_separation
import lib.data.front_probability
def main():
def add_args(parser):
parser.add_argument('-threshold_N', type=float, default=15.)
parser.add_argument('-threshold_S', type=float, default=6.)
args = lib.get_args(['region', 'days',
'scale', 'number', 'fixes'],
add_args)
sst = lib.data.ostia.get_data(args)
hi = lib.data.hi.get_data(args)
sep = lib.data.SN_separation.get_data(args)
sep = lib.data.SN_separation.smooth(sep, time_step=8)
sst, hi, sep = xr.align(
sst, hi, sep,
exclude=['lat', 'lon'], join='inner')
HI = lib.data.hi.apply_coef(hi, lib.data.hi.get_coef(args))
zone_S = sst.analysed_sst > sep
frt = (HI.where(zone_S) > args['threshold_S'])
frt = frt.where(zone_S, HI > args['threshold_N'])
frt = frt.where(np.isfinite(sst.analysed_sst))
tot = (frt * 1.).mean('time')
var_name = 'p_frt'
tot = tot.to_dataset(name=var_name)
args['period'] = 'total'
outdir = lib.data.front_probability.get_root(args)
lib.check_output_dir(outdir)
year = '2000'
tot.to_netcdf(path.join(outdir, year, year + '0101.nc'),
encoding={var_name: {'zlib': True}})
return tot
if __name__ == '__main__':
tot = 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