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

This study references numerous geospatial layers, including polygons, points, and rasters. These layers are preprocessed here and saved in the data folder for reproducibility and code efficiency. Below, vector-format layers (points, etc.) are saved as GeoJSON files and raster-format layers are saved as GeoTiffs.

library(raster)
library(sf)
library(tidyverse)
library(janitor)
library(ggplot2)
library(dplyr)

Setup study bounds

We setup a 100 km x 100 km grid for this study, with a resolution of 2.5 km x 2.5 km (see the paper).

study_bound_coords <- rbind(
  c(460000, 5980000), c(460000, 6080000), c(560000, 6080000), 
  c(560000, 5980000), c(460000, 5980000)
  )

# make a polygon
study_bounds <- list(study_bound_coords) %>%
  st_polygon() %>%
  st_sfc(., crs = 26911) %>%
  st_sf(.)

study_proj = st_crs(study_bounds)$proj4string

st_write(study_bounds, '../data/study_bounds.GeoJSON', delete_dsn = TRUE)

Make a 2500x2500 m raster

We now transform the bounds into a grid using the raster package.

study_raster <- raster(
  extent(study_bounds), 
  resolution = 2500,
  crs = st_crs(study_bounds)$proj4string
  )
study_raster[is.na(study_raster[])] <- 0 

writeRaster(study_raster, '../data/study_raster.tif', overwrite=TRUE)

Creat a Duvernay-wide raster for mapping purposes

Load surface locations

We load the surface hole locations from the AGS well list layer, filter to make sure they are in the study bounds, and combine them with the AGS well list to get both the license, wellname, uwi etc.

Original data can be found here: http://www1.aer.ca/ProductCatalogue/10.html

A copy is saved to data-raw which isn't uploaded to the repository. Generally, anything in data-raw is omitted due to licensing restrictions or politeness of not reproducing an original system of record.

surface_locs = sf::st_read(
  dsn = '../data-raw/sh_shapefile/', layer='ST37_SH_GCS_NAD83') %>% 
  st_transform(., crs = 26911)

study_locs = surface_locs %>% filter(
  st_within(surface_locs, study_bounds,sparse = FALSE)
  ) %>% janitor::clean_names()

well_list = read_tsv(
  '../data-raw/sh_shapefile/WellList.txt', 
  col_names = c('uwi_disp','uwi','flag','name','field','pool','os_area',
                'os_dep','licence', 'status', 'issue_data', 'code','agent', 
                'operator','drill_data', 'td_m','stat_code', 'stat_date', 
                'fluid', 'mode', 'struct', 'scheme', 'scheme_sub')
  )

well_list$uwi_clean = well_list$uwi_disp %>% 
  str_c(1,.) %>%
  str_replace(., "/", "") %>% 
  str_replace(., "/", "0") %>% 
  str_replace_all(., "-", "")

# Merge with surface locations
study_locs = study_locs %>% 
  merge(well_list %>% dplyr::select(licence, uwi = uwi_clean), 
        on='licence', how='left')

st_write(study_locs, '../data/study_surf_locs.GeoJSON', delete_dsn = TRUE)

Load the Duvernay isopach and write the study mask for predictions & figures: http://www.ags.gov.ab.ca/publications/DIG/ZIP/DIG_2008_0124.zip

dv_bounds <- st_read(dsn = '../data-raw/dvisopach/', layer = 'fg1217_py_ll') %>%
  st_transform(., crs = 26911) %>%
  janitor::clean_names() %>%
  dplyr::filter(descr1 != 'Dolomite') %>%
  dplyr::filter(descr1 != 'Mixed dolomite and evaporite') %>%
  dplyr::filter(descr1 != 'Limestone')

st_write(dv_bounds, '../data/entire_dv_bounds.GeoJSON', delete_dsn=TRUE)

study_dv_bounds <- dv_bounds %>%
  st_intersection(study_bounds)

st_write(study_dv_bounds, '../data/study_dv_bounds.GeoJSON', delete_dsn=TRUE)

Load the ground surface from Alberta Geological Survey's 3D Model https://ags.aer.ca/publication/3d-pgf-model-v2

ground_surface <- raster(
  '../data-raw/ags_geomodel/Continuous_Model_Horizons/cmh_01_QNgPg_sediment_top.asc'
  ) %>%
  raster::projectRaster(., crs = study_proj) %>%
  crop(., study_bounds)
names(ground_surface) <- 'ground_elev_masl'

# Verify projections visually
surface_df = ground_surface %>%
  rasterToPoints(., spatial=TRUE) %>%
  data.frame(.)

ggplot() +
  geom_tile(data=surface_df, aes(x=x,y=y,fill=ground_elev_masl)) +
  geom_sf(data=study_dv_bounds, fill=NA, colour='yellow')

writeRaster(ground_surface, '../data/ground_elev_masl.tif', overwrite=TRUE)

Load the Duvernay top from Alberta Geological Survey's 3D Model

dv_top <- raster(
  '../data-raw/ags_geomodel/Discrete_Model_Horizons_TOP/dmh_54_D_Dv_Mu_top.asc'
  ) %>%
  raster::projectRaster(., crs = study_proj) %>%
  crop(., study_bounds)
names(dv_top) <- 'dv_top_masl'

# Verify projections visually
dv_top_df = dv_top %>%
  rasterToPoints(., spatial=TRUE) %>%
  data.frame(.)

ggplot() +
  geom_tile(data=dv_top_df, aes(x=x,y=y,fill=dv_top_masl)) +
  geom_sf(data=study_dv_bounds, fill=NA, colour='yellow')

writeRaster(dv_top, '../data/dv_top_masl.tif', overwrite=TRUE)

Load the Duvernay Base

dv_base <- raster(
  '../data-raw/ags_geomodel/Discrete_Model_Horizons_BASE/dmh_54_D_Dv_Mu_base.asc'
  ) %>%
  raster::projectRaster(., crs = study_proj) %>%
  crop(., study_bounds)
names(dv_base) <- 'dv_base_masl'

# Verify projections visually
dv_base_df = dv_base %>%
  rasterToPoints(., spatial=TRUE) %>%
  data.frame(.)

ggplot() +
  geom_tile(data=dv_base_df, aes(x=x,y=y,fill=dv_base_masl)) +
  geom_sf(data=study_dv_bounds, fill=NA, colour='yellow')

writeRaster(dv_base, '../data/dv_base_masl.tif', overwrite=TRUE)

Get DV Isopach

dv_isopach = dv_top - dv_base
names(dv_isopach) <- 'dv_isopach_m'

dv_isopach_df = dv_isopach %>%
  rasterToPoints(., spatial=TRUE) %>%
  data.frame(.)

ggplot() +
  geom_tile(data=dv_isopach_df, aes(x=x,y=y,fill=dv_isopach_m)) +
  geom_sf(data=study_dv_bounds, fill=NA, colour='yellow')

writeRaster(dv_isopach, '../data/dv_iso_m.tif', overwrite=TRUE)

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