1 Introduction

In this document, we apply the pseudo polygon within mixture Gaussian model on GIA. The idea is that GIA is a stationary porcess on a subset of \(S^2\). The subset is given by removing polygons where the values are certainly zero according to experts from the entire sphere. In previous experiments, we have done a global stationary model, pseudo polygons with point observations at the polygon boundaries. In this document, we use a single pseudo polygon observation to update the process.

We generate dense mesh in the subset and sparse mesh insise the pseudo polygons. Pseudo observations are placed sparsely incide the polygons and along the bouondaries. The correlations length of the pseudo polygons are set to be larege. Why?

  1. process inside these polygons are uniformly zero. So either they are strongly correlated or they are independent but having tiny variance.

  2. we want a sparse mesh for these regions and to have a good approximation requires longer correlation length.

1.1 1 Generate the polygons and mesh

First we generate the mesh for a discreate reparesentation of the GIA process. The mesh is restricted on a subset of the sphere defined by the pseudo-polygons. More details about choosing and generating the polygons can be found here. The following chunk generate the pseudo polygons using the ensemble mean of 13 GIA mode solutions for a given threshold value, say \(0.3\).

## load library and functions
library(INLA)
library(sp); library(GEOmap); library(rgdal)
library(ggplot2); library(grid); library(gridExtra)
source("functions.R")
source("functions-barriers-dt-models-march2017.R")

Next we generate the mesh and separate the triangles in/out-side of the polygons.

#### Dense points over in the subset
fibo_points <- fiboSphere(N = 12960, L0 = TRUE)
pinPoly <- unlist(over(zeroPoly, SpatialPoints(coords = fibo_points), returnList=T))
fibo_inSub<- fibo_points[-pinPoly,]
plot(zeroPoly)
points(fibo_inSub, pch = ".")

#### Sparse points in the polygons
fibo_points <- fiboSphere(N = 500, L0=TRUE)
pinPoly <- unlist(over(zeroPoly, SpatialPoints(coords = fibo_points), returnList=T))
fibo_inPoly<- fibo_points[pinPoly,]
plot(zeroPoly)
points(fibo_inPoly, pch = ".")
points(fibo_inSub, pch = ".")

fibo_points_all <- rbind(fibo_inPoly, fibo_inSub)
mesh_points_xyz <- do.call(cbind, Lll2xyz(lat = fibo_points_all[,2], lon = fibo_points_all[,1]))
mesh <- inla.mesh.2d(loc = mesh_points_xyz, cutoff = 0.01, max.edge = 0.5)
summary(mesh) # give the desired number of vertices and triangles.
## 
## Manifold:    S2
## CRS: N/A
## Vertices:    12981
## Triangles:   25958
## Boundary segm.:  0
## Interior segm.:  0
## xlim:    -0.9973623 1
## ylim:    -0.9997098 0.9993412
## zlim:    -0.9999614 0.9999614

Now separate the triangles by the pseudo-polygons.

mesh <- dt.mesh.addon.posTri(mesh = mesh, globe = TRUE)
Tlonlat <- Lxyz2ll(list(x = mesh$posTri[,1], y = mesh$posTri[,2], z = mesh$posTri[,3]))
Tlonlat$lon <- ifelse(Tlonlat$lon >=0, Tlonlat$lon, Tlonlat$lon + 359)
mesh$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
TinPoly <- unlist(over(zeroPoly, SpatialPoints(coords=mesh$Trill), returnList=T))
TAll <- 1:mesh$t
ToutPoly <- TAll[-TinPoly]
Omega = dt.Omega(list(TinPoly, 1:mesh$t), mesh)
plot(mesh, t.sub = Omega[[2]])

plot(mesh, t.sub = Omega[[1]])

1.2 Set up prior

We set up the SPDE priors for the mixture model and simulate some data using the prior. The correlation length in the zero polygon regions is not sensitive and we set it to be much smaller than the modelling regions.

mu_r <- 500/6371
v_r <- (1000/6371)^2
mu_s <- 20
v_s <- 40^2

## Transform the parameters for the SPDE_GMRF approximation
Tlognorm <- function(mu, v){
  logv <- log(1 + v/mu^2)
  logmu <- log(mu^2) - 0.5*log(mu^2 + v)
  return(c(logmu, logv))
}
trho <- Tlognorm(mu_r, v_r)
tsigma <- Tlognorm(mu_s, v_s)

