Calibration user exp. bug fixed. Doc updated

krige.head.calib outputs changed to a list variable defining the entire mapping method and inputs. Inputting this to krige.head will map the head without any of the previously required handling of calibrated parameter values.  See ?krige.head.calib for details.

.Rproj.user/19E94E50/pcs/source-pane.pper

@@ -1,3 +1,3 @@
 {
-    "activeTab": 2
+    "activeTab": 0
 }

.Rproj.user/19E94E50/pcs/windowlayoutstate.pper

@@ -1,12 +1,12 @@
 {
     "left": {
-        "splitterpos": 313,
+        "splitterpos": 320,
         "topwindowstate": "NORMAL",
         "panelheight": 1079,
         "windowheight": 1117
     },
     "right": {
-        "splitterpos": 490,
+        "splitterpos": 501,
         "topwindowstate": "NORMAL",
         "panelheight": 1079,
         "windowheight": 1117

.Rproj.user/19E94E50/persistent-state

@@ -1,6 +1,6 @@
 build-last-errors="[]"
 build-last-errors-base-dir="C:/Users/tpet0008/Documents/HydroMap/"
-build-last-outputs="[{\"type\":0,\"output\":\"==> Rcmd.exe INSTALL --preclean --no-multiarch --with-keep.source HydroMap\\n\\n\"},{\"type\":1,\"output\":\"* installing to library 'C:/Program Files/R/R-4.0.3/library'\\r\\n\"},{\"type\":1,\"output\":\"* installing *source* package 'HydroMap' ...\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** using staged installation\\r\\n\"},{\"type\":1,\"output\":\"** R\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** data\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** moving datasets to lazyload DB\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"  converting help for package 'HydroMap'\\r\\n\"},{\"type\":1,\"output\":\"    import.DEM                              html  \"},{\"type\":1,\"output\":\"*** installing help indices\\r\\n\"},{\"type\":1,\"output\":\"    finding HTML links ... done\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"    import.pointData                        html  \\r\\n\"},{\"type\":1,\"output\":\"    krige.head                              html  \"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"    krige.head.calib                        html  \"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** building package indices\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\r\\n\"},{\"type\":1,\"output\":\"* DONE (HydroMap)\\r\\n\"},{\"type\":1,\"output\":\"\"}]"
+build-last-outputs="[{\"type\":0,\"output\":\"==> Rcmd.exe INSTALL --preclean --no-multiarch --with-keep.source HydroMap\\n\\n\"},{\"type\":1,\"output\":\"* installing to library 'C:/Program Files/R/R-4.0.3/library'\\r\\n\"},{\"type\":1,\"output\":\"* installing *source* package 'HydroMap' ...\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** using staged installation\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** R\\r\\n\"},{\"type\":1,\"output\":\"** data\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** moving datasets to lazyload DB\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"  converting help for package 'HydroMap'\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** installing help indices\\r\\n\"},{\"type\":1,\"output\":\"    finding HTML links ...\"},{\"type\":1,\"output\":\"    import.DEM                              html  \\r\\n\"},{\"type\":1,\"output\":\"    import.pointData                        html  \\r\\n\"},{\"type\":1,\"output\":\"    krige.head                              html  \"},{\"type\":1,\"output\":\" done\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"    krige.head.calib                        html  \"},{\"type\":1,\"output\":\"\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** building package indices\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\r\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"* DONE (HydroMap)\\r\\n\"},{\"type\":1,\"output\":\"\"}]"
 compile_pdf_state="{\"tab_visible\":false,\"running\":false,\"target_file\":\"\",\"output\":\"\",\"errors\":[]}"
 files.monitored-path=""
 find-in-files-state="{\"handle\":\"\",\"input\":\"\",\"path\":\"\",\"regex\":false,\"ignoreCase\":false,\"results\":{\"file\":[],\"line\":[],\"lineValue\":[],\"matchOn\":[],\"matchOff\":[],\"replaceMatchOn\":[],\"replaceMatchOff\":[]},\"running\":false,\"replace\":false,\"preview\":false,\"gitFlag\":false,\"replacePattern\":\"\"}"

.Rproj.user/19E94E50/sources/prop/0072A15F

@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "202,24",
-    "scrollLine": "184"
+    "cursorPosition": "26,60",
+    "scrollLine": "11"
 }

.Rproj.user/19E94E50/sources/prop/0243BC54

@@ -0,0 +1,6 @@
+{
+    "source_window_id": "",
+    "Source": "Source",
+    "cursorPosition": "26,2",
+    "scrollLine": "0"
+}

.Rproj.user/19E94E50/sources/prop/06EF12A0

@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "4,4",
+    "cursorPosition": "6,90",
     "scrollLine": "0"
 }

