This document contains the code for all the examples and application to GIA in the paper. First we load the packages and functions.

library(sp)
library(INLA)
library(rgeos)
library(fields)
library(rgdal) 

source("functions.R")
source("partition_fun.R")

1 Triangulation on \(S^2\)

For the SPDE approximation, we need to generate triangulation on a sphere. The triangulation algorithm requires some starting points. We first generate the starting points for a evenly spaced triangulation for a stationary process.

1.1 Fibonacci points

The Fibonacci points are evenly distributed on the spehre; hence we use them as the starting points.

## Generates Fibonacci points on the sphere in xyz coords
fibo_xyz <- fiboSphere(N = 10000, LL = FALSE) 
## Project xyz to lon-lat coords
fibo_ll <- GEOmap::Lxyz2ll(list(x = fibo_xyz[,1], y = fibo_xyz[,2], z = fibo_xyz[,3]))
## Plot the points 
plot(fibo_ll$lat ~ fibo_ll$lon, pch = ".", xlab = "longitude", ylab = "latitude")

1.2 Stationary process

Now we generate the Triangulation for a stationary process using the Fibbonacci points generated above.

mesh1 <- inla.mesh.2d(loc = fibo_xyz, cutoff = 0.01, max.edge = 0.5)
summary(mesh1)

## 
## Manifold:    S2
## CRS: N/A
## Vertices:    10000
## Triangles:   19996
## Boundary segm.:  0
## Interior segm.:  0
## xlim:    -0.9996625 0.9996625
## ylim:    -0.9996324 0.9996558
## zlim:    -0.9999 0.9999

plot(mesh1, rgl = TRUE) # plot the mesh in 3D -- the rgl object can not be shown in the html

1.3 Subset process

Now we build the triangulation for a subset process that is defined on the ocean only.

## load the land mask and remove unecessary details
lands <- readOGR(dsn = "ne_110m_land", layer = "ne_110m_land")

## OGR data source with driver: ESRI Shapefile 
## Source: "ne_110m_land", layer: "ne_110m_land"
## with 127 features
## It has 2 fields

landsareas <- sapply(lands@polygons, function(x) slot(x, "area"))
lands2 <- lands[landsareas > 10,]

## Genreate the sub mesh
#-- use the stationary mesh and remove the triangles on the land
mesh2 <- dt.mesh.addon.posTri(mesh = mesh1, globe = TRUE) # add the centroid xyz coords for all the triangles
Tlonlat <- GEOmap::Lxyz2ll(list(x = mesh2$posTri[,1], y = mesh2$posTri[,2], z = mesh2$posTri[,3]))
mesh2$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
# -- Separate triangles in land and ocean
landid <- unlist(over(lands2, SpatialPoints(coords=mesh2$Trill, proj4string = CRS(proj4string(lands2))), returnList=T))
Allid <- 1:mesh2$t
Oceanid <- Allid[-landid]
Omega <- dt.Omega(list(landid, 1:mesh2$t), mesh2)

## Keep triangles in ocean only
mesh_ocean <- mesh.sub(mesh2, Omega, 2)
plot(mesh_ocean, rgl = T)

1.4 Partiton process

The following construct a triangulation for a partition process where the resolution of land is much lower than on the ocean.

## Generate sparse points on land
fibo_ll <- fiboSphere(N = 500, L0=FALSE) # Points in lon-lat coords [-180, 180]
landid <- unlist(over(lands2, SpatialPoints(coords = fibo_ll,  proj4string = CRS(proj4string(lands2))), returnList=T))
fibo_inland<- fibo_ll[landid,]
fibo_inland2 <- do.call(cbind, GEOmap::Lll2xyz(lat=fibo_inland[,2], lon = fibo_inland[,1]))

