Define Pseudo Observation Polygons

Find suitable area for defining pseudo observation polygons

In this document, we use the 14 GIA forward model solutions to find regions where the observations are certain and close to zero. The 13 GIA forward model solutions are given at 1 degree resolution with the mean trend and estimated errors.

In the following we explore the data sets individually as well as in an ensemble approach.

0 Load useful packages and the data

### Load the 14 GIA data 
data_path <- "Z:/WP2-SolidEarth/GIAforwardModels/textfiles/"
file_names <- system(paste("ls", data_path), intern = TRUE)[1:14]
## Remove the no.11 W&O-EGOD
file_names <- file_names[-11]
GIA_name <- unlist(strsplit(file_names, split = ".txt"))
GIA_priors <- list()
for (i in 1:13){
GIA_priors[[i]] <- read.table(paste0(data_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)))/13)

1 Criteria for selecting zero-value regions

To select the reasonable regions where the GIA are zero with high certainty, we first examine the distributions of the GIA values from the 13 forward models and their ensemble mean.

## All 13 datasets together
alltrend <- sapply(GIA_priors, "[[", "trend")
range(alltrend)

## [1] -8.389504 25.188101

boxplot.stats(alltrend)$conf

## [1] -0.1514417 -0.1492397

## the esemble mean
range(GIA_ensemble$trend)

## [1] -5.718985 13.684467

boxplot.stats(GIA_ensemble$trend)$conf

## [1] -0.1440926 -0.1371992

Plot the original 13 dataset and the ensemble means.

## combine the 13 lists to one dataframe
alltrend <- data.frame(alltrend)
names(alltrend) <- paste0("trend", 1:13)
alltrend$trend_ens <- GIA_ensemble$trend
alltrend$x_center <- GIA_priors[[1]]$x_center
alltrend$y_center <- GIA_priors[[1]]$y_center

library(sp)
raw_trend <- alltrend
coordinates(raw_trend) <- c("x_center", "y_center")
gridded(raw_trend) <- TRUE

brks <- seq(-8.5, 25.5, 2)
colpal <- colorRamps::matlab.like(length(brks) - 1)
spplot(raw_trend,  colorkey = list(at = brks), col.regions=colpal)

Rescale the colorbar to show more details for values within [-2, 2].

plot_rescale <- function(breaks, abs = FALSE, title = NULL){
  n.brks <- length(breaks)
  if(abs){
    trends<- lapply(abs(alltrend[,1:14]), function(x) cut(x, breaks))
  }else{
     trends<- lapply(alltrend[,1:14], function(x) cut(x, breaks))
  }
  new_trend <- do.call(data.frame, trends)
  new_trend$x_center <-alltrend$x_center
  new_trend$y_center <- alltrend$y_center
  
  colpal <- colorRamps::matlab.like(n.brks-1)
  coordinates(new_trend) <- c("x_center", "y_center")
  gridded(new_trend) <- TRUE
  spplot(new_trend,  col.regions=colpal, main = list(title))
  
}


breaks <- c(-8.5, -2, -1, -0.5, -0.1, 0, 0.1, 0.5, 1, 2, 25.5)
plot_rescale(breaks)

Remove signs and rescale the colorbar down to details within 0.5 in absolute value.

breaks <- c( 0, 0.1, 0.2, 0.5, 1, 25.5)
plot_rescale(breaks, abs = TRUE)

Down to the value of 0.5, most of the plots are very similar to each other except for trend3. They include most of the zero value regions. For a final comparison, we plot the zero value regions with threshold 0.1, 0.2, 0.3 and 0.5.

## 0.1
breaks <- c(0, 0.1, 25.5)
plot_rescale(breaks, abs = TRUE, title = "threshold = 0.1")

## 0.2
breaks <- c(0, 0.2, 25.5)
plot_rescale(breaks, abs = TRUE, title = "threshold = 0.2")

## 0.3
breaks <- c(0, 0.3, 25.5)
plot_rescale(breaks, abs = TRUE, title = "threshold = 0.3")

## 0.5
breaks <- c(0, 0.5, 25.5)
plot_rescale(breaks, abs = TRUE, title = "threshold = 0.5")