.Rproj.user/19E94E50/sources/prop/57B75005

@@ -0,0 +1,6 @@
+{
+    "source_window_id": "",
+    "Source": "Source",
+    "cursorPosition": "2,2",
+    "scrollLine": "0"
+}

.Rproj.user/19E94E50/sources/prop/AF65AB58

@@ -2,5 +2,5 @@
     "source_window_id": "",
     "Source": "Source",
     "cursorPosition": "109,71",
-    "scrollLine": "88"
+    "scrollLine": "0"
 }

.Rproj.user/19E94E50/sources/prop/BD952B7F

@@ -1,5 +1,5 @@
 {
-    "cursorPosition": "2,12",
+    "cursorPosition": "2,16",
     "scrollLine": "0",
     "Source": "Source"
 }

.Rproj.user/19E94E50/sources/prop/D3423677

@@ -1,6 +1,6 @@
 {
     "source_window_id": "",
     "Source": "Source",
-    "cursorPosition": "292,2",
-    "scrollLine": "0"
+    "cursorPosition": "93,8",
+    "scrollLine": "69"
 }

.Rproj.user/19E94E50/sources/prop/INDEX

@@ -4,6 +4,8 @@ C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.allData.R="E8D41511"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.objFunc.R="D3423677"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.pointData.R="00FD9ABF"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.smoothedDEM.R="FBFA4863"
+C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.variogram.R="0243BC54"
+C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fget.variogramModel.R="57B75005"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fkrige.head.R="0072A15F"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fkrige.head.calib.R="06EF12A0"
 C%3A%2FUsers%2Ftpet0008%2FDocuments%2FHydroMap%2FR%2Fkrige.head.crossval.R="AF65AB58"

.Rproj.user/19E94E50/sources/per/t/4C0D952C → .Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C

@@ -3,7 +3,7 @@
     "path": "C:/Users/tpet0008/Documents/HydroMap/R/krige.head.R",
     "project_path": "R/krige.head.R",
     "type": "r_source",
-    "hash": "2758867917",
+    "hash": "1487510892",
     "contents": "",
     "dirty": false,
     "created": 1683426418056.0,
@@ -12,15 +12,15 @@
     "properties": {
         "source_window_id": "",
         "Source": "Source",
-        "cursorPosition": "202,24",
-        "scrollLine": "184"
+        "cursorPosition": "26,60",
+        "scrollLine": "11"
     },
     "folds": "",
-    "lastKnownWriteTime": 1683775888,
+    "lastKnownWriteTime": 1683860925,
     "encoding": "UTF-8",
     "collab_server": "",
     "source_window": "",
-    "last_content_update": 1683775888647,
+    "last_content_update": 1683860925272,
     "read_only": false,
     "read_only_alternatives": []
 }

.Rproj.user/19E94E50/sources/per/t/4C0D952C-contents → .Rproj.user/19E94E50/sources/session-4dda64d8/4C0D952C-contents

@@ -10,18 +10,12 @@
 #' these features can be individually controlled by the user.
 #'
 #' Importantly, if the mapping parameters are not specified by the user, then this function estimates the parameters using
-#' a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. The optimisation by default
-#' includes the variogram parameters (e.g. range, sill and nugget) and the search parameters for local kriging (e.g. radius, minimum and
-#' maximum number of observations to use). Optimising these parameters is not common in kriging. It is done herein because trials for Victoria,
-#' Australia, showed that calibrating these parameters produced significantly lower cross-validation errors (i.e. the error in predicting the observations
-#' removed from the optimisation) compared to the standard approach of graphical estimation from an experimental variogram. The optimisation is
-#' numerically challenging and the following factors should be considered before use:
-#'
-#' \itemize{
-#'  \item{Optimisation of the parameters \code{mrvbf.pslope}, \code{mrvbf.ppctl} and \code{smooth.std} often required the creating of raster grids for every parameter combination. To ease the computation burden, these parameters should be treated as discrete, not continuous, numbers.}
-#'  \item{The optimisation package \code{rgeoud} is used herein. For control the optimisation process, consider directly using \code{\link{krige.head.calib}}.}
-#'  \item{Trials have established default calibration parameters and settings that were effective for Victoria, Australia. There is no guarantee they will be effective for other regions.}
-#' }
+#' a mixed data-type (i.e. real and integer parameters) split-sample maximum likelihood global optimisation. See \code{\link{krige.head.calib}} for details.
+#' 
+#' However, it is advised to separately calibrate the mapping parameters using \code{\link{krige.head.calib}}. The output from 
+#' the calibration can then be passed into \code{krige.head} via the input \code{calibration.results}. This will use the calibrated setting 
+#' and the training observed data set to map the head. To use all of the observed data, or alternate data, input \code{data}. 
+#' See \code{\link{krige.head.calib}} Examples for details.
 #'
 #' In using this function, the primary user decisions are:
 #' \itemize{
