Untitled

library(geoR)

kappa <- 2       # Smoothness parameter for Matérn
tau2 <- 1        # Partial sill
phi <- 120         # Range
sigma2 <- 1      # Nugget

# Define the Matérn function
matern_function <- function(h, tau2, phi, kappa, sigma2) {
  # Matérn covariance function with vectorized handling
  cov <- ifelse(h == 0, 
                sigma2 + tau2, 
                sigma2 + tau2 * (2^(1 - kappa) / gamma(kappa)) * (h / phi)^kappa * besselK(h / phi, kappa))
  return(cov)
}

set.seed(123)

grid_coords <- expand.grid(x = seq(0, 120, length.out=15),
                           y = seq(0, 120, length.out = 15))

grid_coords <- cbind(runif(250), runif(250))

# Simulate Gaussian random field using the Matérn model
simulated_field <- grf(n=250, grid = grid_coords, cov.model="matern", 
                       cov.pars=c(sigma2, phi), kappa=kappa, nugget=tau2)

vario1<-variog(simulated_field,option="cloud")

vario1$max.dist

plot(vario1,xlab="Distance (h), km")

vario2 <- variog(simulated_field, uvec=seq(0, vario1$max.dist, l=20), option = 'bin', estimator.type = "modulus")


plot(vario2)

matern_semivariogram <- function(h, tau2, phi, sigma2, kappa) {
  # Bessel function part of the Matérn function
  bessel_part <- (1 - ((sqrt(2 * kappa) * (h / phi))^ (kappa)) * (besselK(sqrt(2 * kappa) * (h / phi), kappa)))
  # Semivariogram formula
  semivariance <- sigma2 * (1 / (gamma(kappa) * 2^(kappa - 1))) * bessel_part
  return(semivariance)
}

kappa <- 1.5     # Smoothness
tau2 <- 0.01     # Partial sill
phi <- 120       # Range
sigma2 <- 0.1  # Smoothness parameter

# Generate distances (h) for the semivariogram
h_vals <- seq(0, vario2$max.dist, length.out = 100)

# Calculate semivariance for each distance h
semivariance_vals <- matern_semivariogram(h_vals, tau2, phi, sigma2, kappa)

# Plot the Matérn semivariogram curve
curve(matern_semivariogram(x, tau2, phi, sigma2, kappa), 
      from = 0, to = vario2$max.dist, col = "blue", 
      ylab = "Semivariance", xlab = "Distance (h)", 
      ylim = c(0, 0.6), add = TRUE)

vfit_ols = variofit(vario2, ini.cov.pars=c(0.1, 10), nugget = 0.05, fix.nugget = FALSE, fix.kappa = F, cov.model = 'matern', weights = 'equal')
lines(vfit_ols, col='red', lwd=1.5)
summary(vfit_ols)
vfit_wls = variofit(vario2, ini.cov.pars=c(0.01, 10), nugget = 0.05, fix.nugget = FALSE,cov.model = 'matern', weights = 'cressie')
lines(vfit_ols, col='purple', lwd = 1.5)
mlfit_exp = likfit(simulated_field, ini.cov.pars = c(0.1, 120), nugget = 0.5, fix.nugget = FALSE, cov.model = 'matern', lik.method = 'ML', trend = '1st')
summary(mlfit_exp)
remlfit_matern = likfit(simulated_field, ini.cov.pars = c(0.1, 120), nugget = 0.05, fix.nugget = F, cov.model = 'matern', lik.method = 'REML')
lines(remlfit_matern)
lines(mlfit_exp)