0 Introduction

In this document, we apply the pseudo polygon method on GIA. The idea is that GIA is a stationary porcess on a subset of \(S^2\). The subset is given by removing polygons where the values are certainly zero according to experts from the entire sphere.

1 Generate the polygons and mesh

First we generate the mesh for a discreate reparesentation of the GIA process. The mesh is restricted on a subset of the sphere defined by the pseudo-polygons. More details about choosing and generating the polygons can be found here. The following chunk generate the pseudo polygons using the ensemble mean of 13 GIA mode solutions for a given threshold value, say \(0.3\).

### Load the GIA forward solutions
if(Sys.info()["sysname"] == "Windows"){
  GIA_path <- "Z:/WP2-SolidEarth/BHMinputs/GIA/" 
}else if(Sys.info()["sysname"] == "Ubuntu"){
GIA_path <- "~/globalmass/WP2-SolidEarth/GIAforwardModels/textfiles/" #(on ubuntu desktop )
}else{
GIA_path <- "/./projects/GlobalMass/WP2-SolidEarth/GIAforwardModels/textfiles/" #(on server)
}

file_names <- system(paste("ls", GIA_path), intern = TRUE)
## Remove unwanted files: no.1-4 ice6g old, n0.6 isce6g NA,  no.8 s&s1, no.12 SVvLALT, no.13 W&O-EGOD and no.17 readme files
file_names <- file_names[-c(1:4,6, 8,  12, 13, 17)]
GIA_name <- unlist(strsplit(file_names, split = ".txt"))
GIA_priors <- list()
for (i in 1:8){
GIA_priors[[i]] <- read.table(paste0(GIA_path, file_names[i]), header = T)
}

## create the ensemble data sets
GIA_ensemble <- GIA_priors[[1]]
GIA_ensemble$trend <- rowMeans(sapply(GIA_priors, "[[", "trend"))
GIA_ensemble$std <- sqrt((rowSums(sapply(GIA_priors, function(x) (x$std)^2)))/8)

## from Spatialpixels to raster to polygons
library(raster)
## Loading required package: sp
GIA_sp <- GIA_ensemble[, c(2,3,4,5)]
GIA_sp$inPoly1 <- abs(GIA_sp$trend) < 0.3 & abs(GIA_sp$std) < 0.4
GIA_sp$inPoly2 <- abs(GIA_sp$trend) < 0.3 
coordinates(GIA_sp) <- c("x_center", "y_center")
gridded(GIA_sp) <- TRUE
spplot(GIA_sp, c("inPoly1", "inPoly2"))

GIA_ras <- raster(GIA_sp, layer = 3) # 3 is the column no of inPoly
zeroPoly <- rasterToPolygons(GIA_ras, fun=function(inPoly) {inPoly == TRUE}, dissolve=TRUE )
## Loading required namespace: rgeos
plot(zeroPoly, col = "blue", main = "zero regions -- threshold = 0.3")

## Remove polygons that are too small
zeroPolys <- zeroPoly@polygons[[1]]@Polygons
polyareas <- sapply(zeroPolys, function(x) x@area)
polyholes <- sapply(zeroPolys, function(x) x@hole)

zeropolys2 <- zeroPolys[polyareas > 200 ] 
polyareas2 <- sapply(zeropolys2, function(x) x@area)
polyholes2 <- sapply(zeropolys2, function(x) x@hole)

zeroPoly2 <- zeroPoly
zeroPoly2@polygons[[1]]@Polygons <- zeropolys2
  
plot(zeroPoly2, col = "blue",  main = "zero regions -- threshold = 0.3 (refined) ")

Next we generate the mesh and separate the triangles in/out-side of the polygons. We spread fibonacci points on the sphere as initial locations of the triangle vertices. The number of the points controls the resolution. The goal is to produce a map with approximately one degree resolution, hence we also need a similar resolution for the triangulation. Not 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 Fibonacci points on the sphere
fiboSphere2 <- function(N = 1000L, L0 = FALSE) {
  ## Reference (note that points generated from 2D are slightly different from 3D)
  ## Measurement of Areas on a Sphere Using Fibonacci and Latitude–Longitude Lattices (2010)
  phi <- (sqrt(5) + 1) / 2 - 1 # golden ratio
  ga <- phi * 2 * pi           # golden angle

  i <- seq(-N, N)
  P <- 2 * N + 1
  lat <- asin(2*i / P) * 180 / pi
  if(L0){
  lon <- ((2 * pi * i / phi) %% pi) * 360 / pi
  }else{
    lon <- ((2 * pi * i / phi) %% pi) * 360 / pi - 180
    }
  cbind(lon = lon, lat = lat)
}