## Combine the sparse points in land and dense points in ocean to build the mesh
mesh3 <- inla.mesh.2d(loc = rbind(mesh_ocean$loc, fibo_inland2), cutoff = 0.01, max.edge = 0.5)
mesh3 <- dt.mesh.addon.posTri(mesh = mesh3, globe = TRUE) # Add on mesh$posTri: contains the positions of the triangles
Tlonlat <- GEOmap::Lxyz2ll(list(x = mesh3$posTri[,1], y = mesh3$posTri[,2], z = mesh3$posTri[,3]))
mesh3$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)

## Separate triangles in land and ocean
TinPoly <- unlist(over(lands2, SpatialPoints(coords=mesh3$Trill, proj4string = CRS(proj4string(lands2))), returnList=T))
TAll <- 1:mesh3$t
ToutPoly <- TAll[-TinPoly]
Omega <- dt.Omega(list(TinPoly, 1:mesh3$t), mesh3)
plot(mesh3, t.sub = Omega[[2]], col = "lightblue", rgl = TRUE )
plot(mesh3, t.sub = Omega[[1]], col = "yellow",  rgl = TRUE, add = TRUE) # Figure 2 shows a snapshot of the 3d plot

2 Modelling Non-stationarity

We proposed the subset and partion models to handle the non-stationarity. In the following we show two toy examples from simulations.

2.1 Subset model

The following codes generate the example in Section 3.1.

set.seed(18)
## Define the regions
loc.bnd <- matrix(c(0,0, 5,0, 5,5, 0,5), 4, 2, byrow=TRUE)
segm.bnd <- inla.mesh.segment(loc.bnd)
out.Poly <- SpatialPolygons(list(Polygons(list(Polygon(loc.bnd)), ID = "out")))
loc.int <- matrix( c(0.87,0.27, 0.13,0.13, 0.26,0.24, 0.07,0.97, 0.18,0.92, 0.31,0.32), byrow=T, 6, 2)*3 + 1
int.bound <- inla.nonconvex.hull(loc.int, convex=0.6)
int.Poly <- SpatialPolygons(list(Polygons(list(Polygon(int.bound$loc)), ID = "in")))
locPoly <- spsample(as(int.Poly, "SpatialLines"), n = 100, type = "regular")
segm.int <- inla.mesh.segment(locPoly)
## Sample dense and evenly spaced points in the entrie region
loc0 <- spsample(out.Poly, n = 500, type = "hexagonal", offset = c(0, 0))

## Triangulation for the entire regions
mesh4a <- inla.mesh.2d(loc = loc0,   boundary = segm.bnd, offset = c(0.5,0.5), cutoff = 0.2, max.edge = 0.5)
mesh4b <- inla.mesh.2d(loc = loc0,   boundary = segm.bnd, interior = list(segm.int), offset = c(0.5,0.5),
                      cutoff = 0.2, max.edge = 0.5)
## separate the triangles
mesh4b <- dt.mesh.addon.posTri(mesh4b)
TinP <- unlist(over(int.Poly, SpatialPoints(mesh4b$posTri), returnList=T))
Omega <- dt.Omega(list(TinP, 1:mesh4b$t), mesh4b)
Omega.SP <- dt.polygon.omega(mesh4b, Omega)

## Geneate subset mesh
mesh_sub <- mesh.sub(mesh4b, Omega, 2, sphere = FALSE)

## Create the precision matrix
## - The parameters
lsigma0 <- log(2)
lrho0 <- log(1)
lkappa0 <- log(8)/2 - lrho0
ltau0 <- 0.5*log(1/(4*pi)) - lsigma0 - lkappa0
## - The stationary SPDE model
spde1 <- inla.spde2.matern(mesh_sub, B.tau = matrix(c(ltau0, -1, 1),1,3), B.kappa = matrix(c(lkappa0, 0, -1), 1,3),
                  theta.prior.mean = c(0,0), theta.prior.prec = c(1,1))
Q1 <- inla.spde.precision(spde1, c(0,0))
## - The subset SPDE model
spde2 <- inla.spde2.matern(mesh4a, B.tau = matrix(c(ltau0, -1, 1),1,3), B.kappa = matrix(c(lkappa0, 0, -1), 1,3),
                           theta.prior.mean = c(0,0), theta.prior.prec = c(1,1))
