rogiersbart / rmodflow Goto Github PK

When using the function rmf_convert_grid_to_xyz witj i and j arguments, I don't get correct x, y position. Here is the code example

library(tidyverse)
library(RMODFLOW)
library(sf)
prjfile <- here('RMODFLOW/cat-sm_p.prj')
c('+init=epsg:31370','193951.5 211439.7','357.2') %>% 
  write_lines(prjfile)
prj_p <- prjfile %>% 
  rmf_read_prj()

dis_p <- rmf_create_dis(nlay = 6,
                        nrow = 82,
                        ncol = 165,
                        itmuni = 4,
                        lenuni = 2,
                        delr = rep(50,165),
                        delc = rep(50,82))
xyl <- expand.grid(cumsum(dis_p$delr)-dis_p$delr/2,sum(dis_p$delc)-(cumsum(dis_p$delc)-dis_p$delc/2)) %>% 
  as_tibble() %>% 
  magrittr::set_colnames(c('x','y')) %>% 
  mutate(id=1:(dis_p$nrow*dis_p$ncol),i=rep(1:dis_p$nrow,each=dis_p$ncol),j=rep(1:dis_p$ncol,dis_p$nrow))
xyg_p <- rmf_convert_grid_to_xyz(x=xyl$x,y=xyl$y,prj=prj_p) %>% 
  as_tibble() %>% 
  mutate(id=1:(dis_p$nrow*dis_p$ncol)) %>%
  inner_join(xyl %>% select(id,i,j),by="id") %>% 
  st_as_sf(coords=c('x','y'),crs=31370) %>% 
  st_join(sm_bound,join = st_intersects) %>% 
  mutate(active=ifelse(is.na(id.y),0,1)) %>% 
  select(id=id.x,i,j,active)
xyg_p %>% write_sf(here('data/shapes/grid_parent_point.gpkg'))

Now, when running

