In ScottHMcKean/duvernay_geomech_model: A Geostatistical Model of Geomechanical Properties in the Duvernay Play

Geostatistical Cosimulation

This document goes over development of a LMC and cosimulation of stress data

# geospatial packages
library(sf)
library(gstat)
library(raster)
library(sp)
library(duvernaygeomechmodel)

Load spatial dataframes

study_bounds = st_read('../data/study_dv_bounds.GeoJSON')
study_raster = raster::raster('../data/study_raster.tif')

# make a prediction grid
st_grid <- rasterToPoints(study_raster, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")

Load v, g, and lmc_in values

vgm_in = readRDS('../output/geostatistical_models/stress_vgmobj.rds')
v = vgm_in$v
g = vgm_in$g
lmc_in = vgm_in$lmc_in
norm_objs = vgm_in$norm_objs
data_objs = vgm_in$data_objs

Make predictions

pred_spdf = predict(g, st_grid, nsim=1000)

Back transform variables

sv_invert_pred <- function(x){
  predict(norm_objs$sv, newdata=x, inverse=TRUE)
}

pp_invert_pred <- function(x){
  predict(norm_objs$pp, newdata=x, inverse=TRUE)
}

sh_invert_pred <- function(x){
  predict(norm_objs$sh, newdata=x, inverse=TRUE)
}

# invert from normalized values
trans = pred_spdf
trans@data <- pred_spdf@data %>%
  mutate(across(contains('sv'),sv_invert_pred)) %>%
  mutate(across(contains('pp'),pp_invert_pred)) %>%
  mutate(across(contains('sh'),sh_invert_pred))

saveRDS(trans, '../output/geostatistical_models/stress_pred_trans.rds')

Summary statistics

# get rowwise summary statistics for the values (for mapping)
summary = trans@data %>%
  rowwise() %>%
  summarize(
    sv_mean = mean(c_across(contains('sv'))),
    sv_sd = sd(c_across(contains('sv'))),
    sh_mean = mean(c_across(contains('sh'))),
    sh_sd = sd(c_across(contains('sh'))),
    pp_mean = mean(c_across(contains('pp'))),
    pp_sd = sd(c_across(contains('pp')))
  ) %>%
  mutate(x = pred_spdf@coords[,1]) %>%
  mutate(y = pred_spdf@coords[,2])

write_csv(summary, '../output/geostatistical_models/stress_summary_grid.csv')

Plots

ggplot() +
  geom_tile(data = summary, aes(x=x, y=y, fill = pp_sd)) +
  geom_sf(data = study_bounds, fill=NA, colour='yellow', size=1.5) +
  coord_sf(datum=st_crs(26911)) +
  theme_minimal() +
  ggsave('../output/geostatistical_models/stress_pp_sd.pdf')

ScottHMcKean/duvernay_geomech_model documentation built on Sept. 12, 2022, 10:08 a.m.