Experiment 2c (b): parameter partition model
Z Sha
19 April 2018
1 Introduction
Now we allow dependency between land and ocean and m
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 (note the an ocean model is added back to the GRACE 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.
Now we do some sanity check on the GRACE and GIA data. First the GRACE and GIA data should both sum to zero. Then we remove the GIA signal from GRACe and check that also sum to zero.
## GRACE sum to zero -- equal area -- simple average
mean(grace_bdv$trend)## [1] 0.05943476## GIA_ewh sum to zero
gia_sp <- gia_ewh
coordinates(gia_sp) <- c("lon", "lat")
gridded(gia_sp) <- TRUE
gia_sp$areas <- geosphere::areaPolygon(as(gia_sp, "SpatialPolygons"))/(1000^2)
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
proj4string(gia_sp)<- proj4string(grace_bdv) <- CRS("+proj=longlat")
gridded(gia_sp) <- TRUE
giaG <- over(grace_bdv, as(gia_sp, "SpatialPolygons"))
grace_df$trendgia <- grace_bdv$trendgia <- grace_bdv$trend - gia_sp$trend[giaG]
mean(grace_bdv$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.
grace_coords <- grace_bdv@coords
grace_coords[,1] <- ifelse(grace_coords[,1]> 180, grace_coords[,1]-360, grace_coords[,1])
glid <- unlist(over(newLand, SpatialPoints(coords= grace_coords, proj4string = CRS(proj4string(newLand))), returnList = TRUE))
goid <- (1:nrow(grace_coords))[-glid]
v_ocean <- gstat::variogram(trendgia ~ 1, data = grace_bdv[goid,])
v_land <- gstat::variogram(trendgia ~ 1, data = grace_bdv[-goid,])
plot(v_ocean)
plot(v_land)
# -- the range for ocean is likely to be about 1200km
# -- the variance is likely to be around (9mm/yr)^2
# -- the range for land is likely to be 600km
# -- the variance is (19mm/yr)^2For 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.
## Priors for mass
# -- Beta distribution for the range
# -- set the mode and let beta = 2*alpha
# -- mode = (alpha - 1)/(alpha + beta -2)
m_rg1 <- 1200/6371 # mode of rho_ocean
m_rg2 <- 600/6371 # mode of rho_land
# -- transform to beta parameters
Tbeta <- function(mode, beta.s){
alpha <- (2*mode-1)/(mode*(1+beta.s)-1)
beta <- beta.s*alpha
return(c(alpha, beta))
}
trhog1 <- Tbeta(mode = m_rg1, beta.s = 3)
trhog2 <- Tbeta(mode = m_rg2, beta.s = 6)
## Plot the prior
xx <- seq(0,1, 0.01)
yy1 <- dbeta(xx, trhog1[1], trhog1[2])
yy2 <- dbeta(xx, trhog2[1], trhog2[2])
plot(xx, yy1, type = "l")
plot(xx, yy2, type = "l")
# -- initial values for rho_ocean and rho_land
lrhog1 <- log(m_rg1)
lrhog2 <- log(m_rg2)
## -- Gamma distribution for the variance
## -- set mean and variance
## -- mean = alpha/beta, variance = alpha/beta^2
m_sg1 <- 9 # mean of sigma_ocean
v_sg1 <- (9)^2/2 # variance of sigma_ocean
m_sg2 <- 19 # mean of sigma_land
v_sg2 <- (19)^2/2 # variance of sigma_land
TGamma <- function(mean, variance){
beta <- mean/variance
alpha <- mean^2/variance
return(c(alpha, beta))
}
tsigmag1 <- TGamma(mean = m_sg1, variance = v_sg1)
tsigmag2 <- TGamma(mean = m_sg2, variance = v_sg2)
xx <- seq(0, 50, 1)
yy3 <- dgamma(xx, shape = tsigmag1[1], rate = tsigmag1[2])
yy4 <- dgamma(xx, shape = tsigmag2[1], rate = tsigmag2[2])
plot(xx, yy3, type = "l")
plot(xx, yy4, type = "l")
#-- initial values for sigma_Ocean and sigma_land
lsigmag1 <- log(m_sg1)
lsigmag2 <- log(m_sg2)
## Priors for steric
mu_r2 <- 2000/6371 # mean of rho
v_r2 <- 1 # vague variace for rho
mu_s2 <- 8 # mean of sigma
v_s2 <- 16^2 # vaue variance for sigma
## -- steric
trho2 <- Tlognorm(mu_r2, v_r2)
tsigma2 <- Tlognorm(mu_s2, v_s2)
lsigma2 <- tsigma2[1]
theta2_ss <- tsigma2[2]
lrho2 <- trho2[1]
theta2_rs <- trho2[2]
lkappa2 <- log(8)/2 - lrho2
ltau2 <- 0.5*log(1/(4*pi)) - lsigma2 - lkappa24 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. In general, the triangles in the mesh should have similar size and shape and the size of the triangles should be no larger than the spatial correlation length.
In our exploratory analysis, the correlation length is quite long, therefore if the purpose is estimating the hyper-parameters, then we can choose to have a sparse mesh with resolution no large than the correlation length. Here we could use say a triangle size to have approximately 5 degree resolution so that the edge length of the triangles are smaller than 500km on average which is much smaller than 2000km.
However, Since our purpose is to produce a map of the process at high resolution, we must have a mesh with at least the desired prediction resolution. For a 1 degree resolution map, we will need about \(360 \times 180 = 64800\) triangles corresponding to number of vertices of \(v \approx (f+2)*2/5= 25921\). When using the Fibonacci grid points in long lat coordinates as the starting points for generating the mesh, the number of points is about half of the triangle vertices.
## 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.9999667The above triangulation can be used for the mass process; however since the steric process is only defined on the ocean, we need to remove the triangles in the land regions.
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)
TinLand <- unlist(over(newLand, SpatialPoints(coords=mesh0$Trill, proj4string = CRS(proj4string(newLand))), returnList=T))
TAll <- 1:mesh0$t
TinOcean <- TAll[-TinLand]
Omega = dt.Omega(list(TinOcean, 1:mesh0$t), mesh0)
mesh_ocean <- mesh.sub(mesh0, Omega, 1)
mesh_land <- mesh.sub(mesh0, Omega, 2)
summary(mesh_ocean)##
## Manifold: S2
## CRS: N/A
## Vertices: 22622
## Triangles: 42855
## 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: 9822
## Triangles: 17141
## Boundary segm.: 0
## Interior segm.: 0
## xlim: -0.85839 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.
## mass
massQ <- pp.create.Q(mesh = mesh0, subdomain = Omega, sphere = TRUE, initial.theta = c(lrhog1, lrhog2, lsigmag1, lsigmag2))
## Transform the parameters for the SPDE_GMRF approximation
prior <- list(sigma = rbind(tsigmag1, tsigmag2), rho = rbind(trhog1, trhog2))
log.prior <- pp.create.priors2(prior.param = prior)
mass_spde <- pp.inla.model(Q = massQ, log.prior=log.prior)
## steric
steric_spde <- inla.spde2.matern(mesh_ocean, 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)))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)\]
The following chunk builds the matrices and stacks for the above model.
## The data
ygrace <- grace_bdv$trendgia
yaltv<- alt_data2$trend_vlm
## The A matrice for linking observation and processes
## Link GRACE to mass
grace_loc <- do.call(cbind, Lll2xyz(lat = grace_bdv@coords[,2], lon = grace_bdv@coords[,1]))
A_grace_m <- 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, loc = alt_loc)
A_altv_m <- inla.spde.make.A(mesh = mesh0, loc = alt_loc)
## The A matrices for predict ssh
A_ssh_s <- Diagonal(n= steric_spde$n.spde, x =1)
A_ssh_m <- inla.spde.make.A(mesh = mesh0, loc = mesh_ocean$loc)
## The errors
prec_scale <- c(1/grace_bdv$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_m),
effects = list(list(mass = 1:mesh0$n)),
remove.unused = FALSE, 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)),
remove.unused = FALSE, 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)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
## 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_m), 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)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: 522.431912169461 km"## [1] "The estimated marginal variance for mass in Ocean is: 7.03656082037188 mm/yr"## [1] "The estimated correlation length for mass in Land is: 109.239978019774 km"## [1] "The estimated marginal variance for mass in Land is: 34.5315043855623 mm/yr"## [1] "The estimated correlation length for steric is: 306.242518844074 km"## [1] "The estimated marginal variance for steric is: 7.81002223660618 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.
## 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.007110702
## Ocean | 1.090186
## Land | -2.705778