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)
# 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)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 MB16, which is defined in the paths_list in the previous chapter Define Paths.
Copy the content of the DWC cloud folder dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005 to C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005 as shown below:
#> $tdir
#> [1] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z"
#>
#> $cloud_dir
#> [1] "dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005"
#>
#> $local_dir
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005"
#>
#> $download_dir
#> [1] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z"
#>
#> $modflow_dir
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005"
#>
#> $model_name
#> [1] "MB16"
#>
#> $mfout
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16"
#>
#> $mf_cbc
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.cbc"
#>
#> $mf_nam
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.nam"
#>
#> $mf_bas
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bas"
#>
#> $mf_dis
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.dis"
#>
#> $text_heads
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.fhd"
#>
#> $text_drawdown
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.fdn"
#>
#> $binary_cell.by.cell
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.cbc"
#>
#> $binary_heads
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bhd"
#>
#> $binary_drawdown
#> [1] "C:/kwb/projects/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bdn"
paths <- paths_localYou can download the required model files from the DWC cloud if you are a registered user with access to the folder dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005
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 belowFor 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
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))
#> Listing
#> Listing dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005
mf_files
#> file isdir etag lastmodified fileid
#> 1 MB16.bas FALSE 093903a16e8c3ab272797c668353ef35 2021-03-04 13:44:40 212524
#> 2 MB16.bcf FALSE 9acd14724ce9a5cc54de48521b007f7c 2021-03-04 13:45:58 212582
#> 3 MB16.bhd FALSE b0127f18ceefcdd20be207967dc1faba 2021-03-04 10:09:48 212327
#> 4 MB16.cbc FALSE 3d41e58d4a0fe2c06f9a3961ceff7377 2021-03-04 13:52:24 212420
#> 5 MB16.dis FALSE 569dcc1d7791c49e8322649a3b4e31b8 2021-03-04 13:43:34 212534
#> 6 MB16.evt FALSE 3031bc06858c650f5e0e1321e4973a8e 2021-03-04 13:47:58 212360
#> 7 MB16.gmg FALSE d4c09e5e0a48b28c0c5c04e3a7812834 2021-03-04 13:44:42 212258
#> 8 MB16.lst FALSE 239130516d5832c76785f5c14b34e6f1 2021-03-04 13:52:32 212590
#> 9 MB16.mpbas FALSE bbb090ff9645c018d6f16961d0cd7fe9 2021-03-04 13:48:32 212369
#> 10 MB16.mplst FALSE e664a1f05d7e84ff640ac775b7137a04 2021-03-04 13:59:22 212279
#> 11 MB16.nam FALSE 3350173dd364b9b8975e81e010ed83a1 2021-03-04 13:47:58 212262
#> 12 MB16.oc FALSE 2915cae1fa39614223619a5cd8101c56 2021-03-04 13:44:42 212268
#> 13 MB16.rch FALSE 068b37cdace113838de0c890c68ea350 2021-03-04 13:47:22 212326
#> 14 MB16.riv FALSE 2d7c2534f2673e2cb1f3cde33d803d44 2021-03-04 13:47:06 212301
#> 15 MB16.wel FALSE a8b3e2bf6ef21948ca90738f2c0c70d8 2021-03-04 13:45:58 212297
#> permissions size
#> 1 RMGDNVW 121481684
#> 2 RMGDNVW 178045898
#> 3 RMGDNVW 17108832
#> 4 RMGDNVW 60192496
#> 5 RMGDNVW 106860720
#> 6 RMGDNVW 35609356
#> 7 RMGDNVW 387
#> 8 RMGDNVW 275034643
#> 9 RMGDNVW 26524721
#> 10 RMGDNVW 1095
#> 11 RMGDNVW 533
#> 12 RMGDNVW 333
#> 13 RMGDNVW 11869894
#> 14 RMGDNVW 3396530
#> 15 RMGDNVW 492728
#> href
#> 1 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bas
#> 2 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bcf
#> 3 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bhd
#> 4 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.cbc
#> 5 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.dis
#> 6 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.evt
#> 7 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.gmg
#> 8 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.lst
#> 9 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.mpbas
#> 10 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.mplst
#> 11 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.nam
#> 12 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.oc
#> 13 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.rch
#> 14 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.riv
#> 15 /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.wel
kwb.nextcloud::download_files(href = mf_files$href,
target_dir = paths$download_dir)
#> Splitting paths ... ok. (0.00s)
#> Removing the first 11 path segments ... ok. (0.00s)
#> Putting path segments together ... ok. (0.00s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bas ... ok. (16.92s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bcf ... ok. (23.55s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.bhd ... ok. (2.94s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.cbc ... ok. (9.18s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.dis ... ok. (15.07s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.evt ... ok. (5.90s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.gmg ... ok. (0.42s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.lst ... ok. (38.65s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.mpbas ... ok. (4.14s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.mplst ... ok. (0.46s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.nam ... ok. (0.40s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.oc ... ok. (0.47s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.rch ... ok. (2.31s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.riv ... ok. (0.86s)
#> Downloading /remote.php/dav/files/mrustl/dwc/Work-Packages/WP3_Public_Involvement/AR4GW/ModelData/MF2005/MB16.wel ... ok. (0.92s)
#> [1] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.bas"
#> [2] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.bcf"
#> [3] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.bhd"
#> [4] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.cbc"
#> [5] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.dis"
#> [6] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.evt"
#> [7] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.gmg"
#> [8] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.lst"
#> [9] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.mpbas"
#> [10] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.mplst"
#> [11] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.nam"
#> [12] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.oc"
#> [13] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.rch"
#> [14] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.riv"
#> [15] "/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T//RtmpeYyw8Z/MB16.wel"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))
#> * Downloading 'https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh' ...
#> * Installing Miniconda -- please wait a moment ...
#> * Miniconda has been successfully installed at '/Users/runner/Library/r-miniconda'.
#> python: /Users/runner/Library/r-miniconda/envs/flopy/bin/python
#> libpython: /Users/runner/Library/r-miniconda/envs/flopy/lib/libpython3.9.dylib
#> pythonhome: /Users/runner/Library/r-miniconda/envs/flopy:/Users/runner/Library/r-miniconda/envs/flopy
#> version: 3.9.2 | packaged by conda-forge | (default, Feb 21 2021, 05:02:20) [Clang 11.0.1 ]
#> numpy: /Users/runner/Library/r-miniconda/envs/flopy/lib/python3.9/site-packages/numpy
#> numpy_version: 1.20.2
#>
#> NOTE: Python version was forced by use_python function
#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)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -1367.3788 -3.2607 -0.0040 -2.5185 0.0034 1521.7942
str(extended_budget)
#> List of 1
#> $ time_1.0:List of 3
#> ..$ Qx_ext: num [1:8, 1:656, 1:816] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ Qy_ext: num [1:8, 1:657, 1:815] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ Qz_ext: num [1:9, 1:656, 1:815] 0 0 0 0 0 0 0 0 0 0 ...
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)
#> List of 1
#> $ time_365.0: num [1:8, 1:656, 1:815] 1e+30 1e+30 1e+30 1e+30 1e+30 ...
summary(dat_heads$time_365.0)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -1.000e+30 3.200e+01 3.500e+01 1.646e+29 5.700e+01 1.000e+30
summary(dat_heads$time_365.0)[c(1,6)]
#> Min. Max.
#> -1e+30 1e+30
# 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)
#> num [1:7, 1:656, 1:815] 0 0 0 0 0 0 0 0 0 0 ...
summary(gradients)
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#> -1.3690 0.0000 0.0000 0.0143 0.0002 2.5653 1102Model 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
Install flopy development version from GitHub https://github.com/modflowpy/flopy
### 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"))
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)
effective_porosity <- c(0.3, rep(0.0001, 7))
n_layers <- mf$nlay
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")