Replace deprecated alloc / free + modern repo structure

README.md

@@ -0,0 +1,56 @@
+# TOPMODEL R Package
+
+An R implementation of the hydrological model TOPMODEL, based on the 1995 FORTRAN version by Keith Beven. This package provides a set of hydrological functions for rainfall-runoff modeling and catchment analysis.
+
+## Installation
+
+This package is currently not available on CRAN as it is in maintenance mode. You can install it from GitHub using:
+
+```r
+devtools::install_github("ICHydro/topmodel")
+```
+
+## Main Features
+
+- **Rainfall-runoff modeling**: Simulate discharge from precipitation and evapotranspiration data
+- **Topographical analysis**: Calculate topographic indices and flow delay functions from digital elevation models
+- **Sensitivity analysis**: Explore parameter sensitivity using Monte Carlo sampling
+- **Uncertainty analysis**: GLUE (Generalized Likelihood Uncertainty Estimation) framework for prediction uncertainty
+
+## Quick Example
+
+```r
+library(topmodel)
+
+# Load example data into global environment
+data(huagrahuma)
+list2env(huagrahuma, envir = .GlobalEnv)
+
+# Run the model
+Qsim <- topmodel(
+  parameters,
+  topidx,
+  delay,
+  rain,
+  ETp
+)
+
+# Evaluate performance
+NSeff(
+  Qobs,
+  Qsim
+)
+```
+
+## Background
+
+TOPMODEL is a physically-based, variable contributing area model of basin hydrology that uses topographic indices to represent the spatial variability of hydrological processes. The model was originally developed by Beven and Kirkby (1979) and has been widely used in hydrological research and applications.
+
+## References
+
+- Beven, K. J., Kirkby, M. J. (1979). A physically based variable contributing area model of basin hydrology. *Hydrological Sciences Bulletin*, 24, 43-69.
+- Beven K, Lamb R, Quinn P, Romanowicz R, Freer J (1995). TOPMODEL. In: Singh VP (Ed), *Computer Models of Watershed Hydrology*. Water Resources Publications, Colorado. pp. 627-668.
+
+## Documentation
+
+For detailed examples and documentation, see the package help files. A complete example workflow is available in `inst/examples/Full Run.R`.

examples/Full Run.R → inst/examples/Full Run.R

@@ -1,4 +1,3 @@
-
 ## install and load the required packages:
 
 install.packages("topmodel")
@@ -17,11 +16,11 @@ library(Hmisc)
 # Take for instance this minimalistic DEM, saved in a test file called "DEM.txt"
 # Values outside the catchment are given the value -9999 (this can be any other value that, obviously, does not occur in the DEM values):
 
--9999 -9999 828.9 835.6 -9999
-818.3 826.0 830.7 834.5 836.0
-817.1 824.0 825.2 833.3 836.9
-816.5 820.0 824.1 330.8 -9999
-810.7 815.6 822.2 -9999 -9999
+# -9999 -9999 828.9 835.6 -9999
+# 818.3 826.0 830.7 834.5 836.0
+# 817.1 824.0 825.2 833.3 836.9
+# 816.5 820.0 824.1 330.8 -9999
+# 810.7 815.6 822.2 -9999 -9999
 
 # This file can be imported and processed in R with:
 
@@ -30,7 +29,7 @@ DEM <- as.matrix(DEM)
 
 # Remove the values outside the catchment:
 
-DEM[DEM==-9999] <- NA
+DEM[DEM == -9999] <- NA
 
 # You may want to plot the DEM to see whether everything looks OK:
 
@@ -40,21 +39,21 @@ image(DEM)
 # Here we use the DEM from Huagrahuma as an example:
 
 data(huagrahuma.dem)
-DEM <- sinkfill(huagrahuma.dem, cellsize=25, degree=0.1)
-topindex <- topidx(DEM, resolution=25)
+DEM <- sinkfill(huagrahuma.dem, cellsize = 25, degree = 0.1)
+topindex <- topidx(DEM, resolution = 25)
 
 # The values need to be split into a set of classes, since topmodel() is a semidistributed model that lumps hydrologically similar areas into the same hydrological response units.
 # Here we define 16 hydrological response units:
 
-topidx <- make.classes(topindex,16)
+topidx <- make.classes(topindex, 16)
 
 # the delay function is a bit more tricky because this requires cumulative fractions, but you generate it as follows:
 
 n <- 5 # number of classes; a higher number will result in a smoother histogram
