FIX: calc_subbas(): fixed various issues remaining after latest changes

- delineating only one subbasin is possible (there was an error occurring)
- output message regarding the number of subbasins left corrected
- no spurious warning message anymore if number of final subbasin categories and given drainge points differ
- optional column 'subbas_id' in 'drain_points' now corresponds to the final subbasin categories

R/subbasin.R

@@ -25,9 +25,11 @@
 #'      close to boundaries may occur.
 #' @param drain_points \code{SpatialPoints} object containing drainage locations in
 #'      units of and compliant with the projection of your respective GRASS location.
-#'      Can be e.g. be imported from GRASS with
-#'      drain_points = readVECT(vname = "subbas_outlets", layer=1)
-#'      At least the watershed drainage point has to be given. If it contains the attribute 'subbas_id', this will be used for numbering the subbasins.
+#'      Can, e.g., be imported from GRASS with
+#'      drain_points = readVECT(vname = "subbas_outlets", layer=1).
+#'      At least the watershed drainage point has to be given. If it contains column
+#'      'subbas_id' in the attribute table, this will be used for numbering the subbasins.
+#'      IDs of additionally delineated subbasins (if \code{thresh_sub != NULL}) will be appended.
 #' @param river River vector map in GRASS location if available. If set to \code{NULL}
 #'      (default value) river network will be calculated by GRASS function \emph{r.watershed}.
 #' @param flowaccum (optional) Existing raster map of flow accumulation in GRASS location. Must correspond to map generated by \code{r.watershed}.