The map using ensemble mean and thresholds at 0.2 and 0.3 looks reasonable and we use these two in the following processing.

2 Create the zero value polygons.

Now we convert the region defined by the above criteria from pixels to polygons. We will try three threshold values: 0.1, 0.2, 0.3. I would recommend to use 0.2 and 0.3 since they cover more regions and unwanted points can always be removed later.

In the following, we fisrt do the porcessing for threshold value 0.1 and later we do the same for 0.2 and 0.3.

## from Spatialpixels to raster to polygons
library(raster)
GIA_sp <- GIA_ensemble[, c(2,3,4,5)]
GIA_sp$trend01 <- abs(GIA_sp$trend) <0.1
GIA_sp$trend02 <- abs(GIA_sp$trend) <0.2
GIA_sp$trend03 <- abs(GIA_sp$trend) <0.3
coordinates(GIA_sp) <- c("x_center", "y_center")
gridded(GIA_sp) <- TRUE


GIA_ras01 <- raster(GIA_sp, layer = 3) # 3 is the column no of trend01
zeroPoly01 <- rasterToPolygons(GIA_ras01, fun=function(trend01) {trend01 == TRUE}, dissolve=TRUE )

## Loading required namespace: rgeos

plot(zeroPoly01, col = "blue", main = "zero regions -- threshold = 0.1")

Note that there are holes in the large polygons and check if they are correctly labeled.

sapply(zeroPoly01@polygons[[1]]@Polygons, function(x) slot(x, "hole"))

##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [67] FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE
##  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [166] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [188] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE
## [199]  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE
## [210]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
## [221] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

We can convert the spatialPolygonDataFrame object to a data frame that can be plotted by ggplot2.

library(broom)
library(ggplot2)

GIA_zeroDF01 <- tidy(zeroPoly01)

## Regions defined for each Polygons

zeropoly_map01 <- ggplot()+ coord_fixed()+ geom_polygon(data = GIA_zeroDF01, aes(x=long, y=lat, group = group, fill = !hole),colour="black") + ggtitle("Threshold = 0.1")

print(zeropoly_map01)

3 Generate Pseudo observations within the zero polygons

First, generate vertices that evenly distributed over the entire globe. The distance between the vertices is about the same as the correlation length of the of the residual process, which is about 400-500km apart. This is about spreading 3000 points using the Fibonacci spiral on a sphere.

## xyz -- Generate equally distributed points around a sphere -- Fibonacci spiral
fiboSphere1 <- function(n = 1000L) {
  
  phi <- (sqrt(5) + 1) / 2  # golden ratio
  i <- seq(-(n-1), (n-1),2)
  theta <- 2 * pi * i / phi
  sphi <- i/n
  cphi <- sqrt((n+i) * (n-i))/n
  
  x <- cphi * sin(theta)
  y <- cphi * cos(theta)
  z <- sphi
  return(cbind(x,y,z))
}

##lonlat -- Generate equally distributed points around a sphere -- Fibonacci spiral
fiboSphere2 <- function(n = 1000L) {
  phi <- (sqrt(5) + 1) / 2 - 1 # golden ratio
  ga <- phi * 2 * pi           # golden angle
  
  i <- 1L:n
  lon <- ga * i / (2 * pi)
  lon <- 2 * pi * (lon - floor(lon))
  lon <- ifelse(lon <= pi, lon, lon - 2 * pi)/pi*180
  lat <- asin(-1 + 2 * i / n)/pi*180
  
  cbind(lon = lon, lat = lat)
}

fibopoints <- fiboSphere2(n = 3000) 
fibopoints <- data.frame(fibopoints)
fibopoints$lon <- ifelse(fibopoints$lon < 0, fibopoints$lon+359, fibopoints$lon)


## Plot these points 
world_map <- map_data("world2")
map_fibo <- ggplot() + coord_fixed() + 
  geom_polygon(data=world_map, aes(x=long, y=lat, group=group), colour="black", fill = NA)  + 
  geom_point(data=data.frame(fibopoints), aes(x=lon, y=lat), pch=19, col = "blue", alpha=0.5)  + 
  ggtitle("Points evenly distributed on the sphere")

map_fibo

Then remove points that are not in the zero region polygons.