fibo_points <- fiboSphere2(N = 12960)
#fibo_points <- fiboSphere2(N = 1000)
library(GEOmap);library(INLA)
## Loading required package: Matrix
## This is INLA_17.06.20 built 2017-06-20 03:42:30 UTC.
## See www.r-inla.org/contact-us for how to get help.
if(Sys.info()['sysname'] == "Linux"){INLA:::inla.dynload.workaround()}
mesh_points_xyz <- do.call(cbind, Lll2xyz(lat = fibo_points[,2], lon = fibo_points[,1]))
mesh <- inla.mesh.2d(loc = mesh_points_xyz, cutoff = 0.01, max.edge = 0.5)
summary(mesh) # give the desired number of vertices and triangles.
## 
## Manifold:    S2
## CRS: N/A
## Vertices:    25921
## Triangles:   51838
## Boundary segm.:  0
## Interior segm.:  0
## xlim:    -1 0.9999485
## ylim:    -0.9999871 0.9999871
## zlim:    -0.9999614 0.9999614

Now separate the triangles by the pseudo-polygons.

source("C:/ZSwork/glbm/Experiment1b/mixture_model/BarrierModel/functions-barriers-dt-models-march2017.R")
## rgeos version: 0.3-25, (SVN revision 555)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## 
## Attaching package: 'grid'
## The following object is masked from 'package:GEOmap':
## 
##     explode
## Spam version 2.1-1 (2017-07-02) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## The following object is masked from 'package:INLA':
## 
##     Oral
## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve
## Loading required package: maps
## rgdal: version: 1.2-11, (SVN revision 676)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/R_mylibs/R/win-library/3.4/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/R_mylibs/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5
mesh <- dt.mesh.addon.posTri(mesh = mesh, globe = TRUE)
# - Add on mesh$posTri
# - - contains the positions of the triangles

## checking which mesh triangles are inside the land
## First convert xyz to lonlat
Tlonlat <- Lxyz2ll(list(x = mesh$posTri[,1], y = mesh$posTri[,2], z = mesh$posTri[,3]))
Tlonlat$lon <- ifelse(Tlonlat$lon >=0, Tlonlat$lon, Tlonlat$lon + 359)
mesh$Trill <- cbind(lon = Tlonlat$lon, lat =Tlonlat$lat)
TinPoly <- unlist(over(zeroPoly2, SpatialPoints(coords=mesh$Trill), returnList=T))
TAll <- 1:mesh$t
ToutPoly <- TAll[-TinPoly]
Omega = dt.Omega(list(TinPoly, 1:mesh$t), mesh)

## Plot the result in 3d
plot(mesh, t.sub = Omega[[2]], col = "lightblue", rgl = TRUE )
plot(mesh, t.sub = Omega[[1]], col = "yellow",  rgl = TRUE, add = TRUE)

Need to build an mesh from the subset of triangles that can be used in INLA independently. Find the subset of the vertices, and change all indices, etc.