lsigma0 <- tsigma[1]
theta1_s <- tsigma[2]
lrho0 <- trho[1]
theta2_s <- trho[2]
lkappa0 <- log(8)/2 - lrho0
ltau0 <- 0.5*log(1/(4*pi)) - lsigma0 - lkappa0

1.3 2 Data preparation

1.4 2.1 GIA forward model solution

Load the GIA prior mean and GPS data and do the same prepreation as previous.

1.5 GPS data

Remove GPS data inside the pseudo-polygons.

GPS_inPoly <- unlist(over(zeroPoly, SpatialPoints(coords = cbind(GPSV4b$lon, GPSV4b$lat)), returnList=T))
GPS_All <- 1:nrow(GPSV4b)
GPS_outPoly <- GPS_All[-GPS_inPoly]
plot(GPSV4b[GPS_outPoly,c("lon", "lat")], pch = "+")

GPS_data <- GPSV4b[GPS_outPoly,]
GPS_loc <- do.call(cbind, Lll2xyz(lat = GPS_data$lat, lon = GPS_data$lon))
GPS_sp <- SpatialPoints(data.frame(lon = ifelse(GPS_data$lon>359.5, GPS_data$lon - 360, GPS_data$lon), 
                                   lat = GPS_data$lat), proj4string = CRS("+proj=longlat"))

GPS_idx <- over(GPS_sp, Plist)
GPS_mu <- ice6g$trend[GPS_idx]
GPS_data$trend0 <- GPS_data$trend - GPS_mu

We also add some pseudo observations along the boudaries of the polygons to make smooth transition of the predictions and the mesh nodes incide the polygons. These values are set to be the ice6 values at those locations with opposite signs.

## get the boundary of the polygons
boundlines <- as(zeroPoly, 'SpatialLines') 
obs_bounds <- spsample(boundlines, n = 50, type = "regular") # note points more than specified
## The mesh nodes incide the polygons
obs_inpoly <- SpatialPoints(coords = mesh$Trill[TinPoly,])
obs_pseudo <- rbind(obs_bounds, obs_inpoly)
proj4string(obs_pseudo) <- proj4string(Plist)

## Find the ice6g values
pobs_idx <- over(obs_pseudo, Plist)
GIA_pobs <- ice6g$trend[pobs_idx]

nobsb <-nrow(obs_pseudo@coords)
obs_df <- data.frame(ID = rep("pseudo", nobsb), lon = obs_pseudo@coords[,1], lat = obs_pseudo@coords[,2],
                     trend = rep(0, nobsb), std = rep(0.05,nobsb), trend0 = -GIA_pobs)
obs_xyz <- do.call(cbind, Lll2xyz(lat = obs_pseudo@coords[,2], lon = obs_pseudo@coords[,1]))

GPS_all <- rbind(GPS_data, obs_df)
GPS_all_loc <- rbind(GPS_loc, obs_xyz)

1.6 3 Inference on the subset

Mesh_GIA <- mesh
A_data <- inla.spde.make.A(mesh = Mesh_GIA, loc = GPS_all_loc)
A_pred <- inla.spde.make.A(mesh = Mesh_GIA, loc = rbind(GPS_all_loc, Mesh_GIA$loc))

## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=GPS_all$trend0), A = list(A_data),
                     effects = list(GIA = 1:Mesh_GIA$n), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(A_pred),
                      effects = list(GIA=1:Mesh_GIA$n), tag = "pred")
stGIA <- inla.stack(st.est, st.pred)

## Fix the GPS errors
hyper <- list(prec = list(fixed = TRUE, initial = 0))
prec_scale <- c(1/GPS_all$std^2, rep(1, nrow(A_pred)))


Q.mixture = dt.create.Q(mesh, Omega, fixed.ranges = c(0.01, NA))
# - We fix the barrier range to a different value than we 
#   used for simulations
# - - Why? It does not matter, as long as it is 'small' 
#     the models are very
#     similar
# - - This shows that you do not need to know the 
#     true 'barrier range'!

log.prior = dt.create.prior.log.exp(prior.param = c(1,1))
# - The prior parameters are the lambdas in the exponential 
#   priors for standard deviation and inverse-range

GIA_spde = dt.inla.model(Q = Q.mixture, log.prior=log.prior)

