Terrain attributes can be derived by elevation values on a small neighbourhood of each DEM point.
The major attributes can be calculated from the derivatives of the topographic surface. These derivatives measure the rate of change in elevation in relation with the point location.
In literature are available several methods to compute these derivatives.
In this function the following methods have been implemented:
Evans (1980) - In this method a quadratic surface is fitted on a 3x3 elevation window and the derivatives are calculated using all the 9 cells in the windows.
Zevenbergen and Thorne (1987) - Here a more complex Lagrange polynomial is fitted to the topographic surface.
Moore et al. (1993) and Shary (1995), which are similar to the previous two. For more info refer to Florinsky (1998).
This function works with the package raster.
It takes 3 arguments:
data: the raster object
attr: the terrain attribute to be computed. At the moment the choice is between: "slope", "aspect", "plan.curvature" and "prof.curvature".
method: one of the four method implemented: The choice is between: "evans", "zev.tho", "moore" and "shary".
Ex.
DTM<-raster("c:/raster.asc")
slope<-DEMderiv(DTM,"slope","evans")
The output is a raster object.
This is the code:
DEMderiv Function
This method does not work since the raster packages has been updated. In the latest version, the function focal does not produce a list of the raster values in any given moving window, but a list of weights that sum up to 1. For this reason a call to the function presented here will not provide the results expected.
ReplyDeleteIf someone is willing to try and fix it is encouraged to do so.
Here is a function I used to calculate plan and profile curvature using the raster and sp packages in R. I modified the code provided in this post so the function works as of November 2017. I hope it works for the data of others as well, or at least can help create a function that works for you:
ReplyDeletecurvature <- function(dat, calc = "both"){
# Function to calculate plan and profile curvatures from a DEM.
# Author: Katriina O’Kane, UBC Geography, Vancouver
# dat = DEM in raster form
# calc = what calculation you want to perform, either “plan” for plan curvature, “profile” for profile curvature, or “both” to calculate both plan and profile curvature at the same time.
#Define required vetors
ans.pl <- c()
ans.pr <- c()
ans.x <- c()
ans.y <- c()
#Calculate
for (i in 1:length(dat@data@values)){
a <- adjacent(dat, cells=i, directions=8, pairs=TRUE)
a <- a[,2]
#Define cells
C = res(dat)[1]
Z0 = as.numeric(dat[i])
Z1 = as.numeric(dat[a[1]])
Z2 = as.numeric(dat[a[7]])
Z3 = as.numeric(dat[a[4]])
Z4 = as.numeric(dat[a[2]])
Z5 = as.numeric(dat[a[5]])
Z6 = as.numeric(dat[a[3]])
Z7 = as.numeric(dat[a[8]])
Z8 = as.numeric(dat[a[6]])
#Define values
D = ((Z4 + Z5)/2 - Z0)/C^2
E = ((Z2 + Z7)/2 - Z0)/C^2
f = (Z3 - Z1 + Z6 - Z8)/(4*C^2)
G = (Z5 - Z4) / (2*C)
H = (Z2 - Z7) / (2*C)
#Calculate curvatures
if (calc == "plan") {
ans.pl = c(ans.pl, (2*(D*H^2 + E*G^2 - f*G*H)) / (G^2 + H^2))
ans.x <- c(ans.x, as.numeric(coordinates(dat)[i,1]))
ans.y <- c(ans.y, as.numeric(coordinates(dat)[i,2]))
} else if (calc == "profile") {
ans.pr = c(ans.pr, (-2*(D*G^2 + E*H^2 + f*G*H)) / (G^2 + H^2))
ans.x <- c(ans.x, as.numeric(coordinates(dat)[i,1]))
ans.y <- c(ans.y, as.numeric(coordinates(dat)[i,2]))
} else if (calc == "both") {
ans.pl <- c(ans.pl, (2*(D*H^2 + E*G^2 - f*G*H)) / (G^2 + H^2))
ans.pr <- c(ans.pr, (-2*(D*G^2 + E*H^2 + f*G*H)) / (G^2 + H^2))
ans.x <- c(ans.x, as.numeric(coordinates(dat)[i,1]))
ans.y <- c(ans.y, as.numeric(coordinates(dat)[i,2]))
} else {
stop(return ("Error, did not specify calc properly"))
}
}
#Combine into data frams
if (length(ans.pl) > 0 & length(ans.pr) == 0) {
ans <- data.frame("x" = ans.x, "y" = ans.y, "plan" = ans.pl)
} else if (length(ans.pl) == 0 & length(ans.pr) > 0) {
ans <- data.frame("x" = ans.x, "y" = ans.y, "profile" = ans.pr)
} else if (length(ans.pl) > 0 & length(ans.pr) > 0) {
ans <- data.frame("x" = ans.x, "y" = ans.y, "plan" = ans.pl, "profile" = ans.pr)
}
#Return answer
return(ans)
}
Hi Kaia
ReplyDeleteI'm testing your code but R studio give me immediately an error:
>>: In .doExtract(x, i, drop = drop) : some indices are invalid (NA returned)
What do you think about the message?
Thanks in advance
Simo