## get the subset of the mesh for buiding the Q
mesh.sub <- function(mesh, Omega, i = 2){
  mesh_sub <- mesh
  id_tri <- Omega[[i]] # triangel id in the subset
  id_vet <- unique(as.vector(mesh$graph$tv[Omega[[i]], ])) # vertex id in the subset
  
  ## Now update the index for all
  id_all_v <- mesh$idx$loc
  id_all_v[-id_vet] <- NA
  id_all_v[which(!is.na(id_all_v))] <- 1:length(id_vet)
  
  id_all_t <- 1:mesh$t
  id_all_t[-id_tri] <- NA
  id_all_t[which(!is.na(id_all_t))] <- 1:length(id_tri)
  ## Modify the mesh to be the subset
  
  
  ## get the subset of the graph
  sub_tv <- mesh$graph$tv[id_tri, ]
  sub_vt <- mesh$graph$vt[id_vet,]
  sub_tt <- mesh$graph$tt[id_tri,]
  sub_tti <- mesh$graph$tti[id_tri,]
  sub_vv <- mesh$graph$vv[id_vet,id_vet]
  
  ## Change the triangle indices of the subset
  sub_tv2 <- cbind(id_all_v[sub_tv[,1]], id_all_v[sub_tv[,2]], id_all_v[sub_tv[,3]])
  sub_vt2 <- id_all_t[sub_vt]
  sub_tt2 <- cbind(id_all_t[sub_tt[,1]], id_all_t[sub_tt[,2]], id_all_t[sub_tt[,3]])
  
  ## Now assemble 
  mesh_sub$t <- length(id_tri) 
  mesh_sub$n <- length(id_vet)
  mesh_sub$loc <- mesh_sub$loc[sort(id_vet), ]
  mesh_sub$graph$tv <- sub_tv2
  mesh_sub$graph$vt <- sub_vt2
  mesh_sub$graph$tt <- sub_tt2
  mesh_sub$graph$tti <- sub_tti
  mesh_sub$graph$vv <- sub_vv
  
  mesh_sub$idx$loc <- 1:length(id_vet)
  mesh_sub$posTri <- mesh_sub$posTri[Omega[[i]],]
  mesh_sub$Trill <- mesh_sub$Trill[Omega[[i]],]
  return(mesh_sub)
}

mesh_outPoly <- mesh.sub(mesh, Omega, 2)
plot(mesh_outPoly, rgl = T)
plot(mesh_outPoly)

Build an SPDE model from this mesh and simulate some data. Plot the data and compare them with the subset of the complete mesh.

mu_r <- 500/6371
v_r <- (1000/6371)^2
mu_s <- 20
v_s <- 40^2

## Transform the parameters for the SPDE_GMRF approximation
Tlognorm <- function(mu, v){
  logv <- log(1 + v/mu^2)
  logmu <- log(mu^2) - 0.5*log(mu^2 + v)
  return(c(logmu, logv))
}
trho <- Tlognorm(mu_r, v_r)
tsigma <- Tlognorm(mu_s, v_s)


 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
GIA_spde <- inla.spde2.matern(mesh_outPoly, 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)))
idd <- mesh_outPoly$idx$loc
GIA_spde2 <- 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)))
Q <- inla.spde.precision(spde = GIA_spde, theta = c(0,0))
Q2 <- inla.spde.precision(spde = GIA_spde2, theta = c(0,0))[idd, idd]
u1 <- inla.qsample(n=1, Q=Q, seed = 23)[,1]
u2 <- inla.qsample(n=1, Q=Q2, seed = 23)[,1]


local.plot.field = function(field, mesh, ...){
  proj = inla.mesh.projector(mesh,  projection = "longlat", dims = c(360,180), xlim = c(0, 360), ylim = c(-90, 90))
  field.proj = inla.mesh.project(proj, field)
  image.plot(list(x = proj$x, y=proj$y, z = field.proj),
              ...)
}

local.plot.field(u1, mesh = mesh_outPoly, main="The true (simulated) spatial field", asp = 1, zlim = c(-30, 30))  

local.plot.field(u2, mesh = mesh_outPoly, main="The true (simulated) spatial field", asp = 1, zlim = c(-30, 30))

2 Data preparation

2.1 GIA forward model solution

Load the GIA forward solution and project the values to the mesh grid

if(Sys.info()["sysname"] == "Windows"){
  ice6g <- read.table("Z:/WP2-SolidEarth/BHMinputs/GIA/GIA_Pel-6-VM5.txt", header = T)
}else if(Sys.info()["sysname"] == "Ubuntu"){
  ice6g <- read.table("~/globalmass/WP2-SolidEarth/GIAforwardModels/textfiles/GIA_Pel-6-VM5.txt", header = T)
}else{
  ice6g <- read.table("/./projects/GlobalMass/WP2-SolidEarth/GIAforwardModels/textfiles/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"))
## Convert the Cartesian to longlat with long in [0, 360)
MlocLL <- Lxyz2ll(list(x=mesh_outPoly$loc[,1], y = mesh_outPoly$loc[,2], z = mesh_outPoly$loc[,3]))
MlocLL$lon <- ifelse(MlocLL$lon >= -0.5, MlocLL$lon, MlocLL$lon + 360)
M_sp <- SpatialPoints(data.frame(lon = MlocLL$lon, lat = MlocLL$lat), proj4string = CRS("+proj=longlat")) 
Midx <- over(M_sp, Plist)
GIA_prior <- ice6g$trend[Midx]

