Merge pull request #257 from 4minakov/newtutorial

new geophysics tutorial draft and MT dataset added

vignettes/data/MT_Svalbard_Z.csv

@@ -0,0 +1,37 @@
+T,Code,Lat,Lon,X,Y,Zxx,DZxx
+64,W01,79.439116,13.371895,344352.8496,8829185.331,0.52156-1.2687i,0.054863
+128,W01,79.439116,13.371895,344352.8496,8829185.331,0.97673-0.94615i,0.079632
+256,W01,79.439116,13.371895,344352.8496,8829185.331,1.275-0.54561i,0.10554
+64,W02,79.471702,13.403403,345465.1801,8832708.491,0.36581+0.027603i,0.012205
+128,W02,79.471702,13.403403,345465.1801,8832708.491,0.3344+0.13025i,0.019636
+256,W02,79.471702,13.403403,345465.1801,8832708.491,0.33298+0.071971i,0.016023
+64,W03,79.35887,14.091742,357911.5152,8818456.864,0.67767+0.53303i,0.01449
+128,W03,79.35887,14.091742,357911.5152,8818456.864,0.6028+0.52425i,0.022896
+256,W03,79.35887,14.091742,357911.5152,8818456.864,0.46475+0.40616i,0.029995
+64,W04,79.392861,14.0056,356602.2036,8822435.623,-0.44539-1.2946i,0.012704
+128,W04,79.392861,14.0056,356602.2036,8822435.623,-0.12597-0.99152i,0.016897
+256,W04,79.392861,14.0056,356602.2036,8822435.623,0.13967-0.52094i,0.021547
+64,W05,79.372686,13.864863,353455.5873,8820549.936,0.8007+0.63546i,0.018563
+128,W05,79.372686,13.864863,353455.5873,8820549.936,0.90669+0.78722i,0.025465
+256,W05,79.372686,13.864863,353455.5873,8820549.936,0.54578+0.61109i,0.022635
+64,W06,79.337495,13.893527,353563.6793,8816578.44,0.85439+0.74466i,0.039878
+128,W06,79.337495,13.893527,353563.6793,8816578.44,1.0516+1.1238i,0.054053
+256,W06,79.337495,13.893527,353563.6793,8816578.44,0.49092+0.83855i,0.042233
+64,W07,79.42051667,13.736,351485.8998,8826175.357,0.81914+0.80551i,0.023629
+128,W07,79.42051667,13.736,351485.8998,8826175.357,0.78128+0.68671i,0.030426
+256,W07,79.42051667,13.736,351485.8998,8826175.357,0.51883+0.59256i,0.030384
+64,W08,79.431243,13.699478,350892.5935,8827456.802,0.39831-0.1402i,0.01706
+128,W08,79.431243,13.699478,350892.5935,8827456.802,0.44428-0.09936i,0.022626
+256,W08,79.431243,13.699478,350892.5935,8827456.802,0.42242+0.15965i,0.028102
+64,W09,79.49686111,13.91783333,356218.9673,8834177.87,-0.12304-0.52841i,0.00536
+128,W09,79.49686111,13.91783333,356218.9673,8834177.87,-0.011371-0.37717i,0.0055408
+256,W09,79.49686111,13.91783333,356218.9673,8834177.87,0.040909-0.16656i,0.0060876
+64,W10,79.472456,13.944975,356438.2811,8831406.374,0.19744-0.74573i,0.014357
+128,W10,79.472456,13.944975,356438.2811,8831406.374,0.4019-0.45596i,0.021548
+256,W10,79.472456,13.944975,356438.2811,8831406.374,0.50063+0.068497i,0.027314
+64,W11,79.486643,13.235664,342293.8848,8834811.661,0.82246+2.7602i,0.016108
+128,W11,79.486643,13.235664,342293.8848,8834811.661,0.38066+1.8882i,0.02087
+256,W11,79.486643,13.235664,342293.8848,8834811.661,0.17339+1.0451i,0.027693
+64,W12,79.425669,13.378223,344285.4562,8827679.973,0.69293-0.38905i,0.032038
+128,W12,79.425669,13.378223,344285.4562,8827679.973,0.90861-0.030763i,0.0516
+256,W12,79.425669,13.378223,344285.4562,8827679.973,0.813+0.39348i,0.073448

vignettes/tutorial_geophysics.Rmd

