tag:blogger.com,1999:blog-1442302563171663500.post5276310906213424281..comments2023-10-14T15:58:07.881+02:00Comments on R tutorial for Spatial Statistics: Terrain Attributes with the raster packageFabio Veronesihttp://www.blogger.com/profile/07827549157455488947noreply@blogger.comBlogger3125tag:blogger.com,1999:blog-1442302563171663500.post-27158572380981642018-07-19T11:54:43.344+02:002018-07-19T11:54:43.344+02:00Hi Kaia
I'm testing your code but R studi...Hi Kaia<br /><br /> I'm testing your code but R studio give me immediately an error:<br /><br />>>: In .doExtract(x, i, drop = drop) : some indices are invalid (NA returned)<br /><br />What do you think about the message?<br />Thanks in advance<br />Simo<br />Anonymoushttps://www.blogger.com/profile/05025329292059394293noreply@blogger.comtag:blogger.com,1999:blog-1442302563171663500.post-85107286382758981922018-01-21T21:16:36.934+01:002018-01-21T21:16:36.934+01:00Here is a function I used to calculate plan and pr...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:<br /><br />curvature <- function(dat, calc = "both"){<br /> # Function to calculate plan and profile curvatures from a DEM.<br /> # Author: Katriina O’Kane, UBC Geography, Vancouver<br /> # dat = DEM in raster form<br /> # 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.<br /> <br /> #Define required vetors<br /> ans.pl <- c()<br /> ans.pr <- c()<br /> ans.x <- c()<br /> ans.y <- c()<br /> <br /> #Calculate<br /> for (i in 1:length(dat@data@values)){ <br /> <br /> a <- adjacent(dat, cells=i, directions=8, pairs=TRUE) <br /> a <- a[,2]<br /> <br /> #Define cells<br /> C = res(dat)[1]<br /> Z0 = as.numeric(dat[i])<br /> Z1 = as.numeric(dat[a[1]])<br /> Z2 = as.numeric(dat[a[7]])<br /> Z3 = as.numeric(dat[a[4]])<br /> Z4 = as.numeric(dat[a[2]])<br /> Z5 = as.numeric(dat[a[5]])<br /> Z6 = as.numeric(dat[a[3]])<br /> Z7 = as.numeric(dat[a[8]])<br /> Z8 = as.numeric(dat[a[6]])<br /> <br /> #Define values <br /> D = ((Z4 + Z5)/2 - Z0)/C^2<br /> E = ((Z2 + Z7)/2 - Z0)/C^2<br /> f = (Z3 - Z1 + Z6 - Z8)/(4*C^2)<br /> G = (Z5 - Z4) / (2*C)<br /> H = (Z2 - Z7) / (2*C)<br /> <br /> #Calculate curvatures<br /> if (calc == "plan") {<br /> ans.pl = c(ans.pl, (2*(D*H^2 + E*G^2 - f*G*H)) / (G^2 + H^2))<br /> ans.x <- c(ans.x, as.numeric(coordinates(dat)[i,1]))<br /> ans.y <- c(ans.y, as.numeric(coordinates(dat)[i,2]))<br /> } else if (calc == "profile") {<br /> ans.pr = c(ans.pr, (-2*(D*G^2 + E*H^2 + f*G*H)) / (G^2 + H^2))<br /> ans.x <- c(ans.x, as.numeric(coordinates(dat)[i,1]))<br /> ans.y <- c(ans.y, as.numeric(coordinates(dat)[i,2]))<br /> } else if (calc == "both") {<br /> ans.pl <- c(ans.pl, (2*(D*H^2 + E*G^2 - f*G*H)) / (G^2 + H^2))<br /> ans.pr <- c(ans.pr, (-2*(D*G^2 + E*H^2 + f*G*H)) / (G^2 + H^2))<br /> ans.x <- c(ans.x, as.numeric(coordinates(dat)[i,1]))<br /> ans.y <- c(ans.y, as.numeric(coordinates(dat)[i,2]))<br /> } else {<br /> stop(return ("Error, did not specify calc properly"))<br /> }<br /> }<br /> <br /> #Combine into data frams<br /> if (length(ans.pl) > 0 & length(ans.pr) == 0) {<br /> ans <- data.frame("x" = ans.x, "y" = ans.y, "plan" = ans.pl)<br /> } else if (length(ans.pl) == 0 & length(ans.pr) > 0) {<br /> ans <- data.frame("x" = ans.x, "y" = ans.y, "profile" = ans.pr)<br /> } else if (length(ans.pl) > 0 & length(ans.pr) > 0) {<br /> ans <- data.frame("x" = ans.x, "y" = ans.y, "plan" = ans.pl, "profile" = ans.pr)<br /> }<br /> <br /> #Return answer<br /> return(ans)<br />}<br />KaiaOhttps://www.blogger.com/profile/10761637550723060731noreply@blogger.comtag:blogger.com,1999:blog-1442302563171663500.post-89406141873979620642013-11-22T10:22:13.935+01:002013-11-22T10:22:13.935+01:00This method does not work since the raster package...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.<br />If someone is willing to try and fix it is encouraged to do so.Fabio Veronesihttps://www.blogger.com/profile/07827549157455488947noreply@blogger.com