Possible improvement for better readability of the modules of pyCATHY:

Currently, the run Data Assimilation function controlling the flow of the DA is part of the main CATHY class. For readability, this should be part of a different class i.e. class DA()

Remove the Soboldependency as it as been merged within Scipy.

See tweet

Pamphile Roy (@PamphileRoy) twitteó a las 4:45 p. m. on jue, ene 19, 2023:
🎉 we merged Sobol' sensitivity indices into @SciPy_team! Was a long effort and a great collaboration between maintainers and experts! Thanks @cziscience for funding it. This will facilitate the use of sensitivity analysis by the biomedical community!

https://t.co/O4JeqHiR8s
(https://twitter.com/PamphileRoy/status/1616099611952701445?t=uGDnJnEHcz3Hwjaj519vrQ&s=03)

Convert Saturation cathy <-> Resistivity resipy model

• reader for sw file
• reader resipy object/ file
• Archie law tranformation
• mesh interpolation

Improvement of spatial DA considering localisation

Nice explication from https://doi.org/10.1038/s41598-020-75710-5

"
To enhance EnKF performance during assimilation ensemble inflation and localization are applied. It has been shown by literature85,86 that ensemble-based data assimilation methods are sensitive to the size of ensemble. Generally, a larger number of ensemble members can better span the state-observation space and lead to better results but at the expense of strongly increased computation needs. To address this, ensemble inflation and localization methods are usually used to tackle filter divergent or inaccurate estimation87 for a small ensemble size and to avoid filter inbreeding. Ensemble inflation increases ensemble deviation from the ensemble-mean by applying a small coefficient ([1.1−1.3]
for the parameter and state updates) to ensemble members88. Localization using the Local Analysis (LA) scheme is also applied. It performs by spatially limiting the assimilation process within a certain distance from a grid point10,89. The suggested values (3∘
) by Khaki et al.10 are used as localization radii to achieve the best outcomes using a trial and error.
"

Examples of applications in prep:

• Noble dataset (Peruzzo et al.)
• Re della Pietra dataset
• Bousval dataset
• Synthetic dataset Weill

How?

from pyCATHY.DA.cathy_DA import DA
from pyCATHY.DA import perturbate, normalise, localisation
from pyCATHY.DA.cathy_DA import perturbate_parm, dictObs_2pd

# We create a DA object
simu = DA(dirName=path2prj, prj_name=prj_name_DA)

# We read all the observations (spatially distributed)
read_observations()
# Here is the new function to add: need to localize before creation of the data covariance matrice
localisation.localisation()
make_data_cov()

# Parameters pertubation
perturbate.perturbate()


# run DA simulation
simu.run_DA_sequential()

When downloading the src file from the CATHY gitbucket:

Soil file contains 20 layers while the dem parameters is set to nstr =15.

Problem description

The simulation consists in an infiltration in a confined domain (no flux on all the side boundaries). During DA update, after 3hours of infiltration, the state update generates a positive small value of pressure head (1e-3). CATHY breaks (the initial conditions being wrong) and the DA stops.

Simulation parameters

• Pressure head initial: -5
• Constant and uniform atmbc: 5e-7
• no input of ponding head initial conditions
• analysis method: Data Assimilation

In order to illustrate the concept of the wrapper using notebooks, it is necessary to write a notebook with some very general blocks showing:

• How to read function doc
• How to display CATHY variables
• ...

Up to know the grid3d file is generated while setting IPRT=3 during the run_processor step.
It would be useful to access grid3d without waiting for the process to end (espacially for time consuming simulations).

This should help to define the nodes of interest to generate the vp file by reading the grid3d file and searching for the nearest node from x,y,z coordinates.

We need a container tracking all operations during the Data Assimilation.

Using a dataframe including:

• variable of interest sampling/perturbated
• variable value before update/ after update
• time serie (ENKFT, assimilation times)

For clarity, mapping the model state to the observation should be done outside the CATHY DA class.
The map_states2observations function is currently included in the CATHY DA class.

map_states2observations defined 5 operators H:

• For ERT data, H is Archie law, SWC --> ER0 --> ERapp
• For SWC data, H is a multiplicator (porosity)
• For tensiometer data, no mapping is needed or VGP model, if the model state used, is saturation
• For discharge, no mapping is needed (direct comparison with modelled discharge)
• For scale data: no mapping needed (direct comparison with water balance)

Providing mapping has is own class one could use it its proper mapping relationship. CATHY DA class could inherit from the CATHY mapping class to parse it.

Issue with myst_nb

• When incorporating myst_nb into the Sphinx configuration file, the gallery generated from Sphinx containing examples with pyvista plots fails to display them. This issue is resolved by commenting out myst_nb.

• The second problem concerns notebook examples from literature written in myst format. Disabling the extension prevents them from being properly read. Currently, the notebooks are converted to .rst format using JupyterLab and are not executed with Sphinx.

Need a strategy to create temporary files (in particular .vtk) to save some CPU memory and space during Data Assimilation of ERT data.

ERT data assimilation generate:

• a vtk file for the converted mesh ER attribute
• a vtk file for the converted and interpolated mesh ER attribute

Those file are not used except for checking interpolation quality.

Same apply to the CATHY outputs. Would be interesting to store them in a tmp folder especially when running simulation on a Jupyter Notebook/ binder environement.

Data Assimilation of images

C. Lauvernet et al.
https://hal.inrae.fr/hal-02608128

• Usually, remote sensing data and sequences are under-used, though their content in information is very high (shapes evolution, correlations, . . . )
• HR Images would also help to identify the landscape elements (grass strips, hedges,. . . )
• In classical approaches: uncorrelated noise, diagonal error covariance matrices
• How to provide observation error covariance matrices adapted to spatially correlated errors? [2]
• Focusing on the observations operator description, and distances definition in the DA scheme

• https://www.frontiersin.org/articles/10.3389/frwa.2022.948832/full

[2] Chabot, V. et al., 2015. Accounting for observation errors in image data assimilation. Tellus A; Vol 67 (2015).

TO-DO

• Create a notebooks with some basics fct to use RS data, Observation Error Cov
• [ ]

Version pyCATHYv0.0.1 with CATHY 1.0

The simulation runs until TIME: 1.8026E+04, then stops with an "invalid memory reference.". I checked the cathyH file but i don't see apparent inconsistency with the simulation inputs.
Note that time 1.8026E+04, is the first value in the atmbc file with a negative net Precipitation/ET (i.e. -2.212e-08). Changing the value to 0 makes the simulation run until next time...

Attached (https://we.tl/t-sLFlTAlHI2) are the files to reproduce.

Improvement suggestion for DA

See Botto et al. 2018: "When assimilating multiple variables, proper normalization of the measurement error covariance matrices, anomalies of the simulated data, and innovation vectors were performed, using values of 0.6 m, 0.58, and 4.17 × 10−5 m3 s−1 for pressure head, water content and subsurface outflow, respectively. The normalization ensures that in multivariate assimilation scenarios the covariance matrices in the Kalman gain are not ill-conditioned (Evensen, 2003; Camporese et al., 2009b)."

Examples of applications in prep:

• Farming dataset (Carrera et al.)
• Noble dataset (Peruzzo et al.)
• Re della Pietra dataset
• Bousval dataset
• Synthetic dataset Weill

How?

from pyCATHY.DA.cathy_DA import DA
from pyCATHY.DA import perturbate, normalise
from pyCATHY.DA.cathy_DA import perturbate_parm, dictObs_2pd

# We create a DA object
simu = DA(dirName=path2prj, prj_name=prj_name_DA)

# We read all the observations (multiples types)
read_observations()
# Here is the new function to add: need to normalise before creation of the data covariance matrice
normalise.normalise()
make_data_cov()

# Parameters pertubation
perturbate.perturbate()

# run DA simulation
simu.run_DA_sequential()

Data assimilation should be handled outside Cathy's core. To reach v.0.0.2 milestone we need to implement the EnkF algorithm outside CATHY (in python).

The EnkF will rely on Sakov and Eversen (2009) mathematical formulation. See https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1600-0870.2009.00417.x?saml_referrer

Hi Ben,

Thank you for reducing all the dependencies to just blas and lapack, smooth installation.
On Fedora, I think they are blas-devel and lapack-devel rather than ...-dev, which is the right one on Ubuntu.

Do you think it is worth highlighting that pyCATHY already includes CATHY, and thus does not require a separate installation of CATHY before using pyCATHY as a wrapper?

Best,
Luca

Good practices in soft dev consist in implementing testing functions.

As for now, the wrapper is limited to exploratory testing.
To have a complete set of manual tests, here is a tentative list of the pyCATHY features:

• pycathy.tools
• test that the source files from the bitbucket are correctly fetched
• pycathy.DA.cathy_DA
• test that the DA returns an updated state
• pycathy.ERT.simulateERT
• test that the dataframe contains only positive values of apparent resistivities
• pycathy.plotters.cathy_plots
• test plots features (colorbar, export, ...)

These would be done automatically using pytest.

• Need to update the makefile with test run
• Need to update the requirements.txt with pytest

Milestone set to v.0.0.2