alexander-barth / dataassim.jl

Implementation of various ensemble Kalman Filter data assimilation methods in Julia

DataAssim

The packages implements various data assimilation methods:

• (Extended) Kalman Filter
• Incremental 4D-Var
• Ensemble Square Root Filter (EnSRF)
• Ensemble Square Root Filter with serial processing of the observations (serialEnSRF)
• Ensemble Transform Kalman Filter (ETKF)
• Ensemble Transform Kalman Filter (EAKF)
• Singular Evolutive Interpolated Kalman filter (SEIK)
• Error-subspace Transform Kalman Filter (ESTKF)
• Ensemble Kalman Filter (EnKF)

The Julia code is ported from the Matlab/Octave code generated in the frame of the Sangoma project (http://data-assimilation.net/).

Example

An example of using to package is available as a jupyter-notebook.

Kalman Filter

The example below is the result of the Kalman Filter. The red elipse corresponds to the model forecast (Gaussian probability density function), the blue elipse corresponds to the observations and the purple elipse is the analysis (after assimilation).

illustrate_analysis_model.mp4

Ensemble Transform Kalman Filter

The example below is the result of the Ensemble Transform Kalman Filter. The red dots corresponds to the model ensemble, the blue elipse is the Gaussian probability density function of the observations and the purple dots is the model ensemble after analysis.

illustrate_analysis_model.mp4

Reference

Most of the algorithms are described in the review article:

Sanita Vetra-Carvalho, Peter Jan van Leeuwen, Lars Nerger, Alexander Barth, M. Umer Altaf, Pierre Brasseur, Paul Kirchgessner, and Jean-Marie Beckers. State-of-the-art stochastic data assimilation methods for high-dimensional non-Gaussian problems. Tellus A: Dynamic Meteorology and Oceanography, 70(1):1445364, 2018. doi: 10.1080/16000870.2018.1445364.

dataassim.jl's Issues

Support for non-uniform observations?

Hello!

Before I try to call an external model in C++, I wanted to know if there is a support for non-uniform observations? For example, here's our constraint with a hydrological model. We want to assimilate snow data that are gathered "randomly" by human.

• Our observations are not "periodic" : most years, but not always, observations are available on February 15th, March 1st-15th, April 1st (those date, for a given year, can be displaced by some days. i.e. instead of Feb 15, March1 and 15, it could be Feb 16, Feb 28, March 16, ...)
• Rest of the year, there are no observations
• Some year, we have continuous obs between February and May

Cheers|

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Small question: parallellization

Hi!

I have a small question, if I may.
Have you included any parallelization (over the particles) in this package?

Thanks a lot,
Ignace

Ensemble Kalman Smoother

Do you have knowledge of Ensemble Kalman Smoother? Is is possible to implement that in this framework?

Plans for DataAssim?

Hello!

I'm doing some prospective work on data assimilation. I'm gathering which packages or solutions is available (in any languages actually).

I was wondering two things with respect to DataAssim:

1. Is there a notebook or examples somewhere that demonstrate how to use the package? Perhaps I'm dense but I don't fully grasp the details.
2. Is there a plan to also introduce other types of assimilation such as particle filter? Or a simple Kalman filter (seems the method are only of the ensemble type).

Thanks for your time and any info!

Some help to define a model

Dear Alexander
I am trying to use your library with this model

x(k) = 0.95 x(k-1) + eta(k)
y(k) = x(k) + epsilon(k)

eta and epsilon are noises with std Q and R

I used
gaussian = ModelMatrix(0.95 .* Matrix(I,n,n))

but i did not find how to add the noise eta, the value Q is not used in the FreeRun function.
Thanks for your help and for sharing your work

Wrong type called in example.ipynb

Hi,

I found in this line of example.ipynb

xai, = fourDVar(xi, Pi, ℳ, yo, R, H, nmax, no);

this error

ERROR: LoadError: MethodError: no method matching fourDVar(::Vector{Int64}, ::Matrix{Bool}, ::ModelMatrix{Matrix{Float64}}, ::Matrix{Float64}, ::Matrix{Float64}, ::Matrix{Int64}, ::Int64, ::StepRange{Int64, Int64})

You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.

Closest candidates are:
  fourDVar(::AbstractVector, ::Any, ::Any, **::AbstractVector**, **::AbstractVector**, ::Any, ::Any, ::Any; innerloops, outerloops, tol)

in particular yo and R are defined as AbstractVector in fourDVar function defined in fourDVar.jl
while in example.ipynb they are defined as Matrix{Float64}

m=1;

R = 0.1 * Matrix(I, m, m)

no = 10:5:nmax
yo = randn(m, length(no))

Thank you in advance for your support!