Q2 <- inla.spde.precision(spde2, c(0,0))

## Calculate correlation at given locations
data_loc1 <- matrix(c(2.2, 2.8, 1.3,2.8, 3.1,2.8), byrow = TRUE, ncol = 2)
Amat1 <- inla.spde.make.A(mesh = mesh_sub, loc = data_loc1)
Amat2 <- inla.spde.make.A(mesh = mesh4a, loc = data_loc1)
Sigma_data1 <- cov2cor(Amat1%*%solve(Q1, t(Amat1)))
Sigma_data2 <- cov2cor(Amat2%*%solve(Q2, t(Amat2)))

plot(mesh_sub, asp = 1, main = "")
plot(Omega.SP[[2]], add=T, col='grey')
plot(mesh_sub, asp = 1, main = "", add = T)
points(data_loc1, col = 2, pch = 19)
text(data_loc1, pos = 3, c("A", "B", "C"), col = 2, cex = 2.2)

2.2 Partition model

The following codeds are for the example in Section 3.2.2.

## Sample sparse and evely spaced points in the inner region
loc1 <- spsample(int.Poly, n = 10, type = "regular", offset = c(0, 0))
subidx <- unlist(over(out.Poly, SpatialPoints(mesh_sub$loc[,1:2]), returnList=T))

## Combine the two sample points to generat the triangulation
mesh5 <- inla.mesh.2d(loc = rbind(mesh_sub$loc[subidx,1:2],loc1@coords),  boundary = segm.bnd, interior = list(segm.int), offset = c(0.5, 0.5), cutoff = 0.2, max.edge = 0.5 )

## separate the mesh
mesh5 <- dt.mesh.addon.posTri(mesh5)
TinP <- unlist(over(int.Poly, SpatialPoints(mesh5$posTri), returnList=T))
Omega2 = dt.Omega(list(TinP, 1:mesh5$t), mesh5)
Omega.SP2 = dt.polygon.omega(mesh5, Omega2)

## build the SPDE for the partition model and precision matrix
Q.function = dt.create.Q(mesh5, Omega2)
ranges = c(log(1.5), log(1))
Q3 = Q.function(theta = c(log(5), ranges))

## Calculate correlation at given locations
data_loc2 <- matrix(c(2.2,1.9, 2.2,3, 3.3, 3, 3.3, 1.9 ), byrow = TRUE, ncol = 2)
Amat3 <- inla.spde.make.A(mesh = mesh5, loc = data_loc2)
Sigma_data3 <- cov2cor(Amat3%*%solve(Q3, t(Amat3)))


plot(mesh5, main ="Mesh and Omega", asp = 1)
plot(Omega.SP2[[2]], add=T, col='grey')
plot(mesh5, add=T)
points(data_loc2, col = 2, pch = 19)
text(data_loc2, pos = 3, c("A", "B", "C", "D"), col = 2, cex = 2.2)

3 Application to GIA

The following shows the codes for the applciation to GIA in Section 4. Suppose we have loaded the GIA simulations and the GPS data.

3.1 Build the “Zero Region”

We use 8 existing GIA simulations to determine the “zero region”. Below shows how the zero region is generated.

## create the ensemble data sets
GIA_ensemble <- GIA_sim[[1]][, 2:5]
GIA_ensemble$trend <- rowMeans(sapply(GIA_sim, "[[", "trend"))
GIA_ensemble$std <- sqrt((rowSums(sapply(GIA_sim, function(x) (x$std)^2)))/8)

## Choose pixels smaller than 0.3mm/year as the zero region
GIA_sp <- GIA_ensemble
GIA_sp$zero1 <- abs(GIA_sp$trend) < 0.3
coordinates(GIA_sp) <- c("x_center", "y_center")
gridded(GIA_sp) <- TRUE

## Smooth the polygon by setting standard error smaller than 0.4mm/year
GIA_sp$zero2 <- abs(GIA_sp$trend) < 0.3 & abs(GIA_sp$std) < 0.4

