# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: my_conda_env
#     language: python
#     name: python3
# ---

# # II- Demo. Optimal Interpolation

# An example of simulated SSH data access is provided in the "example_data_access_meom.ipynb" notebook. Here, an example of a mapping technique based on a simple optimal interpolation is proposed. The notebook is structured as follow: 
#
#     1) set optimal interpolation parameters,
#     2) reading of pseudo-observations,
#     3) perform optimal interpolation and,
#     4) save the results (reconstructed SSH field)
#
#
# Here, we assume a vector of observations, noted $y$ defined as:
#
# $$y = H x + \epsilon $$
#
# where $H$ is a linear observation operator between the reconstruction grid space and the observation space
# , $x$ is the state to estimate and $\epsilon$ is an independent observation error.
#
# The optimal interpolation consists in estimating an analysed state $x_{a}$ in combining the available observations to approximate the real state $x$:
#
# $$x_{a} = K y $$
# where $K$ is the weigth matrix defined as:
#
# $$ K = BH^T(HBH^T + R)^{-1} $$
#
# $B$ is the covariance matrix of $x$, and $R$ the covariance matrix of the error vector $\epsilon$ ($^T$ is the transpose operator)

import xarray as xr
import numpy
import warnings
import logging
import sys
import os
warnings.filterwarnings('ignore')

logger = logging.getLogger()
logger.setLevel(logging.INFO)

sys.path.append('..')

from src.mod_oi import *
from src.mod_inout import *

# ### 1) set optimal interpolation parameters

# +
# OI Grid
lon_min = 295.                                           # domain min longitude
lon_max = 305.                                           # domain max longitude
lat_min = 33.                                            # domain min latitude
lat_max = 43.                                            # domain max latitude
time_min = numpy.datetime64('2017-01-01')                # domain min time
time_max = numpy.datetime64('2017-12-31')                # domain max time
dx = 0.2                                                 # zonal grid spatial step (in degree)
dy = 0.2                                                 # meridional grid spatial step (in degree)
dt = numpy.timedelta64(1, 'D')                           # temporal grid step

glon = numpy.arange(lon_min, lon_max + dx, dx)           # output OI longitude grid
glat = numpy.arange(lat_min, lat_max + dy, dy)           # output OI latitude grid
gtime = numpy.arange(time_min, time_max + dt, dt)        # output OI time grid

# OI parameters
Lx = 1.                                                  # Zonal decorrelation scale (in degree)
Ly = 1.                                                  # Meridional decorrelation scale (in degree)
Lt = 7.                                                  # Temporal decorrelation scale (in days)
noise = 0.05                                             # Noise level (5%)
# -

# ### Open your AVISO+ session: fill the ```<AVISO_LOGIN>``` and ```<AVISO_PWD>``` items below

my_aviso_session = rq.Session()
my_aviso_session.auth = ("<AVISO_LOGIN>", "<AVISO_PWD>")
url_alongtrack = 'https://tds.aviso.altimetry.fr/thredds/dodsC/2021a-SSH-mapping-OSE-along-track-data'
url_map = 'https://tds.aviso.altimetry.fr/thredds/dodsC/2021a-SSH-mapping-OSE-grid-data'

# ### 2) reading of pseudo-observations + define output folder

inputs = [f'{url_alongtrack}/dt_gulfstream_alg_phy_l3_20161201-20180131_285-315_23-53.nc', 
          f'{url_alongtrack}/dt_gulfstream_j3_phy_l3_20161201-20180131_285-315_23-53.nc', 
          f'{url_alongtrack}/dt_gulfstream_s3a_phy_l3_20161201-20180131_285-315_23-53.nc',
          f'{url_alongtrack}/dt_gulfstream_h2g_phy_l3_20161201-20180131_285-315_23-53.nc',
          f'{url_alongtrack}/dt_gulfstream_j2g_phy_l3_20161201-20180131_285-315_23-53.nc',
          f'{url_alongtrack}/dt_gulfstream_j2n_phy_l3_20161201-20180131_285-315_23-53.nc'] 

# Define outputs
output_directory = '../results/'
if not os.path.exists(output_directory):
    os.mkdir(output_directory)  
output_oi = f'{output_directory}/OSE_ssh_mapping_BASELINE.nc'

# ### 3) perform optimal interpolation

# %%time
# set OI param & grid
ds_oi1_param = oi_param(Lx, Ly, Lt, noise)
ds_oi1_grid = oi_grid(glon, glat, gtime)
# Read input obs + discard a bit...
coarsening = {'time': 5}
#ds_oi1_obs = read_obs(inputs, ds_oi1_grid, ds_oi1_param, coarsening)
ds_oi1_obs = read_obs_from_aviso(inputs, my_aviso_session, ds_oi1_grid, ds_oi1_param, coarsening)
# Run OI (take 1h on my laptop)
for it in range(len(gtime)):
    oi_core(it, ds_oi1_grid, ds_oi1_param, ds_oi1_obs)

# ### 4) save the results (reconstructed SSH field)

url_ds_mdt = f'{url_map}/mdt.nc'
ds_oi1_grid = reformate_oi_output(ds_oi1_grid, url_ds_mdt, my_aviso_session)
ds_oi1_grid.to_netcdf(output_oi)



StarcoderdataJupyterScriptsDedupFiltered