Updating Global GIA
Z Sha
13 November 2017
1 Introduction
In this study, we apply the Bayesian data assimilation model on updatding the global GIA process by using GPS data. We use the soluton from the ICE6G-VM5 model as our prior mean for the GIA process. The GPS data are processed into yearly vertical bedrock movement from a selected global network.
Here we first assume the updating process is stationary on the sphere; then we remove the stationarity constraint by using the information of expert defined pseudo polygons where the GIA should be zero in theroy. Three different approaches are used in modelling the non-stationarity:
Use pseudo-polygon observations — Polygon observation
Use pseudo-polygons to define the subset where GIA can be modelled as stationary process — Subset model
Use psdeuo-polygons to define a mixture Gaussian process where the updating process has different spatial properties inside and outside of the pseudo-polygons — Mixture model
In the following, we demonstrate the all these four approaches and compare the results.
2 Data and prior
First we load the data and prior. Note that the ICE6G solution is given on a 1 degree grid. For any given location, we define the GIA prior to be the value of the grid in which the location falls; hence we set it as a spatial polygon data set.
## Load packages
library(sp); library(INLA); library(GEOmap)
source("functions.R")
source("functions-barriers-dt-models-march2017.R")## Load data
## 1 Load GIA prior
ice6g <- read.table("Z:/WP2-SolidEarth/BHMinputs/GIA/GIA_Pel-6-VM5.txt", header = T)
polycoords <- ice6g[,c(6:13, 6,7)]
plist <- lapply(ice6g$ID,
function(x) Polygons(list(Polygon(cbind(lon = as.numeric(polycoords[x, c(1,3,5,7,9)]),
lat = as.numeric(polycoords[x, c(2,4,6,8,10)])))), ID = x))
Plist <- SpatialPolygons(plist, proj4string = CRS("+proj=longlat"))
## Note that in this data set the grid boundaries are defined on [-0.5, - 359.5] and the centre points on [0, 359].
## Corresponding transformation is needed in the following in geometery operations.
## 2 Load GPS data
GPSV4b <- read.table("Z:/WP2-SolidEarth/BHMinputs/GPS/GPS_v04b.txt", header = T)2.1 Prior setup for the parameters
Setup the priors for the parameters for the Gaussian process. Assume the prior distributions of both \(\rho\) and \(\sigma^2\) are log normal. The prior mean of the correlation length is set to be 500km based on expert prior opinion on the residual process and the mean for the variance is 20 which is about the range of the GIA values. The variances of both distributions are set to be large for a vague enough prior.
## Priors mean and variance for the parameters: rho and sigma
mu_r <- 500/6371
v_r <- (1000/6371)^2
mu_s <- 20
v_s <- 40^2
## Transform the parameters for the SPDE_GMRF approximation
trho <- Tlognorm(mu_r, v_r)
tsigma <- Tlognorm(mu_s, v_s)2.2 Initial mesh
The GIA process will be approximated on a triangular mesh with piece-wise linear basis functions, we will only use the GIA values at the triangle vertices to represent the processes. The following generate a semi-regular mesh on a sphere with approximately 1 degree resolution and we define the GIA prior to be GIA values at the mesh nodes. The semi-regular mesh is generated by using Fibonacci points as intial locations. The Fibonacci points are approximately evenly distributed on the sphere. The goal is to produce a map with approximately one degree resolution, hence we also need a similar resolution for the triangulation. Note that the triangle edge be no larger than the process correlation length. A one degree grid is about \(110 \times 110 km\) which corresponds to a correlation length about \(110/6371 = 0.017\). The number of triangles should be \(f = 180 * 360 = 64800\) corresponding to number of vertices of \(v \approx (f+2)*2/5= 25921\).
## Generate mesh
fibo_points <- fiboSphere(N = 12960, L0 = TRUE)
fibo_points_xyz <- do.call(cbind, Lll2xyz(lat = fibo_points[,2], lon = fibo_points[,1]))
mesh1 <- inla.mesh.2d(loc = fibo_points_xyz, cutoff = 0.01, max.edge = 0.5)
summary(mesh1)##
## Manifold: S2
## CRS: N/A
## Vertices: 25921
## Triangles: 51838
## Boundary segm.: 0
## Interior segm.: 0
## xlim: -0.9999485 1
## ylim: -0.9999871 0.9999871
## zlim: -0.9999614 0.99996143 Global stationary model
We first model the discrepensies process to be globally stationary.
3.1 Mapping GIA prior and detrend GPS
We find the GIA prior at the mesh vertices and the GPS locations. Then the GPS observations used in the model are detrended by the GIA prior mean.
## Find GIA prior mean at mesh vertices
meshLL <- Lxyz2ll(list(x=mesh1$loc[,1], y = mesh1$loc[,2], z = mesh1$loc[,3]))
meshLL$lon <- ifelse(meshLL$lon >= -0.5, meshLL$lon,meshLL$lon + 360)
mesh_sp <- SpatialPoints(data.frame(lon = meshLL$lon, lat = meshLL$lat), proj4string = CRS("+proj=longlat"))
mesh_idx <- over(mesh_sp, Plist)
GIA_prior <- ice6g$trend[mesh_idx]
## Detrend the GPS data
GPS_data <- GPSV4b
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_mu3.2 Inference by INLA
We set up some defualt priors that will be used by INLA in all approaches.
## For the SPDE model prior used in the default setting of INLA, we need the following transformation
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
hyper <- list(prec = list(fixed = TRUE, initial = 0))
formula = y ~ -1 + f(GIA, model = GIA_spde)Then build the SPDE model and link he data to the process.
## Build the SPDE model
mesh <- mesh1
GIA_spde <- inla.spde2.matern(mesh, 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 process to observations and predictions
A_data <- inla.spde.make.A(mesh = mesh, loc = GPS_loc)
A_pred <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc, mesh$loc))
## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=GPS_data$trend0), A = list(A_data),
effects = list(GIA = 1:GIA_spde$n.spde), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(A_pred),
effects = list(GIA=1:GIA_spde$n.spde), tag = "pred")
stGIA <- inla.stack(st.est, st.pred)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.
## Fix the GPS errors
prec_scale <- c(1/GPS_data$std^2, rep(1, nrow(A_pred)))
res_inla <- inla(formula, data = inla.stack.data(stGIA, spde = GIA_spde), family = "gaussian",
scale =prec_scale, control.family = list(hyper = hyper),
control.predictor=list(A=inla.stack.A(stGIA), compute =TRUE))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_data)]
GIA_idx <- pred_idx[-(1:nrow(GPS_data))]
## GPS
GPS_u <- INLA_pred$sd[GPS_idx]
GPS_pred <- data.frame(lon = GPS_data$lon, lat = GPS_data$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))))
ress1 <- list(res_inla = res_inla, spde = GIA_spde, st = stGIA,
mesh = mesh1, GPS_pred = GPS_pred, GIA_pred = GIA_pred)4 Pseudo polygon observations
4.1 Generate the Pseudo polygon
Here we still assume the process is stationary globally but use polygon observations in the pseudo polygon regions defined by experts opinion. More details about choosing and generating the polygons can be found here. The following chunk load the pseudo polygons using the ensemble mean of 8 GIA mode solutions for a given threshold value, say \(0.3\).
## Load the pseudo polygon
zeroPolygon <- readOGR(dsn = "Z:/WP1-BHM/Experiment1b/shapefiles", layer = "zero03")## OGR data source with driver: ESRI Shapefile
## Source: "Z:/WP1-BHM/Experiment1b/shapefiles", layer: "zero03"
## with 1 features
## It has 1 fields## Remove polygons that are too small
zeroPolys <- zeroPolygon@polygons[[1]]@Polygons
polyareas <- sapply(zeroPolys, function(x) x@area)
polyholes <- sapply(zeroPolys, function(x) x@hole)
zeropolys2 <- zeroPolys[polyareas > 200 ]
zeroPoly <- zeroPolygon
zeroPoly@polygons[[1]]@Polygons <- zeropolys24.2 Separate the mesh triangles
Since we have less interest of the process within he zero polygons, we may reduce the resolution of the mesh inside the polygons to save some computational cost. The following chunck generate a mesh that is sparse inside the polygon and the outsid mesh remains the same resolution as previous.
#### Dense points outside the polygons
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_points2 <- fiboSphere(N = 500, L0=TRUE)
pinPoly <- unlist(over(zeroPoly, SpatialPoints(coords = fibo_points2), returnList=T))
fibo_inPoly<- fibo_points2[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]))
mesh2 <- inla.mesh.2d(loc = mesh_points_xyz, cutoff = 0.01, max.edge = 0.5)
mesh2 <- inla.mesh.2d(loc = mesh2$loc, cutoff = 0.01, max.edge = 0.5)
summary(mesh2) # 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## For this new mesh the GIA prior mean for the mesh vertices need to re-calculate.
meshLL <- Lxyz2ll(list(x=mesh2$loc[,1], y = mesh2$loc[,2], z = mesh2$loc[,3]))
meshLL$lon <- ifelse(meshLL$lon >= -0.5, meshLL$lon,meshLL$lon + 360)
mesh_sp <- SpatialPoints(data.frame(lon = meshLL$lon, lat = meshLL$lat), proj4string = CRS("+proj=longlat"))
mesh_idx <- over(mesh_sp, Plist)
GIA_prior2 <- ice6g$trend[mesh_idx]Now separate the triangles by the pseudo-polygons.
mesh2 <- dt.mesh.addon.posTri(mesh = mesh2, globe = TRUE)
Tlonlat <- Lxyz2ll(list(x = mesh2$posTri[,1], y = mesh2$posTri[,2], z = mesh2$posTri[,3]))
Tlonlat$lon <- ifelse(Tlonlat$lon >=0, Tlonlat$lon, Tlonlat$lon + 359)
mesh2$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
TinPoly <- unlist(over(zeroPoly, SpatialPoints(coords=mesh2$Trill), returnList=T))
TAll <- 1:mesh2$t
ToutPoly <- TAll[-TinPoly]
Omega = dt.Omega(list(TinPoly, 1:mesh2$t), mesh2)
plot(mesh2, t.sub = Omega[[2]])
plot(mesh2, t.sub = Omega[[1]])
4.3 Generate pseudo GPS observations
We also add some pseudo observations along the boudaries and insdie the polygons to make smooth transition and restrict the value of in the polygon. The mesh vertices can be used for thes purposes. The values of the pseudo observations are set to be the ice6 values at those locations with opposite signs. The measurement errors are set to be 0.1 which is about the same as the smallest of the real observations.
## Remove GPS data inside the polygon
GPS_inPoly <- unlist(over(zeroPoly, SpatialPoints(coords = cbind(GPSV4b$lon, GPSV4b$lat)), returnList=T))
GPS_All <- 1:nrow(GPS_data)
GPS_outPoly <- GPS_All[-GPS_inPoly]
plot(GPS_data[GPS_outPoly,c("lon", "lat")], pch = "+")
GPS_data2 <- GPS_data[GPS_outPoly,]
GPS_loc2 <- GPS_loc[GPS_outPoly,]
## Find the mesh nodes incide the polygons
Vll <- Lxyz2ll(list(x = mesh2$loc[,1], y = mesh2$loc[,2], z = mesh2$loc[,3]))
Vll$lon <- ifelse(Vll$lon < 0, Vll$lon + 360, Vll$lon)
Vll <- cbind(Vll$lon, Vll$lat)
VinPoly <- unlist(over(zeroPoly, SpatialPoints(coords=Vll), returnList=T))
obs_inpoly <- Vll[VinPoly,]
obs_pseudo <- SpatialPoints(coords = obs_inpoly)
proj4string(obs_pseudo) <- proj4string(Plist)
## Now assemble the new GPS data
pobs_idx <- over(obs_pseudo, Plist)
GIA_pobs <- ice6g$trend[pobs_idx]
nobsb <-nrow(obs_pseudo@coords)
GPS_data3 <- data.frame(ID = rep("pseudo", nobsb), lon = obs_pseudo@coords[,1], lat = obs_pseudo@coords[,2],
trend = rep(0, nobsb), std = rep(0.1, nobsb), trend0 = -GIA_pobs)
GPS_loc3 <- mesh2$loc[VinPoly,]4.4 INLA inference
For the polygon observations we need to linke it as a block to the process.
## Link the process to observations and predictions
mesh <- mesh2
GIA_spde <- inla.spde2.matern(mesh, 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)))
block <- c(1:nrow(GPS_data2), rep(nrow(GPS_data2)+1, nrow(GPS_data3)))
A_data <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc2, GPS_loc3),
block = block, n.block = nrow(GPS_data2) + 1, block.rescale = "count")
A_pred <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc2, GPS_loc3, mesh$loc))
## Create the estimation and prediction stack
GIA_Poly_obs <- mean(GPS_data3$trend0)
st.est <- inla.stack(data = list(y=c(GPS_data2$trend0, GIA_Poly_obs)), A = list(A_data),
effects = list(GIA = 1:GIA_spde$n.spde), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(A_pred),
effects = list(GIA=1:GIA_spde$n.spde), tag = "pred")
stGIA <- inla.stack(st.est, st.pred)
## Fix the GPS errors
prec_scale <- c(1/GPS_data$std^2, 1/0.01^2, rep(1, nrow(A_pred)))Then we run INLA and assemble results same as previously and store the resutsl in .
5 Modelling subset of sphere
In this approach, the process is only modelled on a subset of the sphere. We use the same mesh outside the zero region as above. The following subtract the sub mesh and do the INLA inference with same prior setting.
## Get the submesh
mesh_outPoly <- mesh.sub(mesh2, Omega, 2)
mesh <- mesh_outPoly
GIA_spde <- inla.spde2.matern(mesh, 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 process to observations and predictions
A_data <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc2, GPS_loc3))
A_pred <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc2, GPS_loc3, mesh$loc))
## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=c(GPS_data2$trend0, GPS_data3$trend0)), A = list(A_data),
effects = list(GIA = 1:GIA_spde$n.spde), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(A_pred),
effects = list(GIA=1:GIA_spde$n.spde), tag = "pred")
stGIA <- inla.stack(st.est, st.pred)
## Fix the GPS errors
prec_scale <- c(1/GPS_data2$std^2, 1/GPS_data3$std^2, rep(1, nrow(A_pred)))Then we run INLA and assemble results same as previously and store the resutsl in .
6 Mixture model
In this approach we model the GIA process as a mixture of two stationary processes on the zero region and region of intrest. Again the mesh is same as above but we build need to build first build the spde model for the non-stationary process.
mesh <- mesh2
Q.mixture = dt.create.Q(mesh, Omega, fixed.ranges = c(0.01, NA))
log.prior = dt.create.prior.log.exp(prior.param = c(1,1))
GIA_spde = dt.inla.model(Q = Q.mixture, log.prior=log.prior)
## Link the process to observations and predictions
A_data <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc2, GPS_loc3))
A_pred <- inla.spde.make.A(mesh = mesh, loc = rbind(GPS_loc2, GPS_loc3, mesh$loc))
## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=c(GPS_data2$trend0, GPS_data3$trend0)), A = list(A_data),
effects = list(GIA = 1:mesh$n), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(A_pred),
effects = list(GIA=1:mesh$n), tag = "pred")
stGIA <- inla.stack(st.est, st.pred)
## Fix the GPS errors
prec_scale <- c(1/GPS_data2$std^2, 1/GPS_data3$std^2, rep(1, nrow(A_pred)))Then we run INLA and assemble results same as previously and store the resutsl in .
7 Compare results
Plot the posteriors of the hyper parameters.
pars_pm1 <- marginal_par(res = ress1)
pars_pm2 <- marginal_par(res = ress2)
pars_pm3 <- marginal_par(res = ress3)
pars_pm4 <- marginal_par(res = ress4, mixture = TRUE)
## The posterior modes
print(paste("The estimated correlation lengths are:",
pars_pm1$rho_mode*6371,
pars_pm2$rho_mode*6371,
pars_pm3$rho_mode*6371,
pars_pm4$rho_mode*6371, sep = " "))## [1] "The estimated correlation lengths are: 634.522472892937 574.499117174802 457.746154128337 303.53569332672"print(paste("The estimated marginal variances are:",
pars_pm1$sigma_mode,
pars_pm2$sigma_mode,
pars_pm3$sigma_mode,
pars_pm4$sigma_mode, sep = " "))## [1] "The estimated marginal variances are: 2.84908127417138 3.76924354258209 2.19338089855634 1.10537485937626"par(mfrow = c(1,2))
## The posterior distributions
plot(pars_pm1$Vmar, type = "l", xlim = c(1, 3.5), ylim = c(0, 17),
ylab = "posterior density", xlab = bquote(bold({sigma^2})))
lines(pars_pm3$Vmar, col = 2)
lines(pars_pm4$Vmar, col = 4)
#lines(pars_pm2$Vmar, col = 4)
legend("topright", c("stationary", "subset", "mixture"), lty = rep(1), col = c(1,2,4))
plot(pars_pm1$Rmar, type = "l", xlim = c(0.04, 0.12), ylim = c(0, 180),
ylab = "posterior density", xlab = bquote(bold({rho})))
lines(pars_pm3$Rmar, col = 2)
lines(pars_pm4$Rmar, col = 4)
#lines(pars_pm2$Rmar, col = 4)
legend("topright", c("stationary", "subset", "mixture"), lty = rep(1), col = c(1,2,4))
Plot the predictions.
ress1$GIA_pred$model <- "stationary"
ress2$GIA_pred$model <- "polygon_Obs"
ress3$GIA_pred$model <- "subset"
ress4$GIA_pred$model <- "mixture"
pred_all <- do.call(rbind, list(ress1$GIA_pred, ress3$GIA_pred, ress4$GIA_pred))
pred_all$model <- factor(pred_all$model, levels = c("stationary", "subset", "mixture"))
## The predicted mean
level_GIA(pred_data = pred_all, var = "mean", style ="zero", layout = c(1, 3), title = "predicted GIA mean field (mm/yr)")
## The predicted uncertainties
level_GIA(pred_data = pred_all, var = "u", style ="positive",layout = c(1, 3), title = "predicted GIA uncertanties (mm/yr)")
## The difference between the predicted mean and prior mean
level_GIA(pred_data = pred_all, var = "diff", style ="zero", layout = c(1, 3), title = "predicted GIA mean field - prior mean (mm/yr)")