## Merge the pixels into polygons
GIA_ras <- raster::raster(GIA_sp, layer = 4) # 4 is the column no of zero2
zeroPoly <- raster::rasterToPolygons(GIA_ras, fun=function(zero2) {zero2 == TRUE}, dissolve=TRUE )

## Remove polygons that are too small
zeroPolys <- zeroPoly@polygons[[1]]@Polygons
polyareas <- sapply(zeroPolys, function(x) x@area)
polyholes <- sapply(zeroPolys, function(x) x@hole)
zeropolys2 <- zeroPolys[polyareas > 200 ] 

zeroPolygon <- zeroPoly
zeroPolygon@polygons[[1]]@Polygons <- zeropolys2
proj4string(zeroPolygon) <- CRS("+proj=longlat")

map2 <- map("world2", interior = FALSE, plot = FALSE)
par(xaxs= "i", yaxs = "r")
plot(map2, type = "l", asp = 1, xlab = "longitude", ylab = "latitude", axes=FALSE)
axis(side = 1, at = c(0, 120, 240, 360))
axis(side = 2, at = c(-90, -45,0, 45, 90))
box()
plot(zeroPolygon, add = T, col = "blue") # Figure 5

3.1.1 Remove GPS stations in the zero region

The GPS observations in the zero region are removed since they do not contribute true process signals.

GPSid <- unlist(over(zeroPolygon, GPSsp, returnList=T))
GPS_data <- GPS[-GPSid,]
GPSsp <- GPSsp[-GPSid,]
GPS_locxyz <- GPS_locxyz[-GPSid,]
## Find the GIA prior at the GPS location and detrend
giaGPS <- over(GPSsp, ice6gsp)$trend
GPS_data$y <- GPS_data$trend - giaGPS

3.2 Set up priors for hyperparameters

Then We set up the priors for the hyperparameters which will be used in all models.

GPSdatasp <- GPS_data
coordinates(GPSdatasp) <- c("lon", "lat")
proj4string(GPSdatasp) <- CRS("+proj=longlat")
v_gia <- gstat::variogram(y ~ 1, data = GPSdatasp, cutoff = 1000)
plot(v_gia)

## Priors mean and variance for the parameters: rho and sigma
mu_r <- 1000/6371
v_r <- (2000/6371)^2
mu_s <- 1.5
v_s <- 3^2

## Transform the parameters for the SPDE_GMRF approximation
trho <- Tlognorm(mu_r, v_r)
tsigma <- Tlognorm(mu_s, v_s)

## set up for the INLA SPDE function
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

3.3 Stationary model

The following codes estimate and predict the stationary model.

## Triangulation for stationary process on sphere
fibo_points <- fiboSphere(N = 3e4, LL = FALSE)
mesh_st <- inla.mesh.2d(loc = fibo_points, cutoff = 0.01, max.edge = 0.5)
summary(mesh_st)

## 
## Manifold:    S2
## CRS: N/A
## Vertices:    30000
## Triangles:   59996
## Boundary segm.:  0
## Interior segm.:  0
## xlim:    -0.999934 0.999934
## ylim:    -0.9999457 0.9998971
## zlim:    -0.9999667 0.9999667

## Build the SPDE model
spde_st <- inla.spde2.matern(mesh = mesh_st, B.tau = matrix(c(ltau0, -1, 1),1,3), B.kappa = matrix(c(lkappa0, 0, -1), 1,3),
                              theta.prior.mean = c(0,0), theta.prior.prec = c(sqrt(1/theta1_s), sqrt(1/theta2_s)))

## Link the GIA process and GPS data
Adata_st <- inla.spde.make.A(mesh = mesh_st, loc = GPS_locxyz)

## Build the data stacks
est_st <- inla.stack(data = GPS_data, A = list(Adata_st),
                     effects = list(GIA = 1:spde_st$n.spde), tag = "est")

