Experiment 2c (a): process partition model
Z Sha
19 April 2018
1 Introduction
We following the same setting in the previous experiment 2 series. In this version, we try to build a non-stationary model for the mass process. We start from partion the process into land and ocean domains and use a process partition model that assumes there is no spatial connection at the boundaries.
2 Data
Load the pre-processed altimetry data, GRACE data, the ocean land mask and the ICE6G-VM5 data. The data are the linear year trend estimates over the period 2005-2015 with their corresponding estimated standard errors.
The Altimetry data used here is the CCI given by WL at 1 degree long-lat grid.
We use BDV’s GRACE equiarea data. The GIA signal has not been removed, so we use the ICE-6G GIA solution for the correction. The ICE-6G solution is converted to the unit of mm/yr in equivalent water height by BDV.
The ocean mask is downloaded from Nature Earth at 110m resolution. This will be used for calculating the Ocean areas separating the triangulations in the ocean and land.
The ICE6G-VM5 data are trend of vertical land motions processed by MS and they are at 1 degree resolution. The data is used to calculated the effect of vlm averaged over the ocean.
## Load and process data
source(paste0(wd, "experiments/Experiment2/Doc/c/Rscripts/exp2c_loadData.R"))After loading the data, we do some sanity check on the GRACE and GIA data. First show the global vlm adjustment. The we check that the GRACE and GIA data should both sum to zero. Finally we check that the GRACE after removing GIA also sum to zero.
print(paste("The vlm adjustment is ", vlm, "mm/yr"))## [1] "The vlm adjustment is -0.220296832981777 mm/yr"## GRACE sum to zero -- equal area -- simple average
mean(grace_sp$trend)## [1] 0.05943476## GIA_ewh sum to zero
sum(gia_sp$trend*gia_sp$areas)/sum(gia_sp$areas)## [1] 0.004798744## Remove GIA from GRACE -- use grace - gia grid value which grace falls in
mean(grace_sp$trendgia)## [1] 0.05368853 Set up priors
Now we do some exploratory analysis which will be helpful in setting up the prior for hyper-parameters \(\rho\) and \(\sigma^2\) for the mass and steric processes.
For the mass process, We learn the initial values from the variogram of the the GRACE data. We learn the parameters for land ocean separately.
For the steric process, we do not have direct observations, but we can coarsely learn the parameters from the residuals of altimetry minus GRACE. This is the same as the previous.
4 Generate triangulations for the SPDE approximation
To build the spde model used for approximating the process, we need to generate a triangular mesh for the process. We have about 60,000 triangles with roughly the same size. The resolution is high enough for the approximation.
Based on this, we separate the triangles in Ocean and land. The Ocean trianglulation will also be used for steric.
summary(mesh0)##
## 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.9999667summary(mesh_ocean)##
## Manifold: S2
## CRS: N/A
## Vertices: 22615
## Triangles: 42695
## Boundary segm.: 0
## Interior segm.: 0
## xlim: -0.999934 0.999934
## ylim: -0.9999457 0.9998971
## zlim: -0.9968333 0.9999667summary(mesh_land)##
## Manifold: S2
## CRS: N/A
## Vertices: 10066
## Triangles: 17301
## Boundary segm.: 0
## Interior segm.: 0
## xlim: -0.9621331 0.9955881
## ylim: -0.9889051 0.9924408
## zlim: -0.9999667 0.99435 Build SPDE approximation for the processes
Now we can build the SPDE approximations for the two processes.
5.1 Separate mass processes
In this first approach, we treat mass as sum of two independent sub processes. This is essentially the same as if we use a process partition model but we will not be able to use the sum to zero constraints for two independent processes in INLA.
6 Link data to process and build stacks
Now we link the data to the process. The data we have and the triangulation are about the same resolution, so we use point to point map here. Denote by \(y_{grace}\) the grace observations for updating the mass process and \(y_{altv}\) the altimetry observations (vlm adjusted) for the sum of mass and steric. Then the linear model can be stacked as \[\left(\begin{array}{c} y_{grace} \\ y_{altv} \end{array}\right) = \left(\begin{array}{cc} A_{grace} & 0 \\ A_{altv1} & A_{altv2} \end{array}\right) \left(\begin{array}{c} X_{mass} \\ X_{steric} \end{array}\right) + \left(\begin{array}{c} e_{grace} \\ e_{altv} \end{array}\right)\]
Similarly, we can still use the form \(y = Ax + e\), with \[y = (y_{grace}, y_{altv})^T, \; X = (X_{mass}, X_{steric})^T, \; e = (e_{grace}, e_{altv})^T\\ A = \left(\begin{array}{cc} A_{grace} & 0 \\ A_{altv1} & A_{altv2} \end{array}\right)\]
7 INLA inference
Now we can run INLA for the Bayesian inference. Do not run on a desktop the process may use up to 64GB memory at peak. We ran this on a server with enough memory.
## Fix altimetry errors as they are known
hyper <- list(prec = list(fixed = TRUE, initial = 0))
## The formular -- we add the constraint that mass change sum to zero
formula = y ~ -1 + f(massO, model = massOcean_spde) + f(massL, model = massLand_spde) +
f(steric, model = steric_spde)
## Run INLA
res_inla <- inla(formula, data = inla.stack.data(stkall), family = "gaussian",
scale =prec_scale,
control.family = list(hyper = hyper), control.compute = list(config = TRUE),
control.predictor=list(A=inla.stack.A(stkall), compute =TRUE),
verbose = TRUE)
saveRDS(res_inla, file = "/./projects/GlobalMass/WP1-BHM/Experiment2c/Exp2c_a.rds")8 Results
8.1 Assemble and save results
Now assemble the inla inference and prediction results for the result from .
8.2 Plot the posteriors of the hyper parameters
## [1] "The estimated correlation length for mass in Ocean is: 722.877336323702 km"## [1] "The estimated marginal standard error for mass in Ocean is: 5.39649549288776 mm/yr"## [1] "The estimated correlation length for mass in Land is: 428.378391025416 km"## [1] "The estimated marginal standard error for mass in Land is: 21.350475875856 mm/yr"## [1] "The estimated correlation length for steric is: 757.437816358629 km"## [1] "The estimated marginal standard error for steric is: 5.70619247624516 mm/yr"8.3 Plot the predictions
8.4 Sanity checks on the updated mass
To check whether the updated mass is reasonable, we calculate the mass trend averaged over Earth sphere, ocean, and land. 
## GRACE predicted trend mean over the corresponding regions:
##
## Regions | SA sum to zero
## ---------------------------------------
## Global | 0.2773815
## Ocean | 1.222844
## Land | -2.0478699 Predict steric on Basin level
To compare our prediction with other existing work, we also predict the averge steric change at the the basin level. This can be done either before or after calling the INLA procedure. Here we use both and also compare the results. The uncertainties will be different due to non-zero correlation.
load(paste0(dd, "WP1-BHM/maps/Ocean/basins.rda"))
pred_basin <- ress_2c_a$pred$basins
allbasinsdf$trend1 <- pred_basin$mean
allbasinsdf$trend2 <- pred_basin$meanpost
allbasinsdf$u1 <- pred_basin$u
allbasinsdf$u2 <- pred_basin$u2
allbasinsdf$u3 <- pred_basin$upost
spplot(allbasinsdf, "trend1", at = c(-2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5), col.regions = rev(c('#d73027','#f46d43','#fdae61','#fee090','#ffffbf','#e0f3f8','#abd9e9','#74add1','#4575b4')),
main = "Basin trend (posterior mean)")
spplot(allbasinsdf, "trend2", at = c(-2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5), col.regions = rev(c('#d73027','#f46d43','#fdae61','#fee090','#ffffbf','#e0f3f8','#abd9e9','#74add1','#4575b4')),
main = "Basin trend (posterior sample mean)")
spplot(allbasinsdf, "u1", at = c(0, 0.04, 0.08, 0.12, 0.16, 0.2),
col.regions = c('#fef0d9','#fdd49e','#fdbb84','#fc8d59','#ef6548','#d7301f','#990000'),
main = "uncertainty (posterior covariance approximation)")
spplot(allbasinsdf, "u2", at = c(0, 0.04, 0.08, 0.12, 0.16, 0.2),
col.regions = c('#fef0d9','#fdd49e','#fdbb84','#fc8d59','#ef6548','#d7301f','#990000'),
main = "uncertainty (independent errors)")
spplot(allbasinsdf, "u3", at = c(0, 0.04, 0.08, 0.12, 0.16, 0.2),
col.regions = c('#fef0d9','#fdd49e','#fdbb84','#fc8d59','#ef6548','#d7301f','#990000'),
main = "uncertainty (posterior sample approximation)")