In KWB-R/dwc.ar4gw: R Package for Preparing Modflow Output Data for Articifical Reality Visualisation

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
is_ghactions <- identical(Sys.getenv("CI"), "true")

Install

Install the R package dwc.ar4gw by running the following code:

### Optionally: specify GitHub Personal Access Token (GITHUB_PAT)
### See here why this might be important for you:
### https://kwb-r.github.io/kwb.pkgbuild/articles/install.html#set-your-github_pat

# Sys.setenv(GITHUB_PAT = "mysecret_access_token")

# Install package "remotes" from CRAN
if (! require("remotes")) {
  install.packages("remotes", repos = "https://cloud.r-project.org")
}

# Install KWB package 'kwb-r/dwc.ar4gw' from GitHub
remotes::install_github("kwb-r/dwc.ar4gw", dependencies = TRUE)

Define paths

# Load the R package
library(dwc.ar4gw)

### define paths
paths_list <- list(
  tdir = tempdir(), 
  cloud_dir = "dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005",
  local_dir = "C:/kwb/projects/<cloud_dir>",
  download_dir = tempdir(),
  modflow_dir = "<root_data>",
  model_name = "MB16",
  mfout = "<modflow_dir>/<model_name>",
  mf_cbc = "<mfout>.cbc",
  mf_nam = "<mfout>.nam",
  mf_bas = "<mfout>.bas",
  mf_dis = "<mfout>.dis",
  text_heads = "<mfout>.fhd",
  text_drawdown = "<mfout>.fdn",
  binary_cell.by.cell = "<mfout>.cbc",
  binary_heads = "<mfout>.bhd",
  binary_drawdown = "<mfout>.bdn"
)

paths <- kwb.utils::resolve(paths_list, root_data = paths_list$download_dir)
paths_local <- kwb.utils::resolve(paths_list, root_data = paths_list$local_dir)

Modflow Files

This tutorial needs a few MODFLOW files which can be provided locally from a predefined folder on your Windows computer or directly downloaded from the DWC Cloud. Within this tutorial the latter approach is used.

Currently all input files in the MODFLOW .nam file are required and two output files:

<model_name>.nam (input: name file)
<model_name>.dis (input: discretisation file)
<model_name>.lst|bas|oc|gmg|bcf|wel|riv|rch|evt (not loaded, but need to be available in model directory)
<model_name>.cbc (output: cell by cell budget fie)
<model_name>.bhd (output: binary heads file)

where for <model_name> we use r paths_list$model_name, which is defined in the paths_list in the previous chapter Define Paths.

Locally

Copy the content of the DWC cloud folder r paths$cloud_dir to r paths_local$local_dir as shown below:

paths_local

paths <- paths_local

Cloud

You can download the required model files from the DWC cloud if you are a registered user with access to the folder r paths$cloud_dir

For doing so follow the steps below:

Open RStudio and run usethis::edit_r_environ()
In the opened window add the required environment variables

NEXTCLOUD_URL = "https://<replace-with-dwc-cloud-url>"
NEXTCLOUD_USER = "<your-dwc-cloud-username>" # your username
NEXTCLOUD_PASSWORD = "your-nextcloud-app-password" ### see details below

For creating <your-nextcloud-app-password>:

go to: https://replace-with-dwc-cloud-url/index.php/settings/user/security
scroll down to create new app password
select a name e.g. r-script and copy the token and replace your-nextcloud-app-password
Finally you need to restart Rstudio and proceed with the code below:

required_mffiles <- paste0(paste0(c("nam", "dis", "lst", "bas", "oc", "gmg",
         "bcf", "wel", "riv", "rch", "evt", "bhd", "cbc"), 
       collapse = "$|"), "$")


# Download .cbc and .bhd and .dis files
mf_files <- kwb.nextcloud::list_files(
  paths$cloud_dir,
  full_info = TRUE) %>%
  dplyr::filter(stringr::str_detect(.data$file,
                                    pattern = sprintf("^%s\\.",  
                                                      paths$model_name))) %>% 
  dplyr::filter(stringr::str_detect(.data$file,
                                    pattern = required_mffiles))

mf_files

kwb.nextcloud::download_files(href = mf_files$href,
                              target_dir = paths$download_dir)

Flopy release

Install flopy

Install latest flopy release from conda https://anaconda.org/conda-forge/flopy-

### create an environment "flopy" for installing "flopy" 
kwb.python::conda_py_install(env_name = "flopy", 
                             pkgs = list(conda = c("flopy", "python", "pyshp"), 
                                         py = NULL))

Use flopy