## Fix the GPS errors
hyper <- list(prec = list(fixed = TRUE, initial = 0))
formula1 = y ~ -1 +  f(GIA, model = spde_st)
prec_scale <- c(1/GPS_data$std^2)

Run the INLA inference.

## Run INLA
res_st <- inla(formula1, data = inla.stack.data(est_st, spde = spde_st), family = "gaussian",
                 scale =prec_scale, control.family = list(hyper = hyper),
                 control.predictor=list(A=inla.stack.A(est_st), compute =TRUE))

## Extract and project predictions
gps_idx <- inla.stack.index(est_st, tag = "est")$data
vt_idx <- inla.stack.index(est_st, tag = "est")$effects

## GIA
GIA_diff <- INLA_pred$mean
GIA_u <- INLA_pred$sd
proj <- inla.mesh.projector(mesh_st, projection = "longlat", dims = c(360,180), xlim = c(0,359), ylim = c(-89.5, 89.5))
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))),
                       u = as.vector(inla.mesh.project(proj, as.vector(GIA_u))))
GIA_pred$mean <- GIA_pred$diff + ice6g[order_ice6g, ]$trend

ress_st <- list(res_inla = res_st, spde = spde_st, st = est_st, 
            mesh = mesh_st, GIA_pred = GIA_pred)

3.4 Subset model

For the subset model, we first construct the triangulations for the subset of interest. The triangles inheritated from the mesh of the stationary model and those with centroid in side the region of interest are kept in the subset model.

mesh_st <- dt.mesh.addon.posTri(mesh = mesh_st, globe = TRUE) # Add on mesh$posTri: contains the positions of the triangles
Tlonlat <- GEOmap::Lxyz2ll(list(x = mesh_st$posTri[,1], y = mesh_st$posTri[,2], 
                                z = mesh_st$posTri[,3]))# First convert xyz to lonlat
Tlonlat$lon <- ifelse(Tlonlat$lon >=0, Tlonlat$lon, Tlonlat$lon + 359)
mesh_st$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
TinPoly <- unlist(over(zeroPolygon, SpatialPoints(coords=mesh_st$Trill, proj4string = CRS(proj4string(zeroPolygon))), returnList=T))
TAll <- 1:mesh_st$t
ToutPoly <- TAll[-TinPoly]
Omega = dt.Omega(list(TinPoly, 1:mesh_st$t), mesh_st)
## Get the sub mesh
mesh_sub <- mesh.sub(mesh_st, Omega, 2)
plot(mesh_sub, rgl = T)

We generate some pseudo observation along the boundary of the polygons.

## get the boundary of the polygons
boundlines <- as(zeroPolygon, 'SpatialLines') 
obs_bounds <- spsample(boundlines, n = 50, type = "regular") # note points more than specified
proj4string(obs_bounds) <- CRS(proj4string(ice6gsp))
ps_trend <- over(obs_bounds, ice6gsp)$trend
obs_boundsxyz<- do.call(cbind, GEOmap::Lll2xyz(lat = obs_bounds@coords[,2], lon = obs_bounds@coords[,1]))

obs_df1 <- data.frame(y = -ps_trend, ID = "pseudo", lon = obs_bounds@coords[,1], lat = obs_bounds@coords[,2],
                     trend = 0, std = 0.05)
GPS_data2 <- rbind(GPS_data, obs_df1)
GPS_locxyz2 <- rbind(GPS_locxyz, obs_boundsxyz)

Now do the similar buiding of model and inference as before.

## Build the SPDE model
spde_sub <- inla.spde2.matern(mesh = mesh_sub, B.tau = matrix(c(ltau0, -1, 1),1,3), B.kappa = matrix(c(lkappa0, 0, -1), 1,3),
                              theta.prior.mean = c(0,0), theta.prior.prec = c(sqrt(1/theta1_s), sqrt(1/theta2_s)))


## Link the GIA process and GPS data
Adata_sub <- inla.spde.make.A(mesh = mesh_sub, loc = GPS_locxyz2)

