Description Usage Arguments Details Value warning Note Author(s) References See Also Examples

Computes a unit vector normal to every grid cell in a digital elevation model.

1 |

`dem` |
Digital Elevation Model, either a matrix or a Raster Layer. |

`dlx` |
Resolution along X axis (longitude). |

`dly` |
Resolution along Y axis (latitude). |

`cArea` |
Logical, if |

The normal vector to the grid cell is the cross product of vectors along the sides of the polygon that form the grid cell. By definition the length of this vector is equal to the area of the polygon.

Returns a 3D matrix with the 2 first dimensions as input dem
and the 3rd dimension of value 3 corresponding to the x, y , z coordinates
of a unit vector perpendicular to every grid cell.
If cArea is `TRUE`

, the result is a 2D matrix with the surface area of every grid cell.

`dlx`

ad `dly`

are assumed to be constant over the extension of the DEM,
therefore the projection should not be `latlong`

.
In this case the resolution is a constant angle, and the equivalent distance
on the surface changes with latitude, giving incorrrect results.

The returned information for every cell is contained by the node at the upperleft corner and the last column and row are undefined. The values given for the last colum and row are a replication of the previous column and row.

Javier G. Corripio

Corripio, J. G.: 2003, Vectorial algebra algorithms for calculating terrain parameters from DEMs
and the position of the sun for solar radiation modelling in mountainous terrain, *International
Journal of Geographical Information Science* 17(1), 1-23.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 | ```
## visualize x, y z components of vector normal to a DEM representing a regular pyramid
m = matrix(0,nrow=100,ncol=100)
for (i in 1:100){ for (j in 1:100){
m[i,j]=50-max(abs(i-50),abs(j-50)) }}
grdm = cgrad(m,1)
xcomponent = grdm[,,1]
ycomponent = grdm[,,2]
zcomponent = grdm[,,3]
image(t(xcomponent[nrow(xcomponent):1,]) ,col=grey(1:10/10))
image(t(ycomponent[nrow(ycomponent):1,]) ,col=grey(1:10/10))
image(t(zcomponent[nrow(zcomponent):1,]) ,col=grey(1:10/10))
## Surface area of every grid cell in a mountain region
## Steep slopes correspond to larger surface area per grid cell
zipfile = tempfile()
download.file("https://meteoexploration.com/R/insol/data/dempyrenees.asc.zip",zipfile)
header = read.table(unz(zipfile,'dempyrenees.asc'),nrows=6)
dem = read.table(unz(zipfile,'dempyrenees.asc'),skip=6)
dem = as.matrix(dem)
unlink(zipfile)
cellsize = header[5,2]
grdarea = cgrad(dem,cellsize,cArea=TRUE)
image(t(grdarea[nrow(grdarea):1,]),col=grey(100:1/100))
## plot the relationship between surface slope and surface area per grid cell:
slpg = slope(cgrad(dem,cellsize),degrees=TRUE)
plot(slpg,grdarea)
## This would be a linear relationship in 2D,
## but in a DEM the slope of the grid cell depends on 4 points in 3D
plot(tan(radians(slpg)),grdarea)
## Not run:
## Use raster for better display and keep the dem projection
require(rgdal)
require(raster)
demfile = tempfile()
download.file("https://meteoexploration.com/R/insol/data/dempyrenees.tif",demfile)
dem = raster(demfile)
plot(dem)
grd = cgrad(dem)
grdarea = cgrad(dem,cArea=TRUE)
rgrdarea = raster(grdarea,crs=projection(dem))
extent(rgrdarea) = extent(dem)
plot(rgrdarea,col = grey(100:1/100))
contour(dem,col='sienna1',lwd=.5,nlevels=30,add=TRUE)
unlink(demfile)
## End(Not run)
``` |

insol documentation built on Jan. 16, 2019, 5:04 p.m.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.