#reticulate::use_miniconda("flopy", required = TRUE)
# reticulate::py_help(object = flopy$utils$postprocessing$get_extended_budget)
extended_budget <- dwc.ar4gw::get_extended_budget(cbcfile = paths$binary_cell.by.cell) 
summary(extended_budget$time_1.0$Qx_ext)

str(extended_budget)

flopy <- dwc.ar4gw::import_flopy()

# reticulate::py_help(object = flopy$utils$binaryfile$HeadFile)
heads <- flopy$utils$binaryfile$HeadFile(filename = paths$binary_heads,
                                         verbose = TRUE)

# reticulate::py_help(object = heads$get_times)
times <- heads$get_times()[[1L]]

dat_heads <- stats::setNames(
  lapply(times, function(time) {
  heads$get_data(totim = time) }),
  nm = sprintf("time_%.1f", times)
  )

str(dat_heads)
summary(dat_heads$time_365.0)
summary(dat_heads$time_365.0)[c(1,6)]
# reticulate::py_help(object = flopy$modflow$mf$Modflow$load)
mf <- flopy$modflow$mf$Modflow$load(f = paths$mf_nam, 
                                    model_ws = dirname(paths$mf_nam),
                                    load_only = "dis", #c("dis", "bas6", "bcf6"),
                                    forgive = TRUE,
                                    verbose = TRUE
                                    )



# reticulate::py_help(object = flopy$utils$postprocessing$get_gradients)
gradients <- flopy$utils$postprocessing$get_gradients(
  heads = dat_heads$time_365.0,
  m = mf,
  nodata = c(-1e+30,1e+30))

str(gradients)
summary(gradients)

Export

Model to shapefile (gets really big !!!!)

mf$export(f = "MB16.shp")

For details on the used flopy functions have a look at the documentation here:

https://flopy.readthedocs.io/en/latest/source/flopy.utils.postprocessing.html

Flopy development

Install flopy development version from GitHub https://github.com/modflowpy/flopy

Install flopy-dev

### create an environment "flopy" for installing "flopy" 
kwb.python::conda_py_install(env_name = "flopy_dev", 
                             pkgs = list(conda = c("python", "pyshp"), 
                                         py = "https://github.com/modflowpy/flopy/zipball/develop/801a41705d8979c6982fb8c2955d56225d971218"))

Use flopy-dev

reticulate::use_miniconda(condaenv = "flopy_dev", required = TRUE)
flopy_dev <- reticulate::import("flopy", convert = TRUE)


kstpkper <- c(0,0)

# reticulate::py_help(object = flopy_dev$utils$postprocessing$get_extended_budget)
budget_array <- flopy_dev$utils$postprocessing$get_extended_budget(
      cbcfile = paths$binary_cell.by.cell, kstpkper = kstpkper) 

# reticulate::py_help(object = flopy_dev$utils$binaryfile$HeadFile)
heads <- flopy_dev$utils$binaryfile$HeadFile(filename = paths$binary_heads,
                                         verbose = TRUE)

# reticulate::py_help(object = heads$get_data)
heads_array <- heads$get_data(kstpkper = kstpkper)

# reticulate::py_help(object = flopy_dev$modflow$mf$Modflow$load)
mf <- flopy_dev$modflow$mf$Modflow$load(f = paths$mf_nam, 
                                    model_ws = dirname(paths$mf_nam),
                                    load_only = c("dis", "bas6", "bcf6"),
                                    forgive = TRUE,
                                    verbose = TRUE
                                    )


# reticulate::py_help(object = flopy$utils$postprocessing$get_specific_discharge)
# Works now with flopy 3.3.4 release candidate: 
# https://github.com/modflowpy/flopy/blob/b6092cbc94039750de309a832c7f554751782ac1/flopy/utils/postprocessing.py#L606
# reticulate::py_help(object = flopy_dev$utils$postprocessing$get_specific_discharge)
specific_discharge <- setNames(flopy_dev$utils$postprocessing$get_specific_discharge(
  vectors = budget_array, 
  model = mf, 
  head = heads_array
  ), 
  nm = c("qx", "qy", "qz")
)
# stats for all layers
summary(specific_discharge$qx)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -3837       0       0       0       0    1340  742928 
summary(specific_discharge$qy)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -3967       0       0       0       0    1744  743864 
summary(specific_discharge$qz)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#   -0.6     0.0     0.0     0.0     0.0     0.2  737240