@@ -274,11 +276,13 @@ calc_subbas <- function(
 
     ### ensure that given drainage points are NOT precisely at cell centres, because this may cause pathologic
     ### cases when snapping to streams (ending up at cell corners instead of cell interior)
-    # add data slot and column subbas_id, if not given
+    # add data slot and columns subbas_id and cat, if not given
     if(!any(slotNames(drain_points) == "data"))
       drain_points <- SpatialPointsDataFrame(drain_points, data=data.frame(subbas_id=1:length(drain_points)))
     if(!any(colnames(drain_points@data) == "subbas_id"))
       drain_points@data <- cbind(drain_points@data, subbas_id=1:length(drain_points))
+    if(!any(colnames(drain_points@data) == "cat"))
+      drain_points@data <- cbind(drain_points@data, cat=1:length(drain_points))
     # force conversion to numeric 
     drain_points@data$subbas_id=as.numeric(as.character(drain_points@data$subbas_id)) 
     if (any(!is.finite(drain_points$subbas_id)))
@@ -406,6 +410,8 @@ calc_subbas <- function(
         outs <- t(outs)
         colnames(outs) <- c("x", "y", "cat")
         outs <- as.data.frame(outs)
+        # subbas_id, make sure they are distinct from drain_points_snap
+        outs <- cbind(outs, subbas_id = seq(1,1e6,by=1)[-drain_points_snap@data$subbas_id][1:nrow(outs)])
         coordinates(outs) <- c("x", "y")
 
         # as SPDF
@@ -414,8 +420,8 @@ calc_subbas <- function(
         # write to GRASS location
         writeVECT(drain_points_calc, paste0(points_processed, "_calc"), ignore.stderr = T, v.in.ogr_flags = "o")
 
-        # df should only contain a cat column
-        drain_points_snap@data <- data.frame(cat=1:nrow(drain_points_snap@data))
+        # drain_points_snap df requires columns 'cat' and 'subbas_id' as drain_points_calc
+        drain_points_snap@data <- drain_points_snap@data[,c("cat", "subbas_id")]
         suppressWarnings(proj4string(drain_points_snap) <- CRS(getLocationProj()))
 
         # merge with existing drain points object (snapped points first as there the outlet is identified)
@@ -432,32 +438,45 @@ calc_subbas <- function(
       warning("Duplicated subbas_id in drainage points. Ignoring IDs, using row numbers instead.")
       drain_points@data$subbas_id=NULL 
     }
+
+    # combined drain points as raster (easier to identify double drain points sharing one raster cell)
+    suppressWarnings(writeVECT(drain_points_snap, paste0(points_processed, "_all_t"), v.in.ogr_flags = "o"))
+    # WINDOWS PROBLEM: delete temporary file otherwise an error occurs when calling writeVECT or readVECT again with the same (or a similar) file name 
+    if(.Platform$OS.type == "windows") {
+      dir_del <- dirname(execGRASS("g.tempfile", pid=1, intern=TRUE, ignore.stderr=T))
+      files_del <- grep(substr(paste0(points_processed, "_all_t"), 1, 8), dir(dir_del), value = T)
+      file.remove(paste(dir_del, files_del, sep="/"))
+    }
+    x <- execGRASS("v.to.rast", input=paste0(points_processed, "_all_t"), output=paste0(points_processed, "_all_t"), use="attr", attribute_column="subbas_id", flags="overwrite", intern=T)
+    x <- execGRASS("r.mapcalc", expression=paste0(points_processed, "_all_t=round(", points_processed, "_all_t)"), flags = "overwrite", intern=T)
+
+    # get coordinates of drain point cells
+    drainp_coords <- execGRASS("r.stats", input = paste0(points_processed, "_all_t"), flags=c("n", "g"), intern=T)
+    drainp_coords <- matrix(as.numeric(unlist(strsplit(drainp_coords, " "))), ncol = 3, byrow = T)
+
     # loop over drainage points of subbasins TODO: This step is slow!
-    for (p in 1:length(drain_points_snap)) {
+    for (p in 1:nrow(drainp_coords)) {
 
       # outlet coordinates
-      outlet_coords <- coordinates(drain_points_snap)[p,]
+      outlet_coords <- drainp_coords[p,c(1,2)]
+
+      # drain points subbasin id (optionally defined in 'drain_points' input)
+      id <- drainp_coords[p,3]
 
       # basin
-      cmd_out <- execGRASS("r.water.outlet", input="drain_t", output=paste0("basin_", p, "_t"), coordinates=outlet_coords, intern = T)
+      cmd_out <- execGRASS("r.water.outlet", input="drain_t", output=paste0("basin_", id, "_t"), coordinates=outlet_coords, intern = T)
 
-      cmd_out = execGRASS("r.stats", input=paste0("basin_", p, "_t"), flag=c("c","n","quiet"), intern = TRUE)
+      cmd_out = execGRASS("r.stats", input=paste0("basin_", id, "_t"), flag=c("c","n","quiet"), intern = TRUE)
       ncells = as.numeric(strsplit(cmd_out, split = " ")[[1]][2])
       if (!is.finite(ncells) | ncells < 100)
-        warning(paste0("Number of cells in calculated catchment ",p," is very low (",ncells,"). Try using a filled DEM and check outlet points."))
-
+        warning(paste0("Number of cells in calculated catchment ",id," is very low (",ncells,"). Try using a filled DEM and check outlet points."))
 
-      # reclass (subbasins gets number of i for crossing later)
-      if (!is.na(drain_points@data$subbas_id[p])) {
-        id = drain_points@data$subbas_id[p]
-      } else { #use specified ID, if available
-        id = p
-      }
-      cmd_out <- execGRASS("r.mapcalc", expression=paste0("basin_recl_", p, "_t = if(basin_", p, "_t,", id, ")"), intern=T)
+      # reclass (for crossing later on)
+      cmd_out <- execGRASS("r.mapcalc", expression=paste0("basin_recl_", id, "_t = if(basin_", id, "_t,", id, ")"), intern=T)
 
     }
 
-    no_catch <- p
+    no_catch <- nrow(drainp_coords)
 
     if(!silent) message(paste("% -> Identified", no_catch, "subbasins."))
     if(!silent) message("% OK")
@@ -469,22 +488,11 @@ calc_subbas <- function(
     if(!silent) message("% Merge calculated catchments...")
 
     # put sub-catchments together
-    #subcatch_rasts <- execGRASS("g.mlist", type="rast", pattern=paste0("basin_recl_[0-9]*_t"), intern=T)
-    subcatch_rasts <- paste0("basin_recl_",1:no_catch, "_t")
+    subcatch_rasts <- paste0("basin_recl_",drainp_coords[,3], "_t")
 
     # if more than one sub-catchment
     if(no_catch > 1) {
 
-      # drain points needed as raster
-      suppressWarnings(writeVECT(drain_points_snap, paste0(points_processed, "_all_t"), v.in.ogr_flags = "o"))
-      # WINDOWS PROBLEM: delete temporary file otherwise an error occurs when calling writeVECT or readVECT again with the same (or a similar) file name 
-      if(.Platform$OS.type == "windows") {
-        dir_del <- dirname(execGRASS("g.tempfile", pid=1, intern=TRUE, ignore.stderr=T))
-        files_del <- grep(substr(paste0(points_processed, "_all_t"), 1, 8), dir(dir_del), value = T)
-        file.remove(paste(dir_del, files_del, sep="/"))
-      }
-      x <- execGRASS("v.to.rast", input=paste0(points_processed, "_all_t"), output=paste0(points_processed, "_all_t"), use="cat", flags="overwrite", intern=T)
-
       # iterate until configuration without 'spurious' sub-catchments is found (if rm_spurious > 0) TODO: This step is slow in case many iterations are needed!
       while (TRUE) {
 
@@ -506,23 +514,6 @@ calc_subbas <- function(
           }
         }
 
-        # x <- execGRASS("r.patch", input=paste(subcatch_rasts, collapse=","), #r.patch does not work under windows
-        #                output=paste0("basin_cross_", j, "_t"), flags = c("overwrite"), intern=T, ignore.stderr=T)
-        # 
-
-        # simply adding up does not do the job as intermediate basins will overwrite headwater basins
-        # x <- execGRASS("g.remove", rast="basin_all_t", intern=T, ignore.stderr=T)
-        # x <- execGRASS("g.copy", rast=paste0(subcatch_rasts[1],",", "basin_all_t"), intern=T, ignore.stderr=T)
-        # 
-        # for (j in 2:length(subcatch_rasts))
-        # {
-        #     execGRASS("r.mapcalculator", amap="basin_all_t", bmap=subcatch_rasts[j], outfile=basin_out,
-        #               formula="if(isnull(A),0,A) + if(isnull(B),0,B)", flags = "overwrite")
-        #     x <- execGRASS("g.copy", rast=paste0(basin_out,",", "basin_all_t"), intern=T, ignore.stderr=T)
-        # } 
-
-
-
         # merge cross products
         cross_rasts <- execGRASS("g.list", type="raster", pattern=paste0("basin_cross_[0-9]*_t"), intern=T)
 
@@ -553,9 +544,16 @@ calc_subbas <- function(
             cmd_cols <- grep("zone|^mean$", cmd_out[[1]])
             basins_points <- do.call(rbind, cmd_out)[-1,cmd_cols, drop=F]
             sub_rm_f <- as.numeric(basins_points[which(as.numeric(basins_points[,1]) %in% sub_rm),2])
-            # remove this temporary map from processing and try again (back to start of while loop)
+            # remove this temporary map from processing and drain points raster map and try again (back to start of while loop)
             subcatch_rasts <- grep(paste0("basin_recl_", sub_rm_f, "_t", collapse="|"), subcatch_rasts, invert = T, value = T)
             x <- execGRASS("g.remove", type="raster", name="basin_all_t", flags = "f", intern=T)
+            tmp_file <- tempfile()
+            stats_t <- as.integer(execGRASS("r.stats", input = paste0(points_processed, "_all_t"), flags=c("n"), intern=T))
+            stats_t <- stats_t[!(stats_t %in% sub_rm_f)]
+            write(paste(paste(sub_rm_f, "NULL", sep = " = ", collapse = "\n"), paste(stats_t, stats_t, sep=" = ", collapse = "\n"), sep="\n"), file=tmp_file)
+            x <- execGRASS("r.reclass", input=paste0(points_processed, "_all_t"), output=paste0(points_processed, "_all2_t"),
+                           rules = tmp_file, flags = "overwrite")
+            x <- execGRASS("r.mapcalc", expression=paste0(points_processed, "_all_t=", points_processed, "_all2_t"), flags="overwrite", intern=T) # convert to regular map
             # update no_catch
             no_catch <- no_catch - length(sub_rm_f)
           } else {
@@ -568,42 +566,27 @@ calc_subbas <- function(
       } # end while-loop
 
       # constrain to catchment of outlet point
-      cmd_out <- execGRASS("r.mapcalc", expression=paste0(basin_out, " = basin_outlet_t * basin_all_t"), intern=T)
+      cmd_out <- execGRASS("r.mapcalc", expression=paste0(basin_out, "_t = basin_outlet_t * basin_all_t"), intern=T)
 
 
     } else { # only one sub-catchment
 
+      cmd_out <- execGRASS("g.copy", raster=paste0("basin_outlet_t,", basin_out, "_t"), intern=T)
+
       if(no_catch == 0)
         stop("Number of identified sub-catchments is zero. Check input data!")
     }
 
-
-
-    # TODO: Only works for special cases where all subbasin-ids were pre-specified. Needs to be generalised!
-    # # reclass automatically-created subbasin-ids to those that were used originally
-    # cmd_out = execGRASS("v.db.addcolumn", map=paste0(points_processed, "_snap"), columns="temp_id integer", ignore.stderr = T, intern = TRUE) #add column for taking up new id
-    # stat = attr(cmd_out, "status")
-    # if (!is.null(stat) && stat== 1)
-    #   stop(paste("Could not add column to", drain_points_snap, ".", sep=" "))
-    # cmd_out = execGRASS("v.what.rast", raster="basin_all_t", map=paste0(points_processed, "_snap"), column="temp_id" ,intern=T, ignore.stderr = T)
-    # stat = attr(cmd_out, "status")
-    # if (!is.null(stat) && stat== 1)
-    #   stop("Could not update column to drain_points_snap.")
-    # drain_points = readVECT(vname = paste0(points_processed, "_snap"), layer=1) #re-import vector layer containing new IDs
-    # rules <- paste0(drain_points$temp_id,"=",drain_points$subbas_id, collapse="\n")
-    # tempfile=tempfile()
-    # write(rules, file=tempfile)
-    # cmd_out = execGRASS("r.reclass", input="basin_all_t", output="basin_all2_t", rules=tempfile, flags="overwrite", intern=T)
-    # unlink(tempfile)
-    # stat = attr(cmd_out, "status")
-    # if (!is.null(stat) && stat== 1)
-    #   stop("Could not reclassify subbasin-map.")
-    # 
-    # execGRASS("g.rename", raster=paste("basin_all2_t", basin_out, sep=","))
-
-
-    # set values of zero to NULL
-    execGRASS("r.null", map=basin_out, setnull="0")
+    # assign correct ids (from 'subbas_id') to basin_out
+    cmd_out <- execGRASS("r.to.vect", input = paste0(points_processed, "_all_t"), output = paste0(points_processed, "_all_t"),
+                         type = "point", column = "subbas_id", flags = "overwrite", intern=T)
+    cmd_out <- execGRASS("v.what.rast", raster=paste0(basin_out, "_t"), map=paste0(points_processed, "_all_t"), column="temp_id" ,intern=T, ignore.stderr = T)
+    drain_points_snap <- readVECT(paste0(points_processed, "_all_t"))
+    dat_rules <- paste(drain_points_snap@data$temp_id, "=", drain_points_snap@data$subbas_id, collapse = "\n")
+    tmp_file <- tempfile()
+    write(dat_rules, file=tmp_file)
+    cmd_out <- execGRASS("r.reclass", input = paste0(basin_out, "_t"), output = paste0(basin_out, "2_t"), rules = tmp_file)
+    x <- execGRASS("r.mapcalc", expression = paste0(basin_out, "=", basin_out, "2_t"), intern = T)
 
     no_cross <- length(execGRASS("r.stats", input=basin_out, flags=c("n"), intern=T, ignore.stderr = T))
     if(no_catch != no_cross) warning(paste0("\nNumber of categories in ", basin_out, " not equal to number of drainage points!\nThis might be because there are drainage points outside the catchment of the defined outlet or due to small inconsistencies between calculated and manually defined (and snapped) drainage points. However, you should check the output with the GRASS GUI and consider the help pages of this function.