Experiment 2a – Update Mass

1 Introduction

2 Update mass

Now we do the same thing for the mass process using the GRACE data.

2.1 Load data

First we load the GRACE data.

grace_data <- read.table(paste0(dd, "WP2-SolidEarth/BHMinputs/GRACE/GRACE_v02.3b_trends_v03.txt"), header = T)
grace_loc <-  read.table(paste0(dd, "WP2-SolidEarth/BHMinputs/GRACE/GRACE_v02.3b_loc_v03.txt"), skip = 1)
n_grace <- ncol(grace_loc)
n_ll <- (n_grace-4)/2
names(grace_loc) <- c("id", "area", "lon_c", "lat_c", paste0(c("lon", "lat"), rep(1:n_ll, each = 2)))

## Create spatial polygons data frame
Polygon_list <- list()
for(i in 1:nrow(grace_loc)){
  lons <- na.omit(as.numeric(c(grace_loc[i, seq(5,ncol(grace_loc), 2)], grace_loc[i, 5])))
  lats <- na.omit(as.numeric(c(grace_loc[i, seq(6,ncol(grace_loc), 2)], grace_loc[i, 6])))
  Polygon_list[[i]] <- Polygon(cbind(lons, lats))
}

Polygons_list <- lapply(1:length(Polygon_list), function(x) Polygons(list(Polygon_list[[x]]), x))
SpPolygon <- SpatialPolygons(Polygons_list, proj4string = CRS("+proj=longlat"))

grace_sp <- SpatialPolygonsDataFrame(SpPolygon,grace_data)
grace_sp$mmweqt <- ifelse(abs(grace_sp$mmweq) > 20, sign(grace_sp$mmweq)*20, grace_sp$mmweq )
spplot(grace_sp, "mmweqt", at = seq(-21, 21, 2), main = "The GRACE super mascons data (mm/yr ewh)")

## note the area calcuated by R is different from those given in the dataset -- ask Maike

2.2 prior set up

Use the polygon data to have rough estimate of the correlation length and variance.

grace_v <- data.frame(mean = grace_data$mmweq, lon = grace_loc$lon_c, lat = grace_loc$lat_c)
coordinates(grace_v) <- c("lon", "lat")
proj4string(grace_v) <- CRS("+proj=longlat")
v2 <- variogram(mean~1, grace_v) 
plot(v2)

## Priors mean and variance for the parameters: rho and sigma
mu_r <- 2500/6371
v_r <- 1
mu_s <- 10 # area scale to about 100km^2 ~ 1 degree resolution
v_s <- 20^2

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

## Build the SPDE model with the prior
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

M_spde <- inla.spde2.matern(mesh0, 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)))

2.3 Generate mesh

We also need a mesh to represent the mass process. For the same reason as discussed in previous section for SSH, we need a mesh with 1 degree resolution. The GRACE data is all over the globe, so we can use the .

2.4 Link data to process

The GRACE data are mass anomalies aggregated into large polygon grid and we call them “super mascons” as they are on larger grid than the original data provided by the GRACE project. The area of the super mascons are provided in the original data. We calculate the areas by using function in R and compare them to the provided ones. There are some difference and need to be check with MS later.

grace_area <- geosphere::areaPolygon(grace_sp)/(1000^2) # km^2
plot(grace_loc$area ~ grace_area, xlim = c(0, 3e5), ylim = c(0, 3e5), main = "R area vs data area")
abline(a = 0, b = 1)

grace_sp$area <-grace_area
spplot(grace_sp, "area", main ="Areas of the grid calculated by R (km^2)")

grace_sp$area2 <- grace_loc$area
spplot(grace_sp, "area2", main = "Areas of the grid given in the data (km^2)")

area_mean <- mean(grace_area)

To link the data to the process, we need to integrate the process value over the corresponding polygon. The integration can be done numerically as weighted sum of the grid values within the polygons. The grid values are represented by the GMRF.

## For each polygon observation we generate the regular spaced grid and the number of grid cell is proportional to the area of the polygon

## Generate the integration grid for each polygons
poly_block <- function(i, dis = 10){
  sp_i <- SpatialPolygons(list(grace_sp@polygons[[i]]), proj4string=CRS("+proj=longlat"))
  area_i <- grace_sp$area[i]
  grid_i <- spsample(sp_i, n = round(area_i/dis^2), type = "regular", offset=c(0.5, 0.5))
  ngrid_i <- length(grid_i)
  grid_xyz <- do.call(cbind, Lll2xyz(lat = grid_i@coords[,2], lon = grid_i@coords[,1]))
  block_i <- rep(i, ngrid_i)
  
  weights <- rep(area_i/ngrid_i, ngrid_i)
  return(list(grid_xyz = grid_xyz, block = block_i, weights = weights, ngrid = ngrid_i))
}

grace_block <- lapply(1:nrow(grace_sp), poly_block, dis = 10)