@@ -0,0 +1,149 @@
+---
+title: "4DModeller for geophysical signals"
+output: html_notebook
+---
+
+## Introduction
+This is a tutorial to apply R-INLA to modeling geophysical data. The data represent time-series at some locations distributed over a line or over an area. There are two possible case studies one is seismic data with seismic stations located along a line. Each station record 4 components of acoustic signal X-Y-Z particle displacement on geophone and presure component on hydrophone sensor. The second case is magnteotelluric data with stations distrubuted over some area. The signal has 4 channels: 2 magnetic and 2 electric field channels. In the dataset only 2-3 stations measure at the same. The spatiotemporal evolution of the field is governed by Maxwell equations. The source of EM field is disturbance of ionosphere due to solar activity, electrical structure of the crust and noise component (cultural EM noise, wind, rain, local conductors). The goal is to describe the source signal components in space in time.  
+The spatiotemporal evolution of the field is governed by wave equation. The small-scale heterogenetities in earth crust produce multiply scattered wavefield, getting more expressed at a later times and called coda wave. The goal is to learn about correlations between signals at diffent stations and from this predict distribution of heterogeneities.     
+
+## Pre-processing and Import data
+
+```{r}
+# Set the path to the CSV file
+data_path <- "data/MT_Svalbard_Z.csv"
+
+# Read the CSV file into a data frame
+d <- read.csv(data_path)
+
+# Display the first 30 rows of the data frame
+print(head(d, 30))
+
+# Create a scatter plot
+plot(d$X, d$Y, pch = 20, main = 'Svalbard MT sites', xlab = 'X', ylab = 'Y')
+
+```
+
+```{r}
+locations = d[, c('Lon', 'Lat')]
+locations = unique(locations)
+names(locations)=c('LONG','LAT')
+```
+
+## Meshing
+```{r}
+
+mesh <- fmesher::fm_mesh_2d(loc.domain = locations,
+                                max.edge = 0.05,
+                                cutoff = 1e-3,
+                                offset=0.1
+                                )
+plot(mesh)
+points(locations, col = "red")
+
+fdmr::plot_mesh(mesh)
+
+```
+
+## Stochastic modeling
+```{r}
+library(INLA)
+
+# Synthetic data generation
+set.seed(123) # For reproducibility
+n <- 100 # Number of time points
+stations <- 12 # Number of stations
+
+# Generate time index
+time_index <- 1:n
+
+# Generate spatial index (station IDs)
+space_index <- rep(1:stations, each = n)
+
+# Simulate some harmonic signals with noise for three stations
+harmonic_data <- data.frame(
+  station = factor(space_index),
+  time = rep(time_index, stations),
+  observation = sin(rep(time_index, stations) * 2 * pi / 50) + 
+                rnorm(n * stations, sd = 0.5) +
+                rep(rnorm(stations, sd = 3), each = n) # Station-specific offset
+)
+
+# Define the model with harmonic terms for time and spatial correlation
+formula <- observation ~ f(station, model = "iid") +
+  f(time, model = "rw1", cyclic = TRUE)
+
+# Fit the model using INLA
+result <- inla(formula, family = "gaussian", data = harmonic_data)
+
+# Display the summary of the results
+summary(result)
+
+# Visualize the fitted values
+plot(harmonic_data$time, harmonic_data$observation, col = harmonic_data$station, pch = 19, cex = 0.5, xlab = "Time", ylab = "Observation", main = "MT signals at 12 stations")
+points(harmonic_data$time, result$summary.fitted.values$mean, pch = 4, cex = 0.7, col = "blue")
+legend("topright", legend = c("Observations", "Fitted"), col = c("black", "blue"), pch = c(19, 4))
+
+
+# Extract the hyperparameters of the spatial field
+spatial_hyperparams <- result$summary.hyperpar
+
+# Print the hyperparameters
+print(spatial_hyperparams)
+
+```
+## Observed time series will come here.. 
+
+
+## This part doesn't work
+
+```{r}
+simdf <- data.frame(
+  ID = d$Code,
+  time = d$T,
+  datn = Re(d$Zxx),
+  E = Re(d$DZxx),
+  LONG = d$Lon,
+  LAT = d$Lat
+)
+
+
+initial_range = 0.1
+spde <- INLA::inla.spde2.pcmatern(
+  mesh = mesh,
+  prior.range = c(initial_range, 0.5),
+  prior.sigma = c(1, 0.01)
+)
+
+sigma0 <- 100
+range0 <- 0.03
+kappa0 <- sqrt(8 / 1) / range0
+tau0 <- 1 / (sqrt(4 * pi) * kappa0 * sigma0)
+
+inla.seed <- sample.int(n = 1E4, size = length(d$T))
+Q <- INLA::inla.spde.precision(spde, theta = c(log(tau0), log(kappa0)))
+x.mat <- matrix(NA, ncol = length(d$T), nrow = mesh$n)
+for (co in 1:ncol(x.mat)) {
+  x.mat[, co] <- INLA::inla.qsample(n = 1, Q)
+}
+
+A <- INLA::inla.spde.make.A(mesh = mesh, loc = locations)
+x.dat <- matrix(NA, ncol = n.time, nrow = n)
+for (t in 1:ncol(x.dat)) {
+  x.dat[, t] <- drop(A %*% x.mat[, t])
+}
+
+
+alpha <- 0.9
+sp.mat <- matrix(NA, ncol = n.time, nrow = n)
+sp.mat[, 1] <- x.dat[, 1]
+for (t in 2:n.time) {
+  sp.mat[, t] <- alpha * sp.mat[, t - 1] + x.dat[, t]
+}
+
+beta0 <- 0.5
+sigma_e <- 0.1
+lin.pred <- beta0 + sp.mat + matrix(rnorm(n * n.time, 0, sigma_e), ncol = n.time)
+
+```
+