eseis (version 0.7.3)

fmi_inversion: Invert fluvial data set based on reference spectra catalogue

Description

The fluvial model inversion (FMI) routine uses a predefined look-up table with reference spectra to identify those spectra that fit the empirical data set best, and returns the corresponding target variables and fit errors.

Usage

fmi_inversion(reference, data, n_cores = 1)

Value

List object containing the inversion results.

Arguments

reference

List containing lists with precalculated model spectra.

data

eseis object or numeric matrix (spectra organised by columns), empiric spectra which are used to identify the best matching target parameters of the reference data set.

n_cores

Numeric value, number of CPU cores to use. Disabled by setting to 1. Default is 1.

Author

Michael Dietze

Details

Note that the frequencies of the empiric and modelled data sets must match.

Examples

## NOTE THAT THE EXAMPLE IS OF BAD QUALITY BECAUSE ONLY 10 REFERENCE 
## PARAMETER SETS AND SPECTRA ARE CALCULATED, DUE TO COMPUATATION TIME
## CONSTRAINTS. SET THE VALUE TO 1000 FOR MORE MEANINGFUL RESULTS.

## create 100 example reference parameter sets
ref_pars <- fmi_parameters(n = 10,
                           h_w = c(0.02, 1.20),
                           q_s = c(0.001, 8.000) / 2650,
                           d_s = 0.01,
                           s_s = 1.35,
                           r_s = 2650,
                           w_w = 6,
                           a_w = 0.0075,
                           f_min = 5,
                           f_max = 80,
                           r_0 = 6,
                           f_0 = 1,
                           q_0 = 10,
                           v_0 = 350,
                           p_0 = 0.55,
                           e_0 = 0.09,
                           n_0_a = 0.6,
                           n_0_b = 0.8,
                           res = 100)

## create corresponding reference spectra
ref_spectra <- fmi_spectra(parameters = ref_pars)

## define water level and bedload flux time series
h <- c(0.01, 1.00, 0.84, 0.60, 0.43, 0.32, 0.24, 0.18, 0.14, 0.11)
q <- c(0.05, 5.00, 4.18, 3.01, 2.16, 1.58, 1.18, 0.89, 0.69, 0.54) / 2650
hq <- as.list(as.data.frame(rbind(h, q)))

## calculate synthetic spectrogram
psd <- do.call(cbind, lapply(hq, function(hq) {

  psd_turbulence <- eseis::model_turbulence(h_w = hq[1],
                                            d_s = 0.01,
                                            s_s = 1.35,
                                            r_s = 2650,
                                            w_w = 6,
                                            a_w = 0.0075,
                                            f = c(10, 70),
                                            r_0 = 5.5,
                                            f_0 = 1,
                                            q_0 = 18,
                                            v_0 = 450,
                                            p_0 = 0.34,
                                            e_0 = 0.0,
                                            n_0 = c(0.5, 0.8),
                                            res = 100, 
                                            eseis = FALSE)$power

  psd_bedload <- eseis::model_bedload(h_w = hq[1],
                                      q_s = hq[2],
                                      d_s = 0.01,
                                      s_s = 1.35,
                                      r_s = 2650,
                                      w_w = 6,
                                      a_w = 0.0075,
                                      f = c(10, 70),
                                      r_0 = 5.5,
                                      f_0 = 1,
                                      q_0 = 18,
                                      v_0 = 450,
                                      x_0 = 0.34,
                                      e_0 = 0.0,
                                      n_0 = 0.5,
                                      res = 100,
                                      eseis = FALSE)$power

  ## combine spectra
  psd_sum <- psd_turbulence + psd_bedload

  ## return output
  return(10 * log10(psd_sum))
}))

graphics::image(t(psd))

## invert empiric data set
X <- fmi_inversion(reference = ref_spectra, 
                   data = psd)

## plot model results
plot(X$parameters$q_s * 2650, 
     type = "l")
plot(X$parameters$h_w, 
     type = "l")