Merge branch 'master' of https://github.com/LandSciTech/roads

R/getGraph.R

@@ -159,3 +159,28 @@ getGraph<- function(sim, neighbourhood,method="old",weightFunction = function(x1
     return(invisible(g))
   }
 }
+
+
+
+#' Grade penalty edge weight function
+#'
+#' Method for calculating the weight of an edge between two nodes
+#' from the value of the input raster at each of those nodes (x1 and x2), designed for DEM input.
+#' This is a simplified version of the grade penalty approach taken by Anderson and Nelson:
+#' doesn't distinguish between adverse and favourable grades.
+#' construction cost values from interior appraisal manual.
+#' ignores (unknown) grade penalties beside roads in order to make do with a single input layer.
+#'
+#' @param x1,x2 Value of the input raster at two nodes. A difference of 1 implies a 100% slope.
+#' @param limit Maximum grade (%) on which roads can be built.
+#' @param penalty Cost increase associated with each additional % increase in road grade.
+#' @noRd
+slopePenaltyFn<-function(x1,x2,limit=10,penalty=504){
+  grade = 100*abs(x1-x2)*(pmin(x1,x2)!=0) #percent slope. if one of the locations is a road, don't apply grade penalty
+  slp = 16178+grade*penalty
+  slp[grade>limit]=NA
+  slp[pmax(x1,x2)==0]=0 # if both 0 then this is an existing road link
+
+  return(slp)
+}
+

tests/testthat/test-AltGraphFns.R

@@ -62,6 +62,7 @@ if(FALSE){
   #octagon methods yield similar but not identical graphs because gdistance::geoCorrection
   #method accounts for lengths of diagonals and other geographic distortions on the grid.
 
+
   #speed/memory benchmarking
   data_path_raw <- "~/Documents/gitprojects/RoadPaper/analysis/data/raw_data/"
   out_path <- "~/Documents/gitprojects/RoadPaper/analysis/figures/"
@@ -82,12 +83,16 @@ if(FALSE){
   neighb  ="octagon"
   bm <- bench::mark(min_iterations = 1, check = FALSE,
                     old_100 = getGraph(landscape_100,neighb,method="old"),
+                    dem_100 = getGraph(landscape_100,neighb,method="old",weightFunction = slopePenaltyFn),
                     gdistance_100 = getGraph(landscape_100,neighb,method="gdistance"),
-                    old_500 = getGraph(landscape_500,"rook",method="old"),
+                    old_500 = getGraph(landscape_500,neighb,method="old"),
+                    dem_500 = getGraph(landscape_500,neighb,method="old",weightFunction=slopePenaltyFn),
                     gdistance_500 = getGraph(landscape_500,neighb,method="gdistance"),
-                    old_1000 = getGraph(landscape_1000,"rook",method="old"),
+                    old_1000 = getGraph(landscape_1000,neighb,method="old"),
+                    dem_1000 = getGraph(landscape_1000,neighb,method="old",weightFunction=slopePenaltyFn),
                     gdistance_1000 = getGraph(landscape_1000,neighb,method="gdistance"),
-                    old_2000 = getGraph(landscape_2000,"rook",method="old"),
+                    old_2000 = getGraph(landscape_2000,neighb,method="old"),
+                    dem_2000 = getGraph(landscape_2000,neighb,method="old",weightFunction=slopePenaltyFn),
                     gdistance_2000 = getGraph(landscape_2000,neighb,method="gdistance")
   )