[R-sig-Geo] Implementation of algorithms

[Ashton Shortridge]

Dear Jorge,

there are a lot of ways to do this. In straight R, here's a function designed 
to work on a matrix of elevation values that I have used for many years. It is 
very similar to the code you provide:

# This function calculates slope for a matrix using the Horn (1981) 
# method, as reported by Burrough & McDonnell (1998). Slope is not 
# calculated for the outermost rows and columns of the input matrix. 
# Parameters: z is an input matrix of elevations, and spacing is 
# the distance in ground units between cell centers.
# Returns a matrix of slope, in percent. 
#        written by Ashton Shortridge, 3-30-2001, 3/2004, 3/2008, 11/2009
# References:
# Burrough, P. A., McDonnell, R. A. (1998) Principles of GIS. Oxford Press: 
Oxford. 
# Horn, B. K. P. (1981) Hill shading and the reflectance map. Proc. IEEE 69(1): 
14-47.

calc.slope <- function (z, spacing)
{
   # Takes a matrix z of elevations. Returns a matrix of slope, in percent.
   rows <- dim(z)[1] - 2
   cols <- dim(z)[2] - 2
   slope.mat <- matrix(nrow=rows+2, ncol=cols+2)
   for(i in 2:rows) {
      for(j in 2:cols) {
         dzdx <- (z[i+1,j+1] + 2*z[i+1,j] + z[i+1,j-1]) - (z[i-1,j+1] + 
2*z[i-1,j] + z[i-1,j-1])
         dzdy <- (z[i+1,j+1] + 2*z[i,j+1] + z[i-1,j+1]) - (z[i+1,j-1] + 
2*z[i,j-1] + z[i-1,j-1])
         dzdx <- dzdx/(8*spacing)
         dzdy <- dzdy/(8*spacing)
         slope.mat[i,j] <- sqrt(dzdx^2 + dzdy^2) * 100
      }
   }
   return(slope.mat)
}



On 11/28/11, Jorge Cabrera, wrote:
> Hello everyone in the forum
> For introducing myself I would say I have a basic knowledge of R.
> Now I am intended to implement a flood algorithm using R. I have some
> MatLab experience doing these, and as an example to explain more or less
> what I want, I have a m code to calculate the slope from a Digital
> elevation model:
> 
> 
> slope=zeros(row,col);
> for i=2:row-1
>     for j=2:col-1
>        dzdx=(raster(i,j+1)-raster(i,j-1))/(2*res);
>        dzdy=(raster(i+1,j)-raster(i-1,j))/(2*res);
>        slope(i,j)=sqrt(dzdx^2+dzdy^2);
>     end
> end;
> 
> The question is to know how to do the similar procedure on R.
> All suggestions are welcome Thanks
> 
> All best,
> 
> Jorge
> PD:I am using R on windows system 64 bits
> 
> 	[[alternative HTML version deleted]]


-----
Ashton Shortridge
Associate Professor			ashton at msu.edu
Dept of Geography			http://www.msu.edu/~ashton
235 Geography Building		ph (517) 432-3561
Michigan State University		fx (517) 432-1671