fibsp <- SpatialPoints(fibopoints, proj4string = CRS("+proj=longlat"))
zeroPolys01 <- SpatialPolygons(zeroPoly01@polygons, proj4string = CRS("+proj=longlat"))
fibin01 <- over(fibsp, zeroPolys01)
fibin01 <- ifelse(is.na(fibin01), FALSE, TRUE)
## plot these points
fibpoint_in01 <- data.frame(fibopoints[fibin01, ])

map_fiboZero01 <- zeropoly_map01 + geom_point(data=fibpoint_in01, aes(x=lon, y=lat), pch=19, 
                                    col = "blue", alpha=0.5)  + ggtitle("fibpoints in the zero polygons -- threshold = 0.1")
map_fiboZero01

4 Exclude the GPS area

Given the zero value region, we also want to exclude the region that include real GPS readings. Now we plot the latest GPS data for inspection.

GPS_data <- read.table("Z:/WP2-SolidEarth/GPS/NGL/BHMinputFiles/GPS_v03d.txt",  header = T)
GPS_data$lon <- ifelse(GPS_data$lon < 0, GPS_data$lon + 359, GPS_data$lon)

world_map <- map_data("world2")
map_GPS <-  ggplot() + coord_fixed() + geom_polygon(data=world_map, aes(x=long, y=lat, group=group), colour="black", fill = NA) + geom_point(data=GPS_data, aes(x=lon, y=lat), pch=19, 
                                    col = "red", fill = "red", alpha=0.7) + ggtitle("Real GPS observations")
map_GPS

We remove points that are too close to the GPS observations (say within the correlation length 400-500km , that is about 5 degree in longlat).

fibin01gps <- apply(fibpoint_in01, 1, function(x) all(apply(GPS_data[,2:3], 1, function(y) dist(rbind(y, x)) > 5)))
fibpoint_in01gps <- fibpoint_in01[fibin01gps,]

map_fibo01 <- zeropoly_map01 + geom_point(data=fibpoint_in01gps, aes(x=lon, y=lat), pch=19, 
                                    col = "blue", alpha=0.5)  +
  geom_point(data=GPS_data, aes(x=lon, y=lat), pch=19, 
                                    col = "red", fill = "red", alpha=0.7) +
  ggtitle("fibpoints in the zero polygons and away from GPS -- threshold = 0.1")
map_fibo01

Finally, assemble the pseudo observation and real observation into the same data set. The pseudo observations are all zero and the related uncertainty is set to be tenth of the threshold values.

nps01 <- nrow(fibpoint_in01gps)
pseudo_data01 <- data.frame(ID = rep("pseudo", nps01), lon = fibpoint_in01gps$lon, lat = fibpoint_in01gps$lat, 
                          trend = 0, std = 0.01)
obs_all01 <- rbind(GPS_data, pseudo_data01)
obs_all01$source <- ifelse(obs_all01$ID == "pseudo", "pseudo", "real")

map_final01 <- ggplot() + coord_fixed() + geom_polygon(data=world_map, aes(x=long, y=lat, group=group), colour="black", fill = NA) + 
  geom_polygon(data = GIA_zeroDF01, aes(x=long, y=lat, group = group), colour="red", fill = NA) +
  geom_point(data=obs_all01, aes(x=lon, y=lat, colour = source), pch=19) + 
  ggtitle("Final map of real and pseudo observations -- threhold = 0.1")

map_final01

Plot with another projection.

library(mapproj)
map_final01 + coord_map(projection = "orthographic")

5 Threshold = 0.2

Now we do the same using threshold value 0.2. First create the admissible polygons.

## Generate the polygons
GIA_ras02 <- raster(GIA_sp, layer = 4) # 4 is the column no of trend02
zeroPoly02 <- rasterToPolygons(GIA_ras02, fun=function(trend02) {trend02 == TRUE}, dissolve=TRUE ) 
plot(zeroPoly02, col = "blue", main = "zero regions -- threshold = 0.2")

## Checking holes
sapply(zeroPoly02@polygons[[1]]@Polygons, function(x) slot(x, "hole"))

##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [34] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [166] FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [177] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [188] FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [199]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [210]  TRUE  TRUE  TRUE  TRUE  TRUE