-delay <- flowlength(huagrahuma.dem)*25 # TODO: add the outlet coordinates; otherwise the flowlength will be calculated to the edge of the map.
+delay <- flowlength(huagrahuma.dem) * 25 # TODO: add the outlet coordinates; otherwise the flowlength will be calculated to the edge of the map.
 delay <- make.classes(delay, n)
-delay <- delay[n:1,]
-delay[,2] <- c(0, cumsum(delay[1:(n-1),2]))
+delay <- delay[n:1, ]
+delay[, 2] <- c(0, cumsum(delay[1:(n - 1), 2]))
 
 ############ PART 1: running the rainfall-runoff model ##############
 
@@ -71,12 +70,12 @@ topidx
 parameters
 rain
 
-plot(rain, type="h")
+plot(rain, type = "h")
 
 ## run the model and visualise the outcome:
 
 Qsim <- topmodel(parameters, topidx, delay, rain, ET0)
-plot(Qsim, type="l", col="red")
+plot(Qsim, type = "l", col = "red")
 points(Qobs)
 
 ## Evaluate the model with a performance metric
@@ -99,7 +98,7 @@ NSeff(Qobs, Qsim)
 ## What value do you get? Do you think this is a good simulation?
 ## Verify by plotting:
 
-plot(Qsim, type="l", col="red")
+plot(Qsim, type = "l", col = "red")
 points(Qobs)
 
 ## Now sample all parameters at random. We take a sample size of 100
@@ -118,7 +117,7 @@ k0 <- runif(n, min = 0, max = 10)
 CD <- runif(n, min = 0, max = 5)
 dt <- 0.25
 
-parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt)
+parameters <- cbind(qs0, lnTe, m, Sr0, Srmax, td, vch, vr, k0, CD, dt)
 
 ## run the model and evaluate with the Nash – Sutcliffe efficiency metric:
 ## Note: the function accepts a table of parameter sets
@@ -129,22 +128,22 @@ max(NS)
 
 ## visualisation of the sensitivity using dotty plots:
 
-plot(lnTe, NS, ylim = c(0,1))
+plot(lnTe, NS, ylim = c(0, 1))
 
 ############ PART 3: GLUE uncertainty analysis ##############
 
 ## choose a behavioural threshold and remove the “bad” parameter sets:
 
-parameters <- parameters[NS > 0.3,]
+parameters <- parameters[NS > 0.3, ]
 NS <- NS[NS > 0.3]
 
 ## generate predictions for the behavioural parameter sets:
 
-Qsim <- topmodel(parameters,topidx,delay,rain,ET0)
+Qsim <- topmodel(parameters, topidx, delay, rain, ET0)
 
 ## (have a look at the predictions for the first time step:)
 
-hist(Qsim[1,])
+hist(Qsim[1, ])
 
 ## construct weights based on the performance measure
 
@@ -154,19 +153,20 @@ weights <- weights / sum(weights)
 ## make prediction boundaries by weighted quantile calculation
 ## (we need the Hmisc package for that)
 
-limits <- apply(Qsim, 1, "wtd.quantile", weights = weights,
-probs = c(0.05,0.95), normwt=T)
+limits <- apply(Qsim, 1, "wtd.quantile",
+  weights = weights,
+  probs = c(0.05, 0.95), normwt = T
+)
 
-plot(limits[2,], type="l")
-points(limits[1,], type="l")
-points(Qobs, col="red")
+plot(limits[2, ], type = "l")
+points(limits[1, ], type = "l")
+points(Qobs, col = "red")
 
 ## how many measurements fall outside of the prediction limits?
 
-outside <- (Qobs > limits[2,]) | (Qobs < limits[1,])
+outside <- (Qobs > limits[2, ]) | (Qobs < limits[1, ])
 summary(outside)
 
 ## width of the prediction boundaries
 
-mean(limits[2,] - limits[1,]) / mean(Qobs, na.rm=T)
-
+mean(limits[2, ] - limits[1, ]) / mean(Qobs, na.rm = T)

topmodel/src/c_topidx.c → src/c_topidx.c