grid_xyz <- do.call(rbind, lapply(grace_block, "[[", "grid_xyz"))
grid_block <- do.call(c, lapply(grace_block, "[[", "block"))
weights <- do.call(c, lapply(grace_block, "[[", "weights"))

A most tricky setting is the scaling factor! Since the projection matrix A approximate map inputs to outputs on different scales, this can make the computation unstable when the difference in the scales are huge. The scale also determines the scale of the marginal variance and predicted uncertainties.

We choose to scale the prediction by \(1e4\) for a \(100km \times 100km\) area resolution. This approximately correspond to a 1 degree resolution and also at the similar scale of the GRACE data polygon size.

A_GRACE_data <- inla.spde.make.A(mesh = mesh0, loc = grid_xyz, block = grid_block,  block.rescale = "count")

Similarly, we create the projection matrix for prediction.

## Same for prediction on a 1 degree resolution grid, we need to know the area of the grid for integration
## generate the prediction grid
gx <- seq(0, 359, 1)
gy <- seq(-89.5, 89.5)
grid_ll <- expand.grid(gx, gy)
pred_data <- data.frame(lon = grid_ll[,1], lat = grid_ll[,2])
coordinates(pred_data) <-c("lon", "lat")
gridded(pred_data) <- TRUE
pred_data <- as(pred_data, "SpatialPolygons")
proj4string(pred_data) <- CRS("+proj=longlat")
areas <- geosphere::areaPolygon(pred_data)/(1000^2)
grid_pred <- do.call(cbind,Lll2xyz(lat = grid_ll[,2], lon = grid_ll[,1]))
A_M_pred <- inla.spde.make.A(mesh = mesh0, loc = grid_pred) 

Next, we create the stack for INLA inference. Since the GIA is given as known in mm/year water equivalence from the ice6g solution, we can substract this part from the GRACE observation first.

2.4.1 Remove the GIA signals from the GRACE

First load the GIA data in the water equivalence height unit. We use the same mesh grid as \(mass\) to reprensent GIA.

ice6g <- read.table(paste0(dd, "WP2-SolidEarth/BHMinputs/GIA/GIA_Pel-6-VM5_ewh.txt"), header = T)
ice6g$x_center <- ifelse(ice6g$x_center < 0, ice6g$x_center+360, ice6g$x_center)
ice6g2<- ice6g[order(ice6g$y_center,ice6g$x_center ),] # re-order the data according to the coordinates

Then with the above data link map, we calculate the GIA contribution to the GRACE data.

gia_loc <- do.call(cbind, Lll2xyz(lat = ice6g$y_center, lon = ice6g$x_center))
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"))
grid_area <- geosphere::areaPolygon(Plist)/(1000^2)
meshLL <- Lxyz2ll(list(x=mesh0$loc[,1], y = mesh0$loc[,2], z = mesh0$loc[,3]))
mesh_sp <- SpatialPoints(data.frame(lon = meshLL$lon, lat = meshLL$lat), proj4string = CRS("+proj=longlat")) 
mesh_idx <- over(mesh_sp, Plist)
GIA_v <- ice6g$trend[mesh_idx]
GIA_grace <- A_GRACE_data %*% GIA_v
grace_sp@data$gia <- as.numeric(GIA_grace)
spplot(grace_sp, "gia", main = "GIA aggreated on the GRACE grid (mm/yr ewh)")

## mass contribution to GRACE
mass_grace <- grace_sp@data$mmweq - GIA_grace
grace_sp@data$mass <- as.numeric(mass_grace)
spplot(grace_sp, "mass", main = "Adjusted mass change agrregated on the GRACE grid (mm/yr ewh)")

Now the GRACE observation is corrected for GIA and we can create the stacks for INLA inference.

## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=grace_sp$mass), A = list(A_GRACE_data),
                     effects = list(M= 1:M_spde$n.spde), tag = "est")
st.pred <- inla.stack(data = list(y=NA), A = list(rbind(A_GRACE_data, A_M_pred)),
                      effects = list(M=1:M_spde$n.spde), tag = "pred")
stM <- inla.stack(st.est, st.pred)

2.5 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))
prec_scale <- c(1/grace_sp$std^2, rep(1, nrow(A_GRACE_data) + nrow(A_M_pred)))

## The formular for modelling the SSH mean
formula = y ~ -1 +  f(M, model = M_spde)

## Run INLA
res_inla <- inla(formula, data = inla.stack.data(stM, spde = M_spde), family = "gaussian",
                 scale =prec_scale, control.family = list(hyper = hyper),
                 control.predictor=list(A=inla.stack.A(stM), compute =TRUE))

2.6 Results

2.6.1 Assemble and save results

Now assemble the inla inference and prediction results.

INLA_pred <- res_inla$summary.linear.predictor
## Extract and project predictions
pred_idx <- inla.stack.index(stM, tag = "pred")$data
idx_grace <- pred_idx[1:nrow(A_GRACE_data)]
idx_grid <- pred_idx[-(1:nrow(A_GRACE_data))]