2.2 GPS data

Load the processed GPS data and remove those inside the pseudo-polygons. Stitch the polygon boundaries with pseudo observations.

if(Sys.info()["sysname"] == "Windows"){
  GPSV4b <- read.table("Z:/WP2-SolidEarth/BHMinputs/GPS/GPS_v04b.txt", header = T)
}else if(Sys.info()["sysname"] == "Ubuntu"){
  GPSV4b <- read.table("~/globalmass/WP2-SolidEarth/GPS/NGL/BHMinputFiles/GPS_v04b.txt", header = T)
}else{
  GPSV4b <- read.table("/./projects/GlobalMass/WP2-SolidEarth/GPS/NGL/BHMinputFiles/GPS_v04b.txt", header = T)
}


GPS_inPoly <- unlist(over(zeroPoly, SpatialPoints(coords = cbind(GPSV4b$lon, GPSV4b$lat)), returnList=T))
GPS_All <- 1:nrow(GPSV4b)
GPS_outPoly <- GPS_All[-GPS_inPoly]
plot(GPSV4b[GPS_outPoly,c("lon", "lat")], pch = "+")

GPS_data <- GPSV4b[GPS_outPoly,]
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"))

Midx2 <- over(GPS_sp, Plist)
GPS_mu <- ice6g$trend[Midx2]
GPS_data$trend0 <- GPS_data$trend - GPS_mu

We also need to add some pseudo observations along the boudaries of the polygons to make smooth transition of the predictions.

## get the boundary of the polygons
boundlines <- as(zeroPoly2, 'SpatialLines') 
obs_bounds <- spsample(boundlines, n = 200, type = "random") # not points more than specified
nobsb <-nrow(obs_bounds@coords) 
obs_df <- data.frame(ID = rep("pseudo", nobsb), lon = obs_bounds@coords[,1], lat = obs_bounds@coords[,2],
                     trend = rep(0, nobsb), std = rep(0.05, nobsb), trend0 = rep(0, nobsb))
obs_xyz <- do.call(cbind, Lll2xyz(lat = obs_bounds@coords[,2], lon = obs_bounds@coords[,1]))

GPS_all <- rbind(GPS_data, obs_df)
GPS_allxyz <- rbind(GPS_loc, obs_xyz)

3 Inference on the subset

## Link the process to observations and predictions
Mesh_GIA <- mesh_outPoly
A_data <- inla.spde.make.A(mesh = Mesh_GIA, loc = GPS_allxyz)

## try predict in zero polygon and see what inla does
zero_lat <- c(-40, 10, 25)
zero_lon <- c(250, 150, 0)
zero_loc <- do.call(cbind, Lll2xyz(lat = zero_lat, lon = zero_lon))

A_pred <- inla.spde.make.A(mesh = Mesh_GIA, loc = rbind(GPS_allxyz, Mesh_GIA$loc, zero_loc))

A_zero <- inla.spde.make.A(mesh = Mesh_GIA, loc = zero_loc)
## Create the estimation and prediction stack
st.est <- inla.stack(data = list(y=GPS_all$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
hyper <- list(prec = list(fixed = TRUE, initial = 0))
formula = y ~ -1 +  f(GIA, model = GIA_spde)
prec_scale <- c(1/GPS_all$std^2, rep(1, nrow(A_pred)))

## Run INLA
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_all)]
GIA_idx <- pred_idx[(nrow(GPS_all)+1):(length(pred_idx) - 3)]
zero_idx <- pred_idx[(length(pred_idx) - 2):length(pred_idx)]

## GPS 
GPS_u <- INLA_pred$sd[GPS_idx]
GPS_pred <- data.frame(lon = GPS_all$lon, lat = GPS_all$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))))
## zero
zero_m <- INLA_pred$mean[zero_idx]
zero_m
## [1] -3.858872e-05 -1.678254e-10 -1.678254e-10
# inla predicts all zero
zero_sd <- INLA_pred$sd[zero_idx]
zero_sd
## [1] 1.6248397442 0.0005530715 0.0005530715
# inla predicts all zero
ress <- list(res_inla = res_inla, spde = GIA_spde, st = stGIA, 
            mesh = Mesh_GIA, GPS_pred = GPS_pred, GIA_pred = GIA_pred)