## Build the data stacks
est_sub <- inla.stack(data = GPS_data2, A = list(Adata_sub),
                     effects = list(GIA = 1:spde_sub$n.spde), tag = "est")

## Fix the GPS errors
hyper <- list(prec = list(fixed = TRUE, initial = 0))
formula2 = y ~ -1 +  f(GIA, model = spde_sub)
prec_scale <- c(1/GPS_data2$std^2)

Run the INLA inference.

## Run INLA
res_sub <- inla(formula2, data = inla.stack.data(est_sub, spde = spde_sub), family = "gaussian",
                 scale =prec_scale, control.family = list(hyper = hyper),
                 control.predictor=list(A=inla.stack.A(est_sub), compute =TRUE))

INLA_pred <- res_sub$summary.random$GIA

## Extract and project predictions
gps_idx <- inla.stack.index(est_sub, tag = "est")$data
vt_idx <- inla.stack.index(est_sub, tag = "est")$effects

## GIA
GIA_diff <- INLA_pred$mean
GIA_u <- INLA_pred$sd
proj <- inla.mesh.projector(mesh_sub, projection = "longlat", dims = c(360,180), xlim = c(0,359), ylim = c(-89.5, 89.5))
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))),
                       u = as.vector(inla.mesh.project(proj, as.vector(GIA_u))))
GIA_pred$mean <- GIA_pred$diff + ice6g[order_ice6g, ]$trend

ress_sub <- list(res_inla = res_sub, spde = spde_sub, st = est_sub, 
            mesh = mesh_sub, GIA_pred = GIA_pred)

3.5 Partition model

The following build the mesh for the partiontion model. The triangles are in finer resoltuion in the region of interest than those in the zero region. We use the mesh from the subset model for the region of interest. For the zero region, we generate about 100 points (119 to be exact) inside the zero polygons.

points0 <- spsample(zeroPolygon, n = 100, type = "Fibonacci") # in lon lat
points0xyz <- do.call(cbind, GEOmap::Lll2xyz(lat = points0@coords[,2], lon = points0@coords[,1])) 
points1xyz <- mesh_sub$loc
mesh_part <- inla.mesh.2d(loc = rbind(points0xyz, points1xyz), cutoff = 0.01, max.edge = 0.5)

## Separate the mesh by regions
mesh_part <- dt.mesh.addon.posTri(mesh = mesh_part, globe = TRUE)
Tlonlat <- GEOmap::Lxyz2ll(list(x = mesh_part$posTri[,1], y = mesh_part$posTri[,2], z = mesh_part$posTri[,3]))
Tlonlat$lon <- ifelse(Tlonlat$lon >=0, Tlonlat$lon, Tlonlat$lon + 359)
mesh_part$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
TinPoly <- unlist(over(zeroPolygon, SpatialPoints(coords=mesh_part$Trill, proj4string = CRS(proj4string(zeroPolygon))), returnList=T))
TAll <- 1:mesh_part$t
ToutPoly <- TAll[-TinPoly]
Omega = dt.Omega(list(TinPoly, 1:mesh_part$t), mesh_part)
plot(mesh_part, t.sub = Omega[[2]], rgl = TRUE)
plot(mesh_part, t.sub = Omega[[1]], rgl = TRUE, col = "blue", add = TRUE)

We use a constraint partition model where we fix the range in the zero region to be 5 and the variance of the zero region to be the same as the region of interest. Psuedo observations are added in the zero region to further constrain the uncertainty in this region. The following generate the pseudo GPS observations.

points1 <- spsample(zeroPolygon, n = 50, type = "Fibonacci") # in lon lat
points1xyz <- do.call(cbind, GEOmap::Lll2xyz(lat = points1@coords[,2], lon = points1@coords[,1])) 

ps_trend <- over(points1, ice6gsp)$trend

obs_df2 <- data.frame(y = -ps_trend, ID = "pseudo", lon = points1@coords[,1], lat = points1@coords[,2],
                     trend = 0, std = 0.1)