## mass
M_m <- INLA_pred$mean[idx_grid] 
M_u <- INLA_pred$sd[idx_grid]
proj <- inla.mesh.projector(mesh0, projection = "longlat", dims = c(360,180), xlim = c(0,359), ylim = c(-89.5, 89.5))
M_grid <- expand.grid(proj$x, proj$y)
M_pred <- data.frame(lon = M_grid[,1], lat = M_grid[,2],
                       mean = M_m,
                       u = M_u)

res_M <- list(res_inla = res_inla, spde = M_spde, st = stM, 
            mesh = mesh0,  M_pred = M_pred, Adata = A_GRACE_data, Apred = A_M_pred)

grace_m <- INLA_pred$mean[idx_grace]
grace_u <- INLA_pred$sd[idx_grace]
grace_sp@data$pred_mean <- grace_m
grace_sp@data$pred_u <- grace_u

save(res_M, grace_sp, file =paste0(dd, "WP1-BHM/Experiment2a/exp2a_M.RData"))

2.6.2 Plot the posteriors of the hyper parameters

load(paste0(dd, "WP1-BHM/Experiment2a/exp2a_M.RData"))
pars_M <- marginal_par(res = res_M, process = "M", plot = TRUE)

## The posterior modes
print(paste("The estimated correlation lengths are:", pars_M$rho_mode*6371,  sep = "  "))

## [1] "The estimated correlation lengths are:  822.056906990587"

print(paste("The estimated marginal variances are:", pars_M$sigma_mode,sep = "  "))

## [1] "The estimated marginal variances are:  186.460541260929"

2.7 Plot the predictions

We first plot the predicted \(X_{mass}\) on a 1 degree grid and then the aggregated \(X_m\) on the same super mascons grid as the GRACE data.

M_pred <- res_M$M_pred

M_pred$mean2 <- ifelse(abs(M_pred$mean) > 19, sign(M_pred$mean)*20, M_pred$mean)
M_pred$u2 <- ifelse(abs(M_pred$u) > 5, 7, M_pred$u)
## plot the mean 
lattice::levelplot(mean2 ~ lon + lat, data = M_pred, aspect = "iso", at = seq(-20, 20, 2),
                     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")
                     },
                     main = "The predicited mass change  (mm/yr ewh)", xlab = "longitude", ylab = "latitude")

## plot the uncertainty
lattice::levelplot(u2 ~ lon + lat, data = M_pred, aspect = "iso", col.regions = topo.colors(20),
                     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")
                     },
                     main = "The predicted uncertainties  (mm/yr ewh)", xlab = "longitude", ylab = "latitude")

grace_sp$pred_mean2 <- ifelse(abs(grace_sp@data$pred_mean) > 20, sign(grace_sp@data$pred_mean)*20, grace_sp@data$pred_mean )
spplot(grace_sp, "pred_mean2", at = seq(-21, 21, 2), main = "The predicted mass (mm/yr ewh)")

spplot(grace_sp, "pred_u", main = "The predicted uncertainties (mm/yr ewh)")

3 Sanity checks on the updated mass

To check whether the updated mass is reasonable, we calculate the integrated mass trend and compare it to existing work.

First we check the global trend of the mass change in mm water equivalence height.

order_id <- order(ice6g$y_center,ice6g$x_center ) # re-order the data according to the coordinates
grid_area2 <- grid_area[order_id]

print(paste0("global average is ", mean(M_pred$mean)))

## [1] "global average is 0.0609729539134967"

print(paste0("weighted global average is ", sum(M_pred$mean*areas/sum(areas))))

## [1] "weighted global average is 0.854202064748128"

Then we check the ocean mass change. We need to first identify which grid belongs to the oceans. The values are similar to and between existing studies.

Ocean <- readOGR(dsn = paste0(dd, "WP1-BHM/maps/ne_110m_ocean"), layer = "ne_110m_ocean")

## 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

coords <- ice6g2[,2:3]
coords[,1] <- ifelse(coords[,1] > 180, coords[,1]-360, coords[,1])
grid_inO <- unlist(over(Ocean, SpatialPoints(coords = coords, proj4string = CRS(proj4string(Ocean))), returnList = TRUE))
print(paste0("average over oceans is ", mean(M_pred$mean[grid_inO])))

## [1] "average over oceans is 1.98792302998572"

print(paste0("weighted average over the oceans is ", sum(M_pred$mean[grid_inO]*areas[grid_inO]/sum(areas))))

## [1] "weighted average over the oceans is 1.5071185557515"

print(paste0("average on land is ", mean(M_pred$mean[-grid_inO])))

## [1] "average on land is -3.80983733262189"

print(paste0("weighted average on land is ", sum(M_pred$mean[-grid_inO]*areas[-grid_inO]/sum(areas))))

## [1] "weighted average on land is -0.65291649100337"

4 Update Steric