Experiment 2c partition process model (c)
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, coastlines and ocean domains, and use a parameter partition model. In particular, we set the correlation at coast line domain to be near zero so that the land and ocean are separate. By using the parameter partition model we can avoid boundary effect. The best way would be allowing the correlaty decay smoothly from center to boundary and we will try to build up this in the next step.
2 Data and priors
Load the data and do some pre-processing and formatting. All the priors needed in the folowing are also set up in the scripts. Some of them are learned from the data and some are given by experts’ opinion.
source(paste0(wd,"experiments/Experiment2/Doc/c/data_exp2c_c.R"))## OGR data source with driver: ESRI Shapefile
## Source: "Z:/WP1-BHM/maps/ne_110m_land", layer: "ne_110m_land"
## with 127 features
## It has 2 fields
## [1] "The vlm adjustment is -0.217233626677839 mm/yr"
## File Z:/WP2-SolidEarth/Bramha/BHMinput/GRACE_trends_Oceanadded_Bramha.nc (NC_FORMAT_CLASSIC):
##
## 4 variables (excluding dimension variables):
## double Lat[r,c]
## double Long[r,c]
## double trend[r,c]
## double uncert[r,c]
##
## 2 dimensions:
## r Size:41168
## c Size:1
3 Generate triangulations for the SPDE approximation
Create triangulation and label the triangles as Land, buffer and Ocean.
## Genereate Fibonacci points on the sphere
fibo_points <- fiboSphere(N = 3e4, LL = FALSE)
mesh0 <- inla.mesh.2d(loc = fibo_points, cutoff = 0.01, max.edge = 1)
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.9999667Partition the triangulation.
mesh0 <- dt.mesh.addon.posTri(mesh = mesh0, globe = TRUE)
Tlonlat <- Lxyz2ll(list(x = mesh0$posTri[,1], y = mesh0$posTri[,2], z = mesh0$posTri[,3]))
mesh0$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
T_sp <- SpatialPoints(coords=mesh0$Trill, proj4string = CRS(proj4string(Land)))
## Devide these points into three groups ( > +- 85)
Npole_id <- which(T_sp@coords[,2] > 85)
Spole_id <- which(T_sp@coords[,2] < -85)
mid_id <- (1:nrow(T_sp@coords))[-c(Npole_id, Spole_id)]
Spole_sp <- T_sp[Spole_id]
mid_sp <- T_sp[mid_id]
## Nothing to be done on the north pole region points since they will all be classified in the Ocean.
## Transform the middle part and South pole region into suitable projection system
mid_sp <- spTransform(mid_sp, CRS("+init=epsg:3857"))
Spole_sp <- spTransform(Spole_sp, CRS("+init=epsg:5481"))
## Create the Land label
TinLand1 <- mid_id[unlist(over(Land1, mid_sp, returnList=T))]
TinLand2 <- Spole_id[unlist(over(Land2, Spole_sp, returnList=T))]
TinLand <- c(TinLand1, TinLand2)
plot(T_sp[TinLand], pch = ".")
## Create the independent zone label
TinID1 <- mid_id[unlist(over(Indzone1, mid_sp, returnList=T))]
TinID2 <- mid_id[unlist(over(Indzone2, Spole_sp, returnList=T))]
TinID <- c(TinID1, TinID2)
plot(T_sp[TinID], pch = ".")
## Finally the Ocean label
TinOcean <- (1:nrow(T_sp@coords))[-c(TinLand, TinID)]
plot(T_sp[TinOcean], pch = ".", col = 2, add = T)
## check plot
plot(T_sp, pch = ".")
plot(T_sp[TinOcean], pch = ".", col = 2, add = T)
plot(T_sp[TinLand], pch = ".", col = 3, add = T)
plot(T_sp[TinID], pch = ".", col = 4, add = T)
Omega = dt.Omega(list(TinOcean,TinLand, TinID), mesh0)
mesh_ocean <- mesh.sub(mesh0, Omega, 1)
mesh_land <- mesh.sub(mesh0, Omega, 2)
mesh_ind <- mesh.sub(mesh0, Omega, 3)
# plot(mesh0, t.sub = Omega[[1]], rgl = TRUE, col = "yellow")
# plot(mesh0, t.sub = Omega[[2]], rgl = TRUE, col = "grey", add = TRUE)
# plot(mesh0, t.sub = Omega[[3]], rgl = TRUE, col = "blue", add = TRUE)
#play3d(spin3d(axis = c(1,1,1)), duration = 8)
#movie3d( spin3d(axis = c(1,1,1)), duration = 8 )
## mesh for steric
TinLands <- unlist(over(Land, SpatialPoints(coords=mesh0$Trill, proj4string = CRS(proj4string(Land))), returnList=T))
TAll <- 1:mesh0$t
TinOceans <- TAll[-TinLands]
Omega_s = dt.Omega(list(TinOceans, 1:mesh0$t), mesh0)
mesh_ocean_s <- mesh.sub(mesh0, Omega_s, 1)4 Build SPDE approximation for the processes
Now we can build the SPDE approximations for the two processes.
massQ <- ppb.create.Q(mesh= mesh0, subdomain = Omega, sphere = TRUE,
initial.theta = c(lrhog1, lrhog2, lsigmag1, lsigmag2),
theta_buffer = c(log(0.02), log(10)))
## Transform the parameters for the SPDE_GMRF approximation
prior <- list(sigma = rbind(tsigmag1, tsigmag2), rho = rbind(trhog1, trhog2))
log.prior <- ppb.create.priors2(prior.param = prior)
mass_spde <- ppb.inla.model(Q = massQ, log.prior=log.prior)
## steric
steric_spde <- inla.spde2.matern(mesh_ocean_s, B.tau = matrix(c(ltau2, -1, 1),1,3), B.kappa = matrix(c(lkappa2, 0, -1), 1,3),
theta.prior.mean = c(0,0), theta.prior.prec = c(sqrt(1/theta2_ss), sqrt(1/theta2_rs)))5 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)\]
The following chunk builds the matrices and stacks for the above model.
## The data
ygrace <- grace_df$trendgia
yaltv<- alt_data2$trend_vlm
## The A matrice for linking observation and processes
## Link GRACE to mass
A_grace <- inla.spde.make.A(mesh = mesh0, loc = grace_loc)
## Link Altimetry and grace to steric and mass
A_altv_s <- inla.spde.make.A(mesh = mesh_ocean_s, loc = alt_loc)
A_altv_m <- inla.spde.make.A(mesh = mesh0, loc = alt_loc)
## The A matrices for predict ssh
A_ssh_s <- A_ssh_m <- Diagonal(n= steric_spde$n.spde, x =1)
A_ssh_m <- inla.spde.make.A(mesh = mesh0, loc = mesh_ocean_s$loc)
## The errors
prec_scale <- c(1/(3*grace_df$std)^2, 1/alt_data2$err_ssh^2, rep(1, nrow(A_ssh_s)))
## Create the mass and steric stacks
stkmass <- inla.stack(data = list(y=ygrace), A = list(A_grace),
effects = list(mass = 1:(mesh0$n)), tag = "mass")
stksteric <- inla.stack(data = list(y=yaltv), A = list(A_altv_m, A_altv_s),
effects = list(list(mass = 1:mesh0$n),
list(steric = 1:steric_spde$n.spde)), tag = "steric")
## Create the stack for predicting SSH
stkssh <- inla.stack(data = list(y=NA), A = list(A_ssh_m, A_ssh_s),
effects = list(list(mass = 1:mesh0$n),
list(steric = 1:steric_spde$n.spde)),
remove.unused = FALSE, tag = "ssh")
## join the two stacks
stkall <- inla.stack(stkmass, stksteric, stkssh, remove.unused = FALSE)6 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
## constraint 1 -- vertices sum to zero
A1 <- matrix(1, nrow = 1, ncol = mesh0$n)
## Constraint 2 -- grace loc sum to zero
A2 <- matrix(colSums(A_grace), nrow = 1)
## Constraint 3 -- longlat grid sum to zero
gridll <- ice6g_vlm@coords
gridxyz <- do.call(cbind, Lll2xyz(lat = gridll[,2], lon = gridll[,1]))
A_grid <- inla.spde.make.A(mesh = mesh0, loc = gridxyz)
weights <- ice6g_vlm$areas/sum(ice6g_vlm$areas)
A3 <- matrix(weights, nrow = 1) %*% A_grid
A <- as.matrix(rbind(A1, A2, A3))
## The formular
formula = y ~ -1 + f(mass, model = mass_spde,
extraconstr = list(A = A, e = c(0, 0, 0))) + 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.predictor=list(A=inla.stack.A(stkall), compute =TRUE), verbose = TRUE)7 Results
7.1 Assemble and save results
Now assemble the inla inference and prediction results for the result from .
7.2 Plot the posteriors of the hyper parameters
## [1] "The estimated correlation length for mass in Ocean is: 493.532810940697 km"## [1] "The estimated marginal variance for mass in Ocean is: 0.627669637836871 mm/yr"## [1] "The estimated correlation length for mass in Land is: 167.374392904312 km"## [1] "The estimated marginal variance for mass in Land is: 2.30334588447717 mm/yr"## [1] "The estimated correlation length for steric is: 1177.63890117039 km"## [1] "The estimated marginal variance for steric is: 4.75269047402543 mm/yr"7.3 Plot the predictions
7.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.
## OGR data source with driver: ESRI Shapefile
## Source: "Z:/WP1-BHM/maps/ne_110m_ocean", layer: "ne_110m_ocean"
## with 2 features
## It has 3 fields## GRACE predicted trend mean over the corresponding regions:
##
## Regions | SA sum to zero
## ---------------------------------------
## Global | 0.0008242815
## Ocean | 1.182118
## Land | -2.904424