update.

SHUD-System · Mar 14, 2023 · b459cfd · b459cfd
1 parent 8028c54
commit b459cfd
Show file tree

Hide file tree

Showing 34 changed files with 284 additions and 182 deletions.
diff --git a/.gitignore b/.gitignore
@@ -8,3 +8,5 @@ src/*.dll
 experiments
 ^.DS_Store
 *.swp
+.DS_Store
+Rplots.pdf
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -32,16 +32,18 @@ Imports: Rcpp,
   abind,
   utils,
   lubridate,
-  interp,
-  geometry,
   doParallel,
+  geometry,
   methods
 LinkingTo: Rcpp, sp
-RoxygenNote: 7.1.2
+RoxygenNote: 7.2.1
 Remotes: shulele/RTriangle/pkg
 Suggests: testthat,
     knitr,
     rmarkdown,
+    terra,
   deldir,
+  interp,
+  whitebox,
   ncdf4
 VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
@@ -37,6 +37,8 @@ export(RiverAtt)
 export(RiverLength)
 export(RiverSlope)
 export(RiverType)
+export(SHUD.MESH)
+export(SHUD.RIVER)
 export(SharedPoints)
 export(SimpleSpatial)
 export(SimplifybyLen)
@@ -50,14 +52,15 @@ export(autoBuildModel)
 export(compareMaps)
 export(correctRiverSlope)
 export(count)
+export(crs.Albers)
+export(crs.Lambert)
 export(crs.long2utm)
 export(crs.long2utmZone)
 export(datafilter.riv)
 export(days_in_year)
 export(extractCoords)
 export(extractRaster)
 export(fdc)
-export(filebackup)
 export(fishnet)
 export(getAquiferDepth)
 export(getArea)
@@ -98,6 +101,7 @@ export(readmesh)
 export(readnc)
 export(readnc.time)
 export(readout)
+export(readout.header)
 export(readpara)
 export(readriv)
 export(readriv.sp)
@@ -135,6 +139,7 @@ export(ts2df)
 export(ts2map)
 export(tsLaiRf)
 export(voronoipolygons)
+export(watershedDelineation)
 export(wb.DS)
 export(wb.all)
 export(wb.ele)

diff --git a/R/Func_landcover.R b/R/Func_landcover.R
@@ -341,6 +341,7 @@ lc.UMD <- function(){
   cn = toupper(c('INDEX', 'LAIMAX', 'RMIN', 'RSREF', 'ALBEDO',
                  'VEGFRAC', 'ROUGH', 'RZD', 'SOILDGRD', 'IMPAF') )
   colnames(dtab) = cn
+  dtab = dtab[, -1 * (2:4)]
   return(dtab)
 }
 
@@ -353,7 +354,7 @@ lc.NLCD <- function(lc){
   dtab = lc.UMD()
   y= t(lc_EQ( t(dtab), lc))
   y[,1] = 1:nrow(y)
-  colnames(y) = toupper(c('INDEX', 'LAIMAX', 'RMIN', 'RSREF', 'ALBEDO',
+  colnames(y) = toupper(c('INDEX', 'ALBEDO',
                           'VEGFRAC', 'ROUGH', 'RZD', 'SOILDGRD', 'IMPAF') )
   rownames(y)=lc
   return(y)
@@ -364,9 +365,7 @@ lc.NLCD <- function(lc){
 #' @return Default land cover parameters
 #' @export
 lc.GLC <- function(){
-  cn = toupper(c('INDEX', 'LAIMAX', 'RMIN', 'RSREF', 'ALBEDO',
-                 'VEGFRAC', 'ROUGH', 'RZD', 'SOILDGRD', 'IMPAF') )
-
+  cn = toupper(c('INDEX', 'ALBEDO', 'VEGFRAC', 'ROUGH', 'RZD', 'SOILDGRD', 'IMPAF') )
   mapid = c(1:11, 1, 12, 14, 12, 1, 13)
   v = lc.UMD()[mapid, ]
   v[, 1] = 1:nrow(v)-1

diff --git a/R/MeshDomain.R b/R/MeshDomain.R
@@ -121,3 +121,11 @@ shud.att <- function(tri, r.soil =NULL, r.geol=NULL, r.lc=NULL, r.forc=NULL,
   att$LAKE = apply.raster(sp.lake, p.centroids, v0=0)
   return(att)
 }
+
+# 
+# pa=shud.att(tri, 
+#             r.soil = r.soil, r.geol = r.geol, 
+#             r.lc = rlc.idx, 
+#             r.forc = sp.forc, 
+#             r.BC = 0, 
+#             sp.lake = sp.lake)
diff --git a/R/ModelConfigure.R b/R/ModelConfigure.R
@@ -31,39 +31,51 @@ shud.ic <- function(ncell, nriv, AqD = 10, stage = 0.1, p1=0.4, p2 = p1,
 #' @return model configureation, Vector
 #' @export
 shud.para <- function( nday = 10){
-  dts = c(
-    paste0('dt_',
-           c(paste0('ye_', c('snow', 'surf', 'unsat', 'gw') ), 
-             paste0('Qe_', c('surf', 'sub') ),
-             paste0('qe_et'), 
-             paste0('qe_', c('prcp', 'infil', 'rech') ),
-             paste0('yr_', c('stage')),
-             paste0('Qr_', c('down', 'surf', 'sub', 'up')),
-             paste0('lake')
-           ))
-  )
-
-  vdt = rep(1440, length(dts))
-
+  dts = c(paste0('dt_',
+                 c(paste0('ye_', c('snow', 'surf', 'unsat', 'gw') ), 
+                   paste0('Qe_', c('surf', 'sub') ),
+                   paste0('qe_et'), 
+                   paste0('qe_', c('prcp', 'infil', 'rech') ),
+                   paste0('yr_', c('stage')),
+                   paste0('Qr_', c('down', 'surf', 'sub', 'up')),
+                   paste0('lake')
+                 ))  )
+  vdt = rep(0, length(dts))
   vn = c('VERBOSE', 'INIT_MODE', 'CloseBoundary', 
          'ASCII_OUTPUT', 'Binary_OUTPUT', 'cryosphere',
          'SpinupDay', 'NUM_OPENMP', 'SCR_INTV',
          'ABSTOL', 'RELTOL', 
          'INIT_SOLVER_STEP', 'MAX_SOLVER_STEP', 'LSM_STEP', 
-         'START', 'END', 
-         dts)
+         'START', 'END',  dts)
   val = c(0,  3, 1,
-        0, 1, 0,
-        0, 8, 1440,
-        1e-4, 1e-4,
-        1e-2, 2, 60,
-        0, nday,
-        vdt
-  )
+          0, 1, 0,
+          0, 8, 1440,
+          1e-4, 1e-4,
+          1e-2, 2, 60,
+          0, nday,
+          vdt )
   val=data.frame(rbind(val))
-  names(val) = toupper( vn )
-  val
+  names(val) = toupper( vn )    
+  val$DT_YE_SNOW = 0
+  val$DT_YE_SURF = 0
+  val$DT_YE_UNSAT = 0
+  val$DT_YE_GW = 1440
+  val$DT_QE_SURF = 0
+  val$DT_QE_SUB = 0
+  val$DT_QE_ET = 1440
+  val$DT_QE_PRCP = 1440
+  val$DT_QE_INFIL = 0
+  val$DT_QE_RECH = 0
+  val$DT_YR_STAGE = 0
+  val$DT_QR_DOWN = 1440
+  val$DT_QR_SURF = 0
+  val$DT_QR_SUB = 0
+  val$DT_QR_UP = 0
+  val$DT_LAKE = 1440
+
+  return(val)
 }
+
 #' Generate the default model calibration
 #' \code{shud.calib} 
 #' @return Default calibration values for model

diff --git a/R/NetCDF.R b/R/NetCDF.R
@@ -120,11 +120,12 @@ readnc<-function(ncid, varid=NULL,  ext = NULL){
 #' @param x  list(x, y, arr) or x coordinates
 #' @param y y coordinates
 #' @param arr array. dim = (nx, ny)
+#' @param Dxy Resolution
 #' @param flip Whether flip the matrix.
 #' @param plot Whether plot the raster.
 #' @return a Raster/RasterStack
 #' @export
-xyz2Raster <- function(x, y=NULL, arr=NULL,
+xyz2Raster <- function(x, y=NULL, arr=NULL,Dxy=NULL,
                        flip=TRUE, plot=TRUE){
   if(is.null(y) | is.null(arr)){
     nc = x;
@@ -134,22 +135,33 @@ xyz2Raster <- function(x, y=NULL, arr=NULL,
   dims = dim(nc$arr)
   ndims = length(dims)
   x = nc$x; y = nc$y
+
+  if(is.null(Dxy) & (length(x) == 1 | length(y) == 1) ){
+    message('Dxy cannot be NULL. Dxy MUST be given.')
+    stop()
+  }
+
   if( ndims > 2){
     # multiple layers
     rl = list()
     for(i in 1:dims[3]){
-      rl[[i]] = xyz2Raster(x = x, y = y, 
+      rl[[i]] = xyz2Raster(x = x, y = y, Dxy=Dxy,
                            arr=matrix(nc$arr[, , i], nrow=dim(nc$arr)[1], ncol=dim(nc$arr)[2]) )
     }
     rs = raster::stack(rl)
   }else{
     # single layer
     # val = matrix(arr, nrow=nrow(arr), ncol=ncol(arr))
     val = arr
-    dx = abs(mean(diff(x))); 
-    if(length(y) ==1) {
-      dy=dx
+    if(!is.null(Dxy)){
+      dx = Dxy[1]
+      if(length(Dxy)==1){
+        dy = Dxy[1]
+      }else{
+        dy = Dxy[2]
+      }
     }else{
+      dx = abs(mean(diff(x)));
       dy = abs(mean(diff(y)))
     }
     nx = length(x);   ny = length(y)

diff --git a/R/autoBuildModel.R b/R/autoBuildModel.R
@@ -15,7 +15,6 @@
 #' @param clean Whether clean the existing model files in output directory.
 #' @param cfg.para  model configuration, parameter
 #' @param cfg.calib model calibration
-#' @param backup Whether to write backup files.
 #' @param rm.outlier Whether to remove the outlier in soil/geol;
 #' @param mf Meltfactor
 #' @param quiet Whether to ask confirmation when multiple outlets exist.
@@ -50,8 +49,7 @@ autoBuildModel <- function(
   cfg.calib = shud.calib(),
   mf = MeltFactor(years = years),
   rm.outlier = TRUE, 
-  quiet = FALSE,
-  backup=TRUE
+  quiet = FALSE
 ){
   msg = paste0('autoBuildModel(', prjname, '):: ')
   message(msg, 'Automatic build the SHUD model.')
@@ -152,7 +150,7 @@ autoBuildModel <- function(
   }
   pa=shud.att(tri, r.soil = xsoil, r.geol = xgeol, r.lc = rlc, r.forc = sp.forc )
 
-  write.forc(forcfiles,  file=fin['md.forc'], startdate = paste0(min(years), '0101'), backup = backup)
+  write.forc(forcfiles,  file=fin['md.forc'], startdate = paste0(min(years), '0101'))
 
   # generate SHUD .riv
   AA = rgeos::gArea(wbd) * 1e-6
@@ -255,24 +253,24 @@ autoBuildModel <- function(
   graphics::legend('top', paste0(lc), col=col, lwd=1)
   dev.off()
   message(msg, 'Write SHUD model input files.')
-  write.tsd(backup=backup,lr$LAI, file = fin['md.lai'])
-  write.tsd(backup=backup,lr$RL, file = fin['md.rl'])
+  write.tsd(lr$LAI, file = fin['md.lai'])
+  write.tsd(lr$RL, file = fin['md.rl'])
 
-  write.tsd(backup=backup,mf, file=fin['md.mf'])
+  write.tsd(mf, file=fin['md.mf'])
 
   # write SHUD input files.
-  write.mesh(backup=backup,pm, file = fin['md.mesh'])
-  write.riv(backup=backup,pr, file=fin['md.riv'])
-  write.ic(backup=backup,pic, file=fin['md.ic'])
+  write.mesh(pm, file = fin['md.mesh'])
+  write.riv(pr, file=fin['md.riv'])
+  write.ic(pic, file=fin['md.ic'])
 
-  write.df(backup=backup,pa, file=fin['md.att'])
-  write.df(backup=backup,prs, file=fin['md.rivseg'])
-  write.df(backup=backup,para.lc, file=fin['md.lc'])
-  write.df(backup=backup,para.soil, file=fin['md.soil'])
-  write.df(backup=backup,para.geol, file=fin['md.geol'])
+  write.df(pa, file=fin['md.att'])
+  write.df(prs, file=fin['md.rivseg'])
+  write.df(para.lc, file=fin['md.lc'])
+  write.df(para.soil, file=fin['md.soil'])
+  write.df(para.geol, file=fin['md.geol'])
 
-  write.config(backup=backup,cfg.para, fin['md.para'])
-  write.config(backup=backup,cfg.calib, fin['md.calib'])
+  write.config(cfg.para, fin['md.para'])
+  write.config(cfg.calib, fin['md.calib'])
   message(msg, 
           'Ncell=', nrow(pm@mesh), '\t',
           'Nriv=', nrow(pr@river), '\t',
@@ -291,5 +289,4 @@ autoBuildModel <- function(
 #                   tol.len = 0,
 #                   tol.riv = 0,
 #                   tol.wb = sqrt(a.max)/3,
-#                   quiet = TRUE,
-#                   backup=FALSE)
+#                   quiet = TRUE)
diff --git a/R/readinput.R b/R/readinput.R
@@ -2,9 +2,10 @@
 #' \code{read.df} 
 #' @param file full path of file
 #' @param text text return from readLines(file)
+#' @param sep seperator of the table.
 #' @return a list of matrix
 #' @export
-read.df <-function(file, text = readLines(file)){
+read.df <-function(file, text = readLines(file), sep='\t'){
   # fn=fin['md.mesh']
   # text=readLines(fn)
   r0 = 1
@@ -18,7 +19,8 @@ read.df <-function(file, text = readLines(file)){
       nr = nrow-2
     }else{
     }
-    xl[[i]] = utils::read.table(text = text[0:nr + 1 + r0], header = T)
+    # print(text[0:nr + 1 + r0])
+    xl[[i]] = utils::read.table(text = text[0:nr + 1 + r0], header = T, sep = sep)
     # xl[[i]] = as.matrix(utils::read.table(text = text[0:nr + 1 + r0], header = T))
     r0 = r0 + nr + 2;
     if(r0 + 1 > nrow){
@@ -27,6 +29,7 @@ read.df <-function(file, text = readLines(file)){
   }
   xl
 }
+
 #' Read the .mesh file
 #' \code{readmesh} 
 #' @param file full path of .mesh file