## plot in ggplot2
GIA_zeroDF02 <- tidy(zeroPoly02)

## Regions defined for each Polygons

zeropoly_map02 <- ggplot()+ coord_fixed()+ geom_polygon(data = GIA_zeroDF02, aes(x=long, y=lat, group = group, fill = !hole),colour="black") + ggtitle("Threshold = 0.2")
print(zeropoly_map02)

Then remove the points outside theses polygons.

fibsp <- SpatialPoints(fibopoints, proj4string = CRS("+proj=longlat"))
zeroPolys02 <- SpatialPolygons(zeroPoly02@polygons, proj4string = CRS("+proj=longlat"))
fibin02 <- over(fibsp, zeroPolys02)
fibin02 <- ifelse(is.na(fibin02), FALSE, TRUE)
## plot these points
fibpoint_in02 <- data.frame(fibopoints[fibin02, ])

map_fiboZero02 <- zeropoly_map02 + geom_point(data=fibpoint_in02, aes(x=lon, y=lat), pch=19, 
                                    col = "blue", alpha=0.5)  + ggtitle("fibpoints in the zero polygons -- threshold = 0.2")
map_fiboZero02

Then remove points that are too close the real GPS observations.

fibin02gps <- apply(fibpoint_in02, 1, function(x) all(apply(GPS_data[,2:3], 1, function(y) dist(rbind(y, x)) > 5)))
fibpoint_in02gps <- fibpoint_in02[fibin02gps,]

map_fibo02 <- zeropoly_map02 + geom_point(data=fibpoint_in02gps, aes(x=lon, y=lat), pch=19, 
                                    col = "blue", alpha=0.5)  +
  geom_point(data=GPS_data, aes(x=lon, y=lat), pch=19, 
                                    col = "red", fill = "red", alpha=0.7) +
  ggtitle("fibpoints in the zero polygons and away from GPS -- threshold = 0.2")
map_fibo02

Finally, assemble the data and plot.

nps02 <- nrow(fibpoint_in02gps)
pseudo_data02 <- data.frame(ID = rep("pseudo", nps02), lon = fibpoint_in02gps$lon, lat = fibpoint_in02gps$lat, 
                          trend = 0, std = 0.2)
obs_all02 <- rbind(GPS_data, pseudo_data02)
obs_all02$source <- ifelse(obs_all02$ID == "pseudo", "pseudo", "real")

write.table(obs_all02, file = "Z:/WP1-BHM/Experiment1b/GIA_RGL/pseudo02.txt", 
            row.names = FALSE, eol = "\r\n")

map_final02 <- ggplot() + coord_fixed() + geom_polygon(data=world_map, aes(x=long, y=lat, group=group), colour="black", fill = NA) + 
  geom_polygon(data = GIA_zeroDF02, aes(x=long, y=lat, group = group), colour="red", fill = NA) +
  geom_point(data=obs_all02, aes(x=lon, y=lat, colour = source), pch=19) + 
  ggtitle("Final map of real and pseudo observations -- threhold = 0.2")

map_final02

map_final02 + coord_map(projection = "orthographic")

6 Threshold = 0.3

Now we do the same using threshold value 0.3. First create the admissible polygons.

## Generate the polygons
GIA_ras03 <- raster(GIA_sp, layer = 5) # 5 is the column no of trend03
zeroPoly03 <- rasterToPolygons(GIA_ras03, fun=function(trend03) {trend03 == TRUE}, dissolve=TRUE ) 
plot(zeroPoly03, col = "blue", main = "zero regions -- threshold = 0.3")

## Checking holes
sapply(zeroPoly03@polygons[[1]]@Polygons, function(x) slot(x, "hole"))

##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE
##  [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [111] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
## [144]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [155]  TRUE  TRUE

## plot in ggplot2
GIA_zeroDF03 <- tidy(zeroPoly03)

## Regions defined for each Polygons

zeropoly_map03 <- ggplot()+ coord_fixed()+ geom_polygon(data = GIA_zeroDF03, aes(x=long, y=lat, group = group, fill = !hole),colour="black") + ggtitle("Threshold = 0.3")
print(zeropoly_map03)