formula <- y ~ -1 + f(GIA, model=GIA_spde)
# - The spatial model component is different from before
# - The rest of the model setup is the same! 
#   (as in the stationary case)
# - - e.g. the inla(...) call below is the same, 
#     only this formula is different

Then we run the INLA model. Note that this will take more than 10min and require memory larger than 32GB. We ran this on a server with 56 cores and 256GB memory.

res_inla <- inla(formula, data=inla.stack.data(stGIA), family = "gaussian",
                   scale =prec_scale, control.family = list(hyper = hyper),
                   control.predictor=list(A = inla.stack.A(stGIA), compute = TRUE))
INLA_pred <- res_inla$summary.linear.predictor

2 Analyse results

Now assemble the inla inference and prediction results.

INLA_pred <- res_inla$summary.linear.predictor

## Extract and project predictions
pred_idx <- inla.stack.index(stGIA, tag = "pred")$data
GPS_idx <- pred_idx[1:nrow(GPS_all)]
GIA_idx <- pred_idx[-c(1:nrow(GPS_all))]

## GPS 
GPS_u <- INLA_pred$sd[GPS_idx] 
GPS_pred <- data.frame(lon = GPS_all$lon, lat = GPS_all$lat, u = GPS_u)

## GIA
GIA_diff <- INLA_pred$mean[GIA_idx] 
GIA_m <- GIA_diff + GIA_prior
GIA_u <- INLA_pred$sd[GIA_idx]
proj <- inla.mesh.projector(Mesh_GIA, projection = "longlat", dims = c(360,180), xlim = c(0,360), ylim = c(-90, 90))
GIA_grid <- expand.grid(proj$x, proj$y)
GIA_pred <- data.frame(lon = GIA_grid[,1], lat = GIA_grid[,2],
                       diff = as.vector(inla.mesh.project(proj, as.vector(GIA_diff))),
                       mean = as.vector(inla.mesh.project(proj, as.vector(GIA_m))),
                       u = as.vector(inla.mesh.project(proj, as.vector(GIA_u))))

ress <- list(res_inla = res_inla, spde = GIA_spde, st = stGIA, 
             mesh = Mesh_GIA, GPS_pred = GPS_pred, GIA_pred = GIA_pred)

2.1 Plot the posteriors of the hyper parameters

theta1 <- res_inla$marginals.hyperpar$`Theta1 for GIA`
theta2 <- res_inla$marginals.hyperpar$`Theta2 for GIA`

Vmar<- inla.tmarginal(exp, theta1)
Rmar <- inla.tmarginal(exp, theta2)

## Find the mode of rho and sigma^2
lrho_mode <- res_inla$summary.hyperpar$mode[2]
lrho_mean <- res_inla$summary.hyperpar$mean[2]
lrho_sd <- res_inla$summary.hyperpar$sd[2]
rho_mode <- exp(lrho_mean - lrho_sd^2)

lsigma_mode <- res_inla$summary.hyperpar$mode[1]
lsigma_mean <- res_inla$summary.hyperpar$mean[1]
lsigma_sd <- res_inla$summary.hyperpar$sd[1]
sigma_mode <- exp(lsigma_mean - lsigma_sd^2)
plot(Vmar, type = "l", main = bquote(bold({sigma^2}("mode") == .(round(sigma_mode, 4)))))

plot(Rmar, type = "l", main = bquote(bold(rho("mode") == .(round(rho_mode, 4)))))

## The estimated correlation length is about 568km
rho_mode*6371
## [1] 377.124

2.2 Plot the predictions

GPS_pred <- ress$GPS_pred
GIA_pred <- ress$GIA_pred

map_prior <- map_res(data = ice6g, xname = "x_center", yname = "y_center", fillvar = "trend", 
                      limits = c(-7, 22), title = "Prior GIA mean field")

## Plot the GIA predicted mean
map_GIA <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "mean", 
                   limits = c(-7, 22), title = "Predicted GIA")

## Plot the GIA difference map
map_diff <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "diff", 
                   limits = c(-8, 8), title = "GIA difference: Updated - Prior")

## Plot the GIA difference map
map_sd <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "u", 
                    colpal = colorRamps::matlab.like(12),  title = "Predicted uncertainties")

## Display 
print(map_prior)

print(map_GIA)

print(map_diff)

print(map_sd)