@@ -10,7 +10,7 @@
 
 void c_topidx(double *inputdem,
 	    int    *inputriver,
-	    int    *nrow, 
+	    int    *nrow,
 	    int    *ncol,
 	    double *ew_res,
 	    double *ns_res,
@@ -21,12 +21,12 @@ void c_topidx(double *inputdem,
   int    nrout,river,not_yet;
   double **dem, **atb, **area, **slope, **rivermap;
   double exclude,dnx,routefac,nslp;
-  int    nsink = 0;
+  // int    nsink = 0;
   double routdem[9], tanb[9];
   double c,dx1,dx2,sum,sumtb;
 
   /* memory allocation */
-          
+
   dem   = (double **) R_alloc(*nrow, sizeof(double *));
   atb   = (double **) R_alloc(*nrow, sizeof(double *));
   area  = (double **) R_alloc(*nrow, sizeof(double *));
@@ -97,7 +97,7 @@ void c_topidx(double *inputdem,
 
     for(j = 0; j < *ncol; j++) {
       for(i = 0; i < *nrow; i++) {
-              
+
         /* skip non catchment cells and cells that are done */
         if((dem[i][j] == exclude) || (atb[i][j] >= ZERO))
 	  continue;
@@ -108,15 +108,15 @@ void c_topidx(double *inputdem,
         if(rivermap[i][j] == 1) river = 1;
         else {
 
-          /* check the 8 flow directions for upslope elements 
+          /* check the 8 flow directions for upslope elements
              without a topidx value */
 
           not_yet = 0;
 
           for(jj=-1; jj < 2; jj++){
             for(ii=-1; ii < 2; ii++){
               if(((i+ii >= 0) && (i+ii < *nrow) && (j+jj >= 0) && (j+jj < *ncol))
-                 && ((ii != 0) || (jj != 0)) 
+                 && ((ii != 0) || (jj != 0))
                  && (dem[i+ii][j+jj] != exclude)) {
                 if((dem[i+ii][j+jj] > dem[i][j]) && (atb[i+ii][j+jj] < ZERO))
                   not_yet = 1;
@@ -129,7 +129,7 @@ void c_topidx(double *inputdem,
           /* if there are no upslope elements without a topidx value,
              start calculations */
 
-          /* find the outflow direction and calculate the sum of weights using 
+          /* find the outflow direction and calculate the sum of weights using
              (tanb*countour length). Contour length = 0.5dx for the cardinal
              direction and 0.354dx for diagonal */
 
@@ -145,7 +145,7 @@ void c_topidx(double *inputdem,
           for(jj=-1; jj < 2; jj++){
             for(ii=-1; ii < 2; ii++){
               if(((i+ii >= 0) && (i+ii < *nrow) && (j+jj >= 0) && (j+jj < *ncol))
-                 && ((ii != 0) || (jj != 0)) 
+                 && ((ii != 0) || (jj != 0))
                  && (dem[i+ii][j+jj] != exclude)) {
                 if((ii == 0) || (jj == 0)) {
                   dnx = dx1;
@@ -172,7 +172,7 @@ void c_topidx(double *inputdem,
 	/* if a sink or a river cell... */
 
         if((nrout == 0) || (river == 1)) {
-	        nsink++;
+	        // nsink++;
           river = 0;
 
           /* assume that there is a channel of length dx running midway through
@@ -184,7 +184,7 @@ void c_topidx(double *inputdem,
           for(jj=-1; jj < 2; jj++){
             for(ii=-1; ii < 2; ii++){
               if(((i+ii >= 0) && (i+ii < *nrow) && (j+jj >= 0) && (j+jj < *ncol))
-                 && ((ii != 0) || (jj != 0)) 
+                 && ((ii != 0) || (jj != 0))
                  && (dem[i+ii][j+jj] != exclude)) {
                 if((ii == 0) || (jj == 0)) dnx = dx1;
                 else dnx = dx2;
@@ -237,7 +237,7 @@ void c_topidx(double *inputdem,
 	for(jj=-1; jj < 2; jj++){
 	  for(ii=-1; ii < 2; ii++){
 	    if(((i+ii >= 0) && (i+ii < *nrow) && (j+jj >= 0) && (j+jj < *ncol))
-	       && ((ii != 0) || (jj != 0)) 
+	       && ((ii != 0) || (jj != 0))
 	       && (atb[i+ii][j+jj] != exclude)) {
 	      if(routdem[nrout] > 0) area[i+ii][j+jj] += c * routdem[nrout];
 	    }
@@ -261,7 +261,3 @@ void c_topidx(double *inputdem,
   }
   return;
 }
-
-
-
-