GPS_data3 <- rbind(GPS_data, obs_df2)
GPS_locxyz3 <- rbind(GPS_locxyz, points1xyz)

The following code do the inference.

## Build the SPDE model
Q.mixture <- dt.create.Q(mesh_part, Omega,  fixed.ranges = c(1, NA))

## Transform the parameters for the SPDE_GMRF approximation
prior <- list(sigma = tsigma, range = matrix(trho, ncol = 2))
log.prior <- dt.create.prior.log.norm(prior.param = prior) 
spde_part <- dt.inla.model(Q = Q.mixture, log.prior=log.prior)


## Link the GIA process and GPS data
Adata_part <- inla.spde.make.A(mesh = mesh_part, loc = GPS_locxyz3)

## Build the data stacks
est_part <- inla.stack(data = GPS_data3, A = list(Adata_part),
                     effects = list(GIA = 1:mesh_part$n),
                     remove.unused = FALSE, tag = "est")


## Fix the GPS errors
hyper <- list(prec = list(fixed = TRUE, initial = 0))
formula3 = y ~ -1 +  f(GIA, model = spde_part)
prec_scale <- c(1/GPS_data3$std^2)

Run the INLA inference.

## Run INLA
res_part <- inla(formula3, data = inla.stack.data(est_part), family = "gaussian",
                 scale =prec_scale, control.family = list(hyper = hyper),
                 control.predictor=list(A=inla.stack.A(est_part), compute =TRUE))

INLA_pred <- res_part$summary.random$GIA

## Extract and project predictions
gps_idx <- inla.stack.index(est_sub, tag = "est")$data
vt_idx <- inla.stack.index(est_sub, tag = "est")$effects

## GIA
GIA_diff <- INLA_pred$mean
GIA_u <- INLA_pred$sd
proj <- inla.mesh.projector(mesh_part, projection = "longlat", dims = c(360,180), xlim = c(0,359), ylim = c(-89.5, 89.5))
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))),
                       u = as.vector(inla.mesh.project(proj, as.vector(GIA_u))))
GIA_pred$mean <- GIA_pred$diff + ice6g[order_ice6g, ]$trend

ress_part <- list(res_inla = res_part, spde = spde_part, st = est_part, 
            mesh = mesh_part, GIA_pred = GIA_pred)

3.6 Plot and compare results from the three models

pars_pm1 <- marginal_par(res = ress_st, process = "GIA")
pars_pm2 <- marginal_par(res = ress_sub, process = "GIA")
pars_pm3 <- marginal_par(res = ress_part, partition = TRUE, theta.names = c("sigma", "rho"))

## The prior densities
sr <- 5^seq(-6, -1, 0.1)
dr <- dlnorm(sr, meanlog = trho[1], sdlog=trho[2])

ss <- seq(1, 2.5, 0.05)
ds <- dlnorm(ss, meanlog = tsigma[1], sdlog=tsigma[2])

par(mfrow = c(1,2))
## The posterior distributions
Smar1 <- inla.tmarginal(fun = sqrt, pars_pm1$Vmar)
Smar2 <- inla.tmarginal(fun = sqrt, pars_pm2$Vmar)
Smar3 <- inla.tmarginal(fun = sqrt, pars_pm3$thetamar[[1]])

plot(ds~ss, type = "l", col = "grey", xlim = c(1.3, 2.4), ylim = c(0, 12), 
     ylab = "posterior density", xlab = bquote(bold({sigma})))
polygon(x = c(min(ss), ss, max(ss)), y = c(0, ds, 0), border = "grey", col = "grey")
lines(Smar1, type = "l", lty = 1, lwd = 2)
lines(Smar2, col = 2, lty = 2, lwd = 2)
lines(Smar3, col = 4, lty = 4, lwd = 2)
legend("topright", c("stationary",  "subset", "partition", "prior"), lty = c(1, 2, 4, 1), lwd = c(rep(2, 3), 5), col = c(1,2,4, "grey"))

plot(dr~sr, type = "l",col = "grey",xlim = c(0.06, 0.15), ylim = c(0, 65),
     ylab = "posterior density", xlab = bquote(bold({rho})))
