"""
Spectral Nudging preprocessings
"""
import logging
from ...log import LoggingLevelContext
with LoggingLevelContext(logging.WARNING):
import xarray as xr
from xoa.coords import get_time
from .io import read_blkdat_input, hycom2nc, read_regional_grid
from . import HYCOM3D_GRID_AFILE, HYCOM3D_BLKDAT_INPUT_FILE
import pandas as pd
from sloop.times import convert_to_julian_day
from sloop.interp import interp_time
import numpy as np
LOGGER = logging.getLogger(__name__)
[docs]def get_demerliac_date(ncin):
ds = xr.open_dataset(ncin)
return get_time(ds)
[docs]def demerliac_sel(ncins, outfile_pattern="{field}_timesel.nc", dt=None):
LOGGER.info("Time selection for demerliac filtering")
LOGGER.info("time delta (dt) in day is mandatory")
if dt==None:
dt = pd.to_timedelta(1.46, 'D')
else:
dt = pd.to_timedelta(dt, 'D')
LOGGER.debug("time delta: {dt}")
ncouts = {}
for ncin in ncins:
LOGGER.info(f"Open {ncin}")
ds = xr.open_dataset(ncin)
ds = ds.sel(time=slice(ds.time.mean()-dt, ds.time.mean()+dt, 1))
field = ncin.split('.')[0]
ncout = outfile_pattern.format(**locals())
ncouts[field] = ncout
LOGGER.info(f"Store in {ncout}")
ds.to_netcdf(ncout)
date_demerliac = convert_to_julian_day(ds.time.mean().values)
return date_demerliac, ncouts
[docs]def spectral_sel(ncins, time, outfile_pattern="{field}-mercator.nc", rank="0"):
LOGGER.info("Time selection for spectral filtering on nesting input files")
LOGGER.info("Read blkdat.input")
dsb = read_blkdat_input(HYCOM3D_BLKDAT_INPUT_FILE.format(rank=rank))
nx = int(dsb.idm)
ny = int(dsb.jdm)
nz = int(dsb.kdm)
ncouts = {}
for ncin in ncins:
if ncin.endswith('.b'):
LOGGER.info(f"Processing on {ncin}")
nc = hycom2nc(nx, ny, nz,
HYCOM3D_GRID_AFILE.format(rank=rank),
ncin)
ds = xr.open_dataset(nc)
ds = interp_time(ds, [str(time.values[0])], preserve_encoding=True)
field = (ncin.split('.')[0]).split('-')[0]
ncout = outfile_pattern.format(**locals())
ncouts[field] = {"mercator": ncout}
ds = ds.fillna(0)
LOGGER.info(f"Store results in {ncout}")
ds.to_netcdf(ncout)
return ncouts
[docs]def ncdiff(varname, nc1, nc2, outfile_pattern='{varname}-spn.nc', rank='0'):
LOGGER.info("Compute difference between filtered mercator and filtered hycom fields")
ds = read_regional_grid(HYCOM3D_GRID_AFILE.format(rank=rank), grid_loc='p')
ds = ds.rename({"idm":"x", "jdm":"y"})
ds1, ds2 = xr.open_dataset(nc1), xr.open_dataset(nc2)
if varname=="t":
longname = "temp"
elif varname=='s':
longname = "saln"
else:
longname = varname
field = ds2[longname].load().values-ds1[longname].load().values
field[field==0] = np.nan
da = xr.DataArray(field, name=varname,
dims=('time', 'lev', 'y', 'x')
)
da.encoding["_FillValue"]= np.nan
ds = xr.merge([ds, da])
ds = ds.assign_coords(time=ds2.time)
ncout = outfile_pattern.format(**locals())
LOGGER.info(f"Store results in {ncout}")
ds.to_netcdf(ncout)
return ncout