From the current namespace file, it is obvious we are not really using a consistent approach for naming our S3 classes. We have at least these: bas, bcf, bud, cbc, chd, ddn, de4, dis, drn, evt, ghb, gmg, hed, hfb, hob, hpr, huf, kdep, lmt, lpf, lvda, mlt, modflow, nam, nwt, oc, pcg, pval, rch, riv, rmf_2d_array, rmf_3d_array, rmf_4d_array, rmf_list, rmf_package, sip, upw, wel, zon. I'd like to propose adding the rmf_ prefix everywhere for consistency, and minimizing the potential for conflicts with other packages. Also, I don't think we really have methods in use for rmf_package, but I would think having this there for all objects that can be part of the first level of a full model object (currently modflow I believe, but maybe this should be rmf_model?) might be useful at some point (i.e. we always have the option to use S3 inheritance when it is there). Besides that one, I think it might be useful to introduce an rmf_object one, that is added after e.g. rmf_dis and rmf_package, and serves as parent class for all RMODFLOW S3 objects. This would allow for consistent printing for instance (I'm thinking of using rui::inspect() there, as I would prefer moving the current printing methods to summary methods), by defining a single rmf_object S3 method.

I just noted the current development branch contains a change to many (if not all) of the package object structures: The parameters related to dimensions are now exposed as object$dimensions$np instead of object$np for instance. That seems to be quite a breaking change, which I should have noticed before, and is actually breaking some of my RMODFLOW code.

I am prepared to adapt, but we need to discuss this I believe. My philosophy was, any dataset/variable defined in the online guide, which can be easily exposed to the user as object$name should be exposed in such a way. In case of MODFLOW parameters, I believe it indeed becomes too difficult, but plain integers as hob$nh, hob$mobs, hob$maxm should still be available like that if you ask me. This seems much more user friendly.

Any thoughts?

The idea has surfaced several times:

• RMODFLOW.examples in #6, could also serve as a showcase of real-life models
• RMODFLOW.invert in #7
• Testing in #8 could potentially also be tackled with an extension package, as ideally many real-life example models would be available. These will however not be included in the RMODFLOW package itself. Functions for installing (i.e. downloading and putting them in a predictable location) them could be part of RMODFLOW.examples, and RMODFLOW.test could install all these with .onLoad() for instance)

Reading/writing rmf_list structures has been updated in commits cebe60d and 15d3ad8 (develop branch) using the capabilities of readr. This resulted in a drastic speed-up when reading/writing packages containing MODFLOW list structures. Similarly, we should look at speeding-up reading/writing arrays using readr.

RMODFLOW (development branch) produces erroneous input to hob package and subsequently to the written hob file in case of transient observations. In case there is only one observation in a time series, the nrefsp is not assigned correctly. This leads to consequent errors. Also, when writing the hob package file to disk, the nrefsp remains always the same.

Here is the reproducible example:

library(tidyverse)
library(RMODdev) #The development branch of RMODFLOW
#> ! {RMODdev} is still in its experimental lifecycle stage.�[K
#> ! Use at your own risk, and submit issues here:�[K
#> ! <https://github.com/rogiersbart/RMODdev/issues>�[K
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

loc <- tribble(~name, ~top, ~bottom,  ~x,  ~y,
                "P1",   -2,     - 5, 3.5, 4.2,
                "P2",   -1,     -20, 8.2, 9.3)

ts <- tribble(~name, ~time, ~head,
               "P1",   4.4,  20.7,
               "P2",   0.5,  23.6,
               "P2",   1.5,  23.8,
               "P2",   2.5,  23.8,
               "P2",   3.9,  23.6,
               "P2",   4.8,  23.4)

dis <- rmf_create_dis(nper=5)
hob <- rmf_create_hob(locations = loc,time_series = ts,dis = dis)
hob #Note nrefsp for P1 observation that equals 5, should be 1
#> RMODdev Head-Observation Package object with: 
#> 6 head observations of which 5 are multilayer observations 
#> Maximum number of layers used for multilayer observations is 2 
#> 
#> Observed values and simulated equivalents are written to the head-prediction file on unit number 665 
#> Simulated equivalents of dry cells are assigned a value of -888 in the head-prediction file 
#> Input and output data are not printed to the listing file 
#> Time-offset multiplier: 1 
#> 
#> Observations overview: 
#>   obsnam layer           pr row column nrefsp irefsp itt toffset  roff   coff
#> 1     P1     1            1  10      1      5      5   1     0.4 0.458 -0.465
#> 2     P2  1, 2 0.473684....  10      1      5      1   1     0.5 0.407 -0.418
#> 3     P2  1, 2 0.473684....  10      1      5      2   1     0.5 0.407 -0.418
#> 4     P2  1, 2 0.473684....  10      1      5      3   1     0.5 0.407 -0.418
#> 5     P2  1, 2 0.473684....  10      1      5      4   1     0.9 0.407 -0.418
#> 6     P2  1, 2 0.473684....  10      1      5      5   1     0.8 0.407 -0.418
#>   hobs
#> 1 20.7
#> 2 23.6
#> 3 23.8
#> 4 23.8
#> 5 23.6
#> 6 23.4

hob %>% rmf_write_hob("temp.hob")
cat(readLines("temp.hob"),sep = "\n") #See line of P1 having -5 instead of 1
#> # MODFLOW Head-Observation Package created by RMODdev, version 0.4.0.9000 
#> # 
#> 6 5 2 665 -888
#> 1
#> P1 1 10 1 -5 0.4 0.458 -0.465 20.7
#> 1
#> P1 5 0.4 20.7
#> P2 1 0.5 23.6
#> P2 2 0.5 23.8
#> P2 3 0.5 23.8
#> P2 4 0.9 23.6
#> P2 -2 10 1 -5 0.8 0.407 -0.418 23.4
#> 1 0.473684210526316 2 0.526315789473684
#> 1
#> P2 5 0.8 23.4
#>    
#>    
#>    
#> 

hob <- rmf_read_hob("temp.hob")
hob #NA's produced
#> RMODdev Head-Observation Package object with: 
#> 6 head observations of which 5 are multilayer observations 
#> Maximum number of layers used for multilayer observations is 2 
#> 
#> Observed values and simulated equivalents are written to the head-prediction file on unit number 665 
#> Simulated equivalents of dry cells are assigned a value of -888 in the head-prediction file 
#> Input and output data are not printed to the listing file 
#> Time-offset multiplier: 1 
#> 
#> Observations overview: 
#>    obsnam layer           pr row column nrefsp irefsp itt toffset  roff   coff
#> 1      P1     1            1  10      1      5      5   1     0.4 0.458 -0.465
#> 2      P2     1            1  10      1      5      1   1     0.5 0.458 -0.465
#> 3      P2     1            1  10      1      5      2   1     0.5 0.458 -0.465
#> 4      P2     1            1  10      1      5      3   1     0.5 0.458 -0.465
#> 5      P2     1            1  10      1      5      4   1     0.9 0.458 -0.465
#> 6      P2  1, 2 0.473684....  10      1      5      5   1     0.8 0.407 -0.418
#> 7    <NA>  1, 2 0.473684....  10      1      5     NA   1      NA 0.407 -0.418
#> 8    <NA>  1, 2 0.473684....  10      1      5     NA   1      NA 0.407 -0.418
#> 9    <NA>  1, 2 0.473684....  10      1      5     NA   1      NA 0.407 -0.418
#> 10   <NA>  1, 2 0.473684....  10      1      5     NA   1      NA 0.407 -0.418
#>    hobs
#> 1  20.7
#> 2  23.6
#> 3  23.8
#> 4  23.8
#> 5  23.6
#> 6  23.4
#> 7    NA
#> 8    NA
#> 9    NA
#> 10   NA

^{Created on 2021-02-12 by the reprex package (v0.3.0)}

Find latest versions of the MODFLOW variants on the USGS webpages, download, unzip to C:/WRDAPP, and enable finding the paths to the executables from R.

Currently we rely on hydroGOF for the model performance evaluation. Weights are however not used there, but they are common enough in groundwater modelling to implement I think. Question is what to do: implement things inside RMODFLOW, or suggest enhancement of hydroGOF maybe? Will submit an issue there to take this up with the author (see hzambran/hydroGOF#11).

Concerning where to get the weights from: MODFLOW-2000 has support for this in the input files. I'm not sure what the 2005 variants do when the extra columns are there however. This has to be checked, and if not possible through the input files, the rmf_optimize API should be extended to allow for specification of the weights I suppose.

ModelMuse allows specifying these weight/statistic parameters, but doens't write them to the input files of the 2005 variants. However, when setting up UCODE, it does include this I believe, so the functionality is there.

Currently we have some functionality for metadata in a separate *.prj projection file, possibly containing crs, origin coordinates, start time and grid rotation. Some of these parameters can be part of dataset 1 in the discretization file, with MODFLOW_OWHM. ModelMuse currently also includes the xy coordinates of the upper and lower, left and right corners, and the grid angle in the dataset 0 in the discretization file, and has been exporting usgs.model.reference files with the crs for quite some time now.

We need to figure out if we still want a separate metadata file (if so, I would opt for yaml, or the ModelMuse usgs.model.reference format), or if we want to include everything in the discretization file. The latter would require enabling processing of the ModelMuse comments in dataset 0, reading usgs.model.reference files if they exist, as well as the MODFLOW_OWHM dataset 1 options (and a rule for handling possible conflicts I guess), and a function to easily inject these details in discretization files that do not contain them.

Currently, RMODFLOW imports the akima package for interpolation to a regular grid. This is only used in rmf_plot when type = 'contour' to interpolate the MODFLOW grid values (which are possibly irregular) to a regular grid so that ggplot2::stat_contour can handle these (see also its documentation).

However, as pointed out here, the akima package uses a non-commercial license. I've tried to replace the akima functionality in RMODFLOW with the interp package but the results are not the same, i.e. using interp seems to (sometimes) break rmf_plot eventhough interp::interp is designed as a replacement for akima::interp.

Any ideas on how to tackle this?

At the moment, almost no RMODFLOW functions have examples provided in their documentation. I think it's important to add these and to also make sure all new functions should include examples.
This will not only make the intent of the function clearer but will also serve as a first step in the testing process since all examples from (exported) functions are run by R CMD check (as called by e.g. devtools::check)

Let's create a place to also discuss this. My point of view at the moment is that RMODFLOW should follow ModelMuse conventions as closely as possible. Main argument for that is that the package should in this case be a more welcoming environment for the ModelMuse user. Given the GUI aspect, I think it is still more likely people are introduced to MODFLOW through ModelMuse than through RMODFLOW, FloPy or mfLab for instance.
There's two sides to this I think:

1. File names: It is very useful we adopt the file extension conventions from ModelMuse I think (i.e. for our function names), as it simplifies figuring out what functions have to be used to handle specific files. I mentioned in #8 we should aim at keeping the rmf_read_fhd and rmf_read_bhd functions for that reason. I would even go further and suggest to remove rmf_read_hed for instance, and replace it by one of these two. There's no reason for introducing a separate RMODFLOW function name not based on any file extension convention.
2. Functionality: Certainly at this stage, RMODFLOW should not aim at doing things that are not part of the ModelMuse scope. So standard calibration through optim based on pval parameters should be there I think. More complex methods (based on pval parameters, or with custom parameters set up in R) should be moved to an extension package.

Currently rmf_install() (seems to be broken, and) puts the codes in C:/WRDAPP. I would like to revert back to installing in the system.file(package = "RMODFLOW") folder to avoid interfering with the user file system, and make the package behave more friendly under other OSs.

Furthermore, I'd like the function to "install" example models as well. I included some ModelMuse examples in the past, but I think it is more appropriate if these are not installed at package installation, but rather when it is asked for explicitly. I also thought about including examples in a separate package (e.g. RMODFLOW.examples), but I'd rather do it through this function, as we could then probably also make use of such example models in the function documentation and vignettes.

The idea is to have for all packages/files from MODFLOW-2005 based codes the following functions:

• rmf_read_* for reading,
• rmf_write_* for writing,
• rmf_create_* for creating the object from scratch, and I was thinking about
• rmf_build_* (or maybe another more appropriate verb?) for building the object from spatial and/or tidy data (e.g. an sf object, raster or data frame)

Note @cneyens mentioned some ideas related to the latter two in cneyens#1.

S3 methods for visualization, printing etc is another thing. But these basics should be there in v0.5.0 I think. An exhaustive list of all possible packages/files and status of the corresponding functions might be useful here.

I'd like to modify the api for executing MODFLOW. Currently, we have rmf_run_modflow(), rmf_run_opt(), rmf_run_sen() and rmf_run_rsm() for executing the model, optimization, sensitivity analysis and response surface mapping respectively. I think replacing the first three with rmf_execute(), rmf_optimize() and rmf_analyze would make things more clear. Furthermore, I think the optimization functionality should be reduced to standard local optimization only (with bounds, so only L-BGFS-B probably). The rest should at one point go to a separate package (e.g. RMODFLOW.invert). This relates to the idea of compatibility with, and limiting RMODFLOW functionality to ModelMuse (and friends ModelMate, ModelViewer, GW_Chart, ListingAnalyst). This package could then also include something like the current rmf_run_rsm(), and much more exotic things probably in the long run.

I'm getting the following errors when reading the attached .dis file. Can you help me? Thanks!

dis <- read_dis(paste0("example", ".dis"))
Error in gsub(paste0("([[:alnum:][:punct:][:space:]]{", nPerNum, "})"), :
invalid regular expression '([[:alnum:][:punct:][:space:]]{NA})', reason 'Invalid contents of {}'
In addition: Warning messages:
1: NAs introduced by coercion
2: In read_modflow_array(dis_lines, 1, dis$nrow, 1, ndim = 1) :
NAs introduced by coercion

Session Info

devtools::session_info()
Session info ----------------------------------------------------------------------------------------------------------------------------
setting value
version R version 3.4.0 (2017-04-21)
system x86_64, linux-gnu
ui RStudio (1.0.143)
language (EN)
collate en_US.UTF-8
tz America/Los_Angeles
date 2017-04-28
Packages --------------------------------------------------------------------------------------------------------------------------------
package * version date source
acepack 1.4.1 2016-10-29 cran (@1.4.1)
akima 0.6-2 2016-12-20 cran (@0.6-2)
animation 2.5 2017-03-30 cran (@2.5)
automap 1.0-14 2013-08-29 cran (@1.0-14)
backports 1.0.5 2017-01-18 cran (@1.0.5)
base64enc 0.1-3 2015-07-28 cran (@0.1-3)
checkmate 1.8.2 2016-11-02 cran (@1.8.2)
class 7.3-14 2015-08-30 CRAN (R 3.2.2)
cluster 2.0.6 2017-03-16 CRAN (R 3.3.3)
colorspace 1.3-2 2016-12-14 cran (@1.3-2)
curl 2.5 2017-04-14 CRAN (R 3.4.0)
data.table 1.10.4 2017-02-01 cran (@1.10.4)
DEoptim 2.2-4 2016-12-19 cran (@2.2-4)
devtools * 1.12.0 2016-12-05 CRAN (R 3.4.0)
digest 0.6.4 2013-12-03 CRAN (R 3.0.2)
directlabels 2017.03.31 2017-04-08 cran (@2017.03)
e1071 1.6-8 2017-02-02 cran (@1.6-8)
FNN 1.1 2013-07-31 cran (@1.1)
foreign 0.8-67 2016-09-13 CRAN (R 3.3.1)
Formula 1.2-1 2015-04-07 cran (@1.2-1)
ggplot2 2.2.1 2016-12-30 cran (@2.2.1)
git2r 0.18.0 2017-01-01 CRAN (R 3.4.0)
gridExtra 2.2.1 2016-02-29 cran (@2.2.1)
gstat 1.1-5 2017-03-12 cran (@1.1-5)
gtable 0.2.0 2016-02-26 cran (@0.2.0)
Hmisc 4.0-2 2016-12-31 cran (@4.0-2)
hms 0.3 2016-11-22 cran (@0.3)
htmlTable 1.9 2017-01-26 cran (@1.9)
htmltools 0.3.5 2016-03-21 CRAN (R 3.4.0)
htmlwidgets 0.8 2016-11-09 cran (@0.8)
httpuv 1.3.3 2015-08-04 CRAN (R 3.4.0)
httr 1.2.1 2016-07-03 CRAN (R 3.4.0)
hydroGOF 0.3-8 2014-02-04 cran (@0.3-8)
hydroPSO 0.3-4 2014-04-13 cran (@0.3-4)
hydroTSM 0.4-2-1 2014-01-23 cran (@0.4-2-1)
intervals 0.15.1 2015-08-27 cran (@0.15.1)
jsonlite 1.4 2017-04-08 CRAN (R 3.4.0)
knitr 1.15.1 2016-11-22 cran (@1.15.1)
lattice 0.20-35 2017-03-25 CRAN (R 3.3.3)
latticeExtra 0.6-28 2016-02-09 cran (@0.6-28)
lazyeval 0.2.0 2016-06-12 cran (@0.2.0)
lhs 0.14 2016-08-09 cran (@0.14)
magrittr 1.5 2014-11-22 cran (@1.5)
Matrix 1.2-8 2017-01-20 CRAN (R 3.3.2)
memoise 1.1.0 2017-04-21 CRAN (R 3.4.0)
mime 0.5 2016-07-07 CRAN (R 3.4.0)
munsell 0.4.3 2016-02-13 cran (@0.4.3)
nnet 7.3-12 2016-02-02 CRAN (R 3.2.3)
plyr 1.8.4 2016-06-08 cran (@1.8.4)
quadprog 1.5-5 2013-04-17 cran (@1.5-5)
R6 2.2.0 2016-10-05 CRAN (R 3.4.0)
RColorBrewer 1.1-2 2014-12-07 cran (@1.1-2)
Rcpp 0.12.10 2017-03-19 cran (@0.12.10)
readr 1.1.0 2017-03-22 cran (@1.1.0)
reshape 0.8.6 2016-10-21 cran (@0.8.6)
rgdal 1.2-7 2017-04-25 CRAN (R 3.4.0)
rgl 0.98.1 2017-03-08 cran (@0.98.1)
RMODFLOW * 0.4.0 2017-04-28 Github (8a9578d)
rpart 4.1-11 2017-04-21 CRAN (R 3.4.0)
scales 0.4.1 2016-11-09 cran (@0.4.1)
shiny 1.0.3 2017-04-26 CRAN (R 3.4.0)
sp 1.2-4 2016-12-22 cran (@1.2-4)
spacetime 1.2-0 2016-09-03 cran (@1.2-0)
stringi 1.1.5 2017-04-07 cran (@1.1.5)
stringr 1.2.0 2017-02-18 cran (@1.2.0)
survival 2.41-3 2017-04-04 CRAN (R 3.3.3)
tibble 1.3.0 2017-04-01 cran (@1.3.0)
withr 1.0.2 2016-06-20 CRAN (R 3.4.0)
xtable 1.8-2 2016-02-05 cran (@1.8-2)
xts 0.9-7 2013-08-29 CRAN (R 3.0.1)
zoo 1.8-0 2017-04-12 cran (@1.8-0)

example.dis.tar.gz

-- Tung

Hey RMODFLOW team,

I'm a big fan of this package. Is there any plan to update it to handle the recent deprecation of RGDAL? The dependency seems to prevent RMODFLOW from being installed on the latest version(s) of R. Thanks!

This concerns rmf_read_rch in the development branch of cneyens repository. When reading in a complex recharge file, the code fails. I think the subroutine rmfi_parse_variables does not read all entries, when INSTANCES are used. Here are the files:

rch_related_files.zip