@@ -30,6 +24,9 @@
 #'  \item{The type of variogram model, defined by the input \code{model}}
 #'  }
 #'
+#' @param \code{calibration.results} List variable returned from \code{\link{krige.head.calib}}. Inputting this applies the calibration results here. All other inputs are ignores, except \code{data}, \code{newdata} and \code{grid}. 
+#' Defauli is \code{NULL}.
+#' 
 #' @param \code{formula} defines the R formula (as a character or formula data type) to be used to interpolate the heads. The left hand side of the formula must be \code{head}.
 #' The right hand side can contain any or all of the following terms: \code{elev} for the land surface elevation;
 #' \code{MrVBF} for the Multiresolution Index of Valley Bottom Flatness as a measure of valley-ness at each DEM grd cell;
@@ -165,17 +162,10 @@
 #' 
 #' # Export grids to an ARCMAP *.asc file
 #' write.asciigrid(heads,'DBNS.asc',3,na.value = -999)
-#' 
-#' # Re-calibrate the parameters and map using the default settings.
-#' # Note, here only those parameters that are a vector are calibrated. This is in addition
-#' # to the variogram parameters. 
-#' varigram.model = vgm(psill=10, model='Mat', range= 10000 , nugget=1, kappa=0.1);
-#' heads <- krige.head(formula=f, grid=DEM, data=obs.data, data.errvar.colname='total_err_var',
-#' nmax=Inf, maxdist=Inf, smooth.std = seq(0.5,1.5,length.out=11), 
-#' model = 'Mat',  fit.variogram.type=1, debug.level=1)
 #'
 #' @export
 krige.head <- function(
+  calibration.results = NULL,  
   formula = as.formula("head ~ elev + MrVBF + MrRTF + smoothing"),
   grid=NULL,
   grid.landtype.colname = NULL,
@@ -231,6 +221,32 @@ krige.head <- function(
   if (!do.grid.est && nsim!=0)
     stop('Point estimation cannot be undertaken for geostatistical simulations. To undertake gridded simulations, set newdata=NULL and nsim>0.')
 
+  # Use calibration results if provided. Else use individual input.
+  if (!is.null(calibration.results)) {
+
+    formula = calibration.results$inputs$formula
+    if (is.null(grid))
+      grid = calibration.results$inputs$grid
+    grid.landtype.colname = calibration.results$inputs$grid.landtype.colname;
+    if (is.null(data))
+      data = calibration.results$inputs$data
+    data.fixedHead = calibration.results$inputs$data.fixedHead
+    if (is.null(newdata))
+      newdata = calibration.results$inputs$newdata
+    data.errvar.colname = calibration.results$inputsdata.errvar.colname
+    model= calibration.results$variogramModel
+    fit.variogram.type = 3; 
+
+    mrvbf.pslope = calibration.results$parameter.Values$mrvbf.pslope
+    mrvbf.ppctl = calibration.results$parameter.Values$mrvbf.ppctl
+    smooth.std = calibration.results$parameter.Values$smooth.std
+    nmax = calibration.results$parameter.Values$nmax
+    nmax.fixedHead = calibration.results$parameter.Values$nmax.fixedHead
+    maxdist = calibration.results$parameter.Values$maxdist
+    nmin = calibration.results$parameter.Values$nmin 
+    omax = calibration.results$parameter.Values$omax
+  }
+
   # Get variable names and check if using MrVBF or MrRTF or DEM
   if (debug.level>0)
     message('... Checking formula terms.');

.Rproj.user/19E94E50/sources/per/t/F6CA6EB6-contents → .Rproj.user/19E94E50/sources/session-4dda64d8/5B688876-contents

@@ -121,140 +121,25 @@ get.objFunc <- function(
       return(Inf)
     } 
   }
-
-  # Build variogram model
-  if (fit.variogram.type == 1) {
 
-    # Check the format of the model input.
-    if (class(model)[1] == "variogramModel") {
-
-      # Get model types
-      variogram.type = as.character(model$model)
-
-      # Remove 'nug' from list of model types and add to 'model' variable
-      filt = variogram.type != 'Nug'      
-      modelType = variogram.type[filt]            
-    } else if (is.character(model)) {
-      modelType = model    
-    } else {
-      stop('When fit.variogram.type==1, the model should be a gstat variogram model or a vector of model types, eg c(\'Mat\')')
-    }
-
-    # Initialise variogram settings
-    variogram.nugg = NULL
-
-    variogram.psill = c()
-
-    variogram.range = c()
-
-    variogram.kappa = c()
-
-    variogram.ang1 = c()
-
-    variogram.anis1 = c()
+  # Build variogram model and turn off future fitting via get.variogram()
+  if (fit.variogram.type == 1) {
+    model = get.variogramModel(params, param.names, model)
+    fit.variogram.type = 3
 
-    nModel = 0
 
+    # Get variogram parameter names and values - for reporting
     # Extract variogram settings.
     for (i in 1:length(param.names)) {
-
       if (!is.na(pmatch('psill', param.names[i])) ||
           !is.na(pmatch('range', param.names[i])) ||
           !is.na(pmatch('kappa', param.names[i])) ||
           !is.na(pmatch('nug', param.names[i]))  ||
           !is.na(pmatch('ang1', param.names[i])) ||
           !is.na(pmatch('anis1', param.names[i]))) {
-
-        if (is.na(pmatch('nug', param.names[i]))) {
-          iModel = as.numeric(substr(param.names[i], nchar(param.names[i]), nchar(param.names[i])))
-
-          # Initialise some extra variogram parameters in case they're not used.        
-          if (iModel>nModel) {
-            variogram.kappa[iModel] = 0.5
-            variogram.ang1[iModel] = 0.0
-            variogram.anis1[iModel] = 1.0
-          }
-
-        }
-
-        if (!is.na(pmatch('nug', param.names[i]))) {
-          variogram.nug = max(0, params[i])
-
-        } else if (!is.na(pmatch('psill', param.names[i]))) {
-          variogram.psill[iModel] = max(0, params[i])
-
-          if (iModel > nModel)
-            nModel = iModel
-
-        } else if (!is.na(pmatch('range', param.names[i]))) {
-          variogram.range[iModel] = max(1e-6, params[i])
-
-          if (iModel > nModel)
-            nModel = iModel
-
-        } else if (!is.na(pmatch('kappa', param.names[i]))) {
-          variogram.kappa[iModel] = max(0, params[i])
-
-          if (iModel > nModel)
-            nModel = iModel
-
-        } else if (!is.na(pmatch('ang1', param.names[i]))) {
-          variogram.ang1[iModel] = params[i]
-
-          if (iModel > nModel)
-            nModel = iModel
-
-        } else if (!is.na(pmatch('anis1', param.names[i]))) {
-          variogram.anis1[iModel] = params[i]
-
-          if (iModel > nModel)
-            nModel = iModel
-
-        }
-
         params2Print[i] = params[i]
       }
     }
-
-    # Build variogram
-    for (i in 1:nModel) {
-
-      #message(paste('... DEBGUGGING: i, model type, nug, psill, range, kappa ',i, model[i], variogram.nug, variogram.psill[i], variogram.range[i], variogram.kappa[i]))
-
-      if (i == 1 && !is.null(variogram.nug)) {
-        vgm.model = vgm(
-          nugget = variogram.nug,
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      } else if (i == 1) {
-        vgm.model = vgm(
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      } else {
-        vgm.model = vgm(
-          psill = variogram.psill[i],
-          range = variogram.range[i],
-          model = modelType[i],
-          add.to = vgm.model,
-          kappa = variogram.kappa[i],
-          anis = c(variogram.ang1, variogram.anis1)
-        )
-      }
-    }
-    model = vgm.model
-
-
-    # Turn off fitting within get.variogram()
-    fit.variogram.type = 3
-
   }
 
   # Call krige.head with the input parameter names.