build_scripts/make_dynatop.R

In waternumbers/dynatop: An Implementation of Dynamic TOPMODEL Hydrological Model in R

## Script for creating and building the dynatop R package
rm(list=ls())
graphics.off()

## path of the package
pacPath <- '.'
Rcpp::compileAttributes(pacPath)
devtools::document(pacPath)
devtools::spell_check(pacPath,use_wordlist=TRUE)
devtools::check(pacPath)

## check documentation build
pkgdown::clean_site(pacPath)
pkgdown::build_site(pacPath)
pkgdown::clean_site(pacPath)

## build, populate drat
## linux
dratPath <- "~/Documents/Software/drat"
buildFile <- devtools::build(pacPath)
install.packages(buildFile)
drat::insertPackage(buildFile,dratPath)

## mac and windows
rhub::validate_email() # for first time that session
pkgName <- sub('\\.tar.gz$', '', basename(buildFile))
## rhub::platforms()[,1] # lists platforms

mch <- rhub::check(path = buildFile,
                   platform = c("macos-highsierra-release-cran","windows-x86_64-release",
                                "macos-m1-bigsur-release"))
ext <- c(".tgz",".zip",".tgz")
outPath <- dirname(buildFile)
for(ii in 1:2){ ## m1 not fixed yet in drat
    tmp <- paste0(pkgName,ext[ii])
    outFile <- file.path(outPath,tmp)
    #download.file(file.path(mch$urls()$artifacts[ii],tmp),outFile)
    drat::insertPackage(outFile,dratPath)
}

## prior to submission
## submit to win-builder (release and devel) and macbuilder
## run locally with R CMD check --as-cran --use-valgrind ...
## run locally with r-devel and ASAN build

###########################################################################
## This converts the data for Swindale into the data object used in the examples.
rm(list=ls())
pacPath <- '../'
devtools::load_all(pacPath)
model <- readRDS("Swindale_exp.rds")

qr <- read.csv( "start=2009-11-18_end=2009-11_4_int=0.25-hours_units=mm.hr-1.tsv",sep="\t")
obs <- as.xts(qr[,c("Flow","Rainfall")],order.by= as.POSIXct(qr[,'Date'],format="%d/%m/%Y %H:%M",tz='GMT'))
## According to original notes and code Flow in cumecs and precip in mm/timestep
obs$Rainfall <- obs$Rainfall/1000 # convert to m/timestep
obs$PET <- evap_est(index(obs),0,5/1000) # in m
names(obs) <- c("flow","precip","pet")
Swindale <- list(model=model,obs=obs)
save("Swindale",file="../data/Swindale.rda")






## ########################################
## This code uses Swindale and calls all the different transmissivity profiles
## with models generated by dynatopGIS
rm(list=ls())
devtools::load_all("../")
data("Swindale");
for(ii in list.files(".",pattern="^Swindale.*\\.rds$")){
    ##assign("last.warning", NULL, envir = baseenv())
    print(ii)
    mdl <- readRDS(ii)
    mdl$precip_input$name <- "Rainfall"
    mdl$pet_input$name <- "PET"
    dt <- dynatop$new(mdl)$add_data(Swindale$obs)#[1:2,,drop=FALSE])
    ##tryCatch({
        dt$initialise()$sim_hillslope()
        print(range(dt$get_mass_errors()[,6]))
    ##},
    ##warning = function(w){print(warnings())},
    ##error = function(e){stop(e)})

    print(warnings())

}

## ##############################
## simple run for debugging
rm(list=ls())
devtools::load_all("../")
data("Swindale");

mdl <- Swindale$model

## temp fix to get correct model form
odfn <- data.frame(name="outlet",id=0,flux="q_sf")
h <- list()
for(ii in 1:nrow(mdl$hru)){
    tmp <- list()
    tmp$id <- as.integer( mdl$hru$id[ii]-1 )
    tmp$states <- setNames(as.numeric(rep(NA,4)), c("s_sf","s_rz","s_uz","s_sz"))
    tmp$properties <- c(width =  mdl$hru$area[ii]/mdl$hru$length[ii], area = mdl$hru$area[ii], gradient = mdl$hru$s_bar[ii])

    ## mdl$hru$width[ii], area = mdl$hru$area[ii], gradient = mdl$hru$s_bar[ii])
    ##if( tmp$properties["width"] ==0 ){ tmp$properties["width"] <- 1 }
    tmp$sf <- list(type = mdl$hru$sf[[ii]]$type,
                   parameters = c("c_sf" = 0.8))
    tmp$sz <- list(type = "bexp",
                   parameters = c(t_0=exp(7.46),m=0.0063,D=0.5)) #c(mdl$hru$sz[[ii]]$param, "D" = 0.05))
    if(tmp$id==0){ tmp$sz$parameters["t_0"] <- 0 }
    tmp$precip <- mdl$hru$precip[[ii]]
    names(tmp$precip) <- c("name","fraction")
    tmp$pet <- mdl$hru$pet[[ii]]
    names(tmp$pet) <- c("name","fraction")
    tmp$rz <- list(type="orig", parameters = c("s_rzmax" = 0.1))#mdl$hru$s_rzmax[ii]))
    tmp$uz <- list(type="orig", parameters = c("t_d" = 8*60*60)) #mdl$hru$t_d[ii]))
    tmp$sf_flow_direction <- list(id = as.integer(mdl$hru$sf_flow_direction[[ii]]$id - 1),
                                  fraction = mdl$hru$sf_flow_direction[[ii]]$frc)
    tmp$sz_flow_direction <- list(id = as.integer(mdl$hru$sz_flow_direction[[ii]]$id - 1),
                                  fraction = mdl$hru$sz_flow_direction[[ii]]$frc)
    tmp$initialisation <- c("s_rz_0" = 0.98, #mdl$hru$s_rz0[ii], "r_uz_sz_0" = mdl$hru$r_uz_sz0[ii])
                            "r_uz_sz_0" = 1.75e-20) #mdl$hru$r_uz_sz0[ii])
    h[[ii]] <- tmp
}

hh <- h
hh[[1]]$id <- hh[[1]]$id + 0
system.time({
    dt <- dynatop$new(h)

    dt$add_data(Swindale$obs)
    dt$initialise()

    dt$sim(odfn)
})

plot(dt$get_output())
plot(Swindale$obs$flow)