#top layer
layer <- 1
summary(as.vector(specific_discharge$qx[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#   -0.6     0.0     0.0     0.0     0.0     0.2  737240
summary(as.vector(specific_discharge$qy[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -1.55   -0.03    0.00   -0.01    0.01    1.84   92983 
summary(as.vector(specific_discharge$qz[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -0.30    0.00    0.00    0.00    0.00    0.09   92155

# layer2
layer <- 2
summary(as.vector(specific_discharge$qx[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -0.23    0.00    0.00    0.00    0.00    0.28   92866
summary(as.vector(specific_discharge$qy[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -0.23    0.00    0.00    0.00    0.00    0.19   92983
summary(as.vector(specific_discharge$qz[layer,,]))
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -0.55    0.00    0.00    0.00    0.00    0.17   92155 


### Calculate flow velocity 
# Set effective porosity for each layer
# layer 1: 0.3 (unconfined)
# layer 2-8: 0.0001 (confined)
n_layers <- mf$nlay
n_layers_confined <- n_layers - 1
effective_porosity <- c(0.3, rep(0.0001, n_layers_confined))

flow_velocity <- setNames(lapply(specific_discharge, function(q) {
  for (layer in seq_len(n_layers)) {
    q[layer,,] <- q[layer,,] / effective_porosity[layer]
  }
  q
}), nm = c("vx", "vy", "vz")
)

# stats for all layers
summary(flow_velocity$vx)
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
#-38372483       -24         0       -29         8  13398450    742928 
summary(flow_velocity$vy)
#     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
#-39667686       -24         0       -44        11  17443782    743864 
summary(flow_velocity$vz)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#-5858.8    -1.2    -0.2    -1.5     0.0  1749.5  737240 

#top layer
layer <- 1
summary(as.vector(flow_velocity$vx[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -5.44   -0.06   -0.01   -0.02    0.01    5.89   92866
summary(as.vector(flow_velocity$vy[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -5.16   -0.08   -0.01   -0.03    0.02    6.13   92983 
summary(as.vector(flow_velocity$vz[layer,,])) 
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#  -1.01    0.00    0.00    0.00    0.00    0.29   92155

# layer2
layer <- 2
summary(as.vector(flow_velocity$vx[layer,,])) 
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
#-2288.47   -42.81    -0.34   -23.69     0.42  2842.85    92866 
summary(as.vector(flow_velocity$vy[layer,,])) 
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
#-2337.17   -43.87    -0.36   -42.00     1.08  1852.84    92983 
summary(as.vector(flow_velocity$vz[layer,,]))
#    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
#-5533.12    -3.01    -1.62    -4.79    -0.66  1749.48    92155


### Print as tidy data.frame 
flow_velocity_vx_long <- dwc.ar4gw::to_long(flow_velocity$vx)
head(flow_velocity_vx_long)
#  layer col row value
#1     1   1   1   NaN
#2     2   1   1   NaN
#3     3   1   1   NaN
#4     4   1   1   NaN
#5     5   1   1   NaN
#6     6   1   1   NaN

flow_velocity_vx_wide <- dwc.ar4gw::to_wide(flow_velocity_vx_long, 
                                            parameter = "vx")
flow_velocity_vx_wide
## A tibble: 534,640 x 10
#   column   row  vx_1  vx_2  vx_3  vx_4  vx_5  vx_6  vx_7  vx_8
#    <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1      1     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 2      2     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 3      3     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 4      4     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 5      5     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 6      6     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 7      7     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 8      8     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
# 9      9     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
#10     10     1   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
## ... with 534,630 more rows

### Plot data of all layers
dwc.ar4gw::plot_data(heads_array, title = "Heads")
dwc.ar4gw::plot_data(flow_velocity$vx, title = "Flow velocity x")
dwc.ar4gw::plot_data(flow_velocity$vy, title = "Flow velocity y")
dwc.ar4gw::plot_data(flow_velocity$vz, title = "Flow velocity z")

KWB-R/dwc.ar4gw documentation built on Dec. 18, 2021, 2:35 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

KWB-R/dwc.ar4gw
R Package for Preparing Modflow Output Data for Articifical Reality Visualisation

In KWB-R/dwc.ar4gw: R Package for Preparing Modflow Output Data for Articifical Reality Visualisation

Install

Define paths

Modflow Files

Locally

Cloud

Flopy release

Install flopy

Use flopy

Export

Flopy development

Install flopy-dev

Use flopy-dev

R Package Documentation

Browse R Packages

We want your feedback!

KWB-R/dwc.ar4gw R Package for Preparing Modflow Output Data for Articifical Reality Visualisation

In KWB-R/dwc.ar4gw: R Package for Preparing Modflow Output Data for Articifical Reality Visualisation

Install

Define paths

Modflow Files

Locally

Cloud

Flopy release

Install flopy

Use flopy

Export

Flopy development

Install flopy-dev

Use flopy-dev

R Package Documentation

Browse R Packages

We want your feedback!

KWB-R/dwc.ar4gw
R Package for Preparing Modflow Output Data for Articifical Reality Visualisation