Plot the posteriors of the hyper parameters.

res_inla <- ress$res_inla
GIA_spde <- ress$spde
pars_GIA <- inla.spde2.result(res_inla, "GIA", GIA_spde, do.transf=TRUE)
theta_mean <- pars_GIA$summary.theta$mean
theta_sd <- pars_GIA$summary.theta$sd

## Find the mode of rho and sigma^2
lrho_mode <- pars_GIA$summary.log.range.nominal$mode
lrho_mean <- pars_GIA$summary.log.range.nominal$mean
lrho_sd <- pars_GIA$summary.log.range.nominal$sd
rho_mode <- exp(lrho_mean - lrho_sd^2)

lsigma_mode <- pars_GIA$summary.log.variance.nominal$mode
lsigma_mean <- pars_GIA$summary.log.variance.nominal$mean
lsigma_sd <- pars_GIA$summary.log.variance.nominal$sd
sigma_mode <- exp(lsigma_mean - lsigma_sd^2)

plot(pars_GIA$marginals.range.nominal[[1]], type = "l",
     main = bquote(bold(rho("mode") == .(round(rho_mode, 4))))) # The posterior from inla output

plot(pars_GIA$marginals.variance.nominal[[1]], type = "l", 
     main = bquote(bold({sigma^2}("mode") == .(round(sigma_mode, 4))))) # The posterior from inla output

## The estimated correlation length is about 568km
rho_mode*6371
## [1] 523.3128

Plot the predictions.

library(ggplot2)
library(grid)
library(gridExtra)

GPS_pred <- ress$GPS_pred
GIA_pred <- ress$GIA_pred

GIA_pred$mean2 <- GIA_pred$mean
GIA_pred$mean2[is.na(GIA_pred$mean2)] <- 0
colpal <- colorRamps::matlab.like(12)

## Plot predicted mean
map_res <- function(data, xname, yname, fillvar, colpal, limits=NULL, title){
  plot_data <- data[c(xname, yname, fillvar)]
  names(plot_data) <- c("X", "Y", "Fill")
  
  beauty <- 
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
          panel.background = element_rect(fill = "white", colour = 'white'), 
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10), 
          axis.text = element_text(size = 10),
          axis.title = element_text(size = 10),
          axis.line = element_line(size = 1),
          plot.title = element_text(hjust = 0.5, size = 15),
          plot.subtitle = element_text(hjust = 0.5, size = 10),
          panel.border = element_blank())
  world_map <- map_data("world2")
  baseworld <- geom_polygon(data = world_map, aes(x=long, y=lat, group=group), colour="grey", fill = NA, alpha = 0.5)  
  
  colbar <- guide_colorbar(barwidth = 2, barheight = 10, label.position = "right", title.position = "bottom")
  
  Map <- ggplot(plot_data) + geom_raster(aes(x = X, y = Y, fill = Fill)) + coord_fixed() + 
    xlab("Longitude") + ylab("Latitude") + 
    scale_x_continuous(limits=c(0,360),  expand = c(0, 0)) + 
    scale_y_continuous(limits=c(-90,90),  expand = c(0, 0)) + 
    scale_fill_gradientn(colors = colpal, name = "mm/yr", limits = limits,  guide = colbar) 
  
  Map <- Map +  baseworld + ggtitle(title) + beauty
  return(Map)
}

map_prior <- map_res(data = ice6g, xname = "x_center", yname = "y_center", fillvar = "trend", 
                     colpal = colpal,  limits = c(-7, 22), title = "Prior GIA mean field")

## Plot the GIA predicted mean
map_GIA <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "mean", 
                   colpal = colpal,  limits = c(-7, 22), title = "Predicted GIA")

map_GIA2 <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "mean2", 
                   colpal = colpal,  limits = c(-7, 22), title = "Predicted GIA patched")

## Plot the GIA difference map
map_diff <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "diff", 
                    colpal = colpal, limits = c(-8, 8), title = "GIA difference: Updated - Prior")

## Plot the GIA difference map
map_sd <- map_res(data = GIA_pred, xname = "lon", yname = "lat", fillvar = "u", 
                    colpal = colpal, limits = c(0, 4), title = "Predicted uncertainties")

 
## Display
print(map_prior)

print(map_GIA)

print(map_GIA2)

print(map_diff)

print(map_sd)