polygon(x = c(min(sr), sr, max(sr)), y = c(0, dr, 0), border = "grey", col = "grey")
lines(pars_pm1$Rmar, type = "l", lty = 1, lwd = 2)
lines(pars_pm2$Rmar, col = 2,lty = 2, lwd = 2)
lines(pars_pm3$thetamar[[2]], lty = 4, col = 4,lwd = 2)

ress_st$GIA_pred$model <- "stationary"
ress_sub$GIA_pred$model <- "subset"
ress_part$GIA_pred$model <- "partition"
pred_all <- do.call(rbind, list(ress_st$GIA_pred,  ress_sub$GIA_pred, ress_part$GIA_pred))
pred_all$model <- factor(pred_all$model, levels = c("stationary", "subset", "partition"))

mean_brk <- c(-7, -5, -3, -1, -0.5,  -0.1,  0.1, 0.5, 1, 3, 5, 10, 16)
pred_all$mean2 <- ifelse(is.na(pred_all$mean), 1e4, pred_all$mean)
pred_all$mean2 <- as.numeric(cut(pred_all$mean2, breaks = c(mean_brk, 1e4)))
colb <- rev(RColorBrewer::brewer.pal(5, "Blues"))
colr <- RColorBrewer::brewer.pal(7, "Reds")
colreg1<- c(colb, "#FFFFFF", colr, "#bebebe")
labs = list(labels = c(mean_brk,"NA"), at = c(seq(0, 12, 1),12.8))    

lattice::levelplot(mean2 ~ lon + lat|model, data = pred_all, aspect = "iso",
                     col.regions= colreg1, at=seq(0,13,1),
                    colorkey = list(col = colreg1, labels = labs ),
                     panel = function(x,y,z,...){
                       lattice::panel.levelplot(x,y,z,...)
                       map2 <- map("world2", interior = FALSE, plot = FALSE)
                       lattice::panel.xyplot(x=map2$x, y=map2$y, type = "l", col = "black")
                     },
                     layout = c(1,3),
                     main = "predicted GIA mean field (mm/yr)", xlab = "longitude", ylab = "latitude")

u_brk <- c(0, 0.2, 0.5, 1,1.5, 2, 2.5, 3, 4, 6)
pred_all$u2 <- ifelse(is.na(pred_all$u), 1e4, pred_all$u)
pred_all$u2 <- as.numeric(cut(pred_all$u2, breaks = c(u_brk, 1e4)))
colreg2 <- c('#fff5eb','#fee6ce','#fdd0a2','#fdae6b','#fd8d3c','#f16913','#d94801','#a63603','#7f2704', '#bebebe')
labs = list(labels = c(u_brk, "NA"), at = c(seq(0, 9), 9.8))

lattice::levelplot(u2 ~ lon + lat|model, data = pred_all, aspect = "iso",
                     col.regions= colreg2, at=seq(0,10,1),
                    colorkey = list(col = colreg2, labels = labs),
                     panel = function(x,y,z,...){
                       lattice::panel.levelplot(x,y,z,...)
                       map2 <- map("world2", interior = FALSE, plot = FALSE)
                       lattice::panel.xyplot(x=map2$x, y=map2$y, type = "l", col = "black")
                     },
                     layout = c(1,3),
                     main = "predicted GIA uncertainties (mm/yr)", xlab = "longitude", ylab = "latitude")

GIA Supplementary

Z Sha

14 February 2018

1 Triangulation on \(S^2\)

1.1 Fibonacci points

1.2 Stationary process

1.3 Subset process

1.4 Partiton process

2 Modelling Non-stationarity

2.1 Subset model

2.2 Partition model

3 Application to GIA

3.1 Build the “Zero Region”

3.1.1 Remove GPS stations in the zero region

3.2 Set up priors for hyperparameters

3.3 Stationary model

3.4 Subset model

3.5 Partition model

3.6 Plot and compare results from the three models