Working with Raster Data

Check projection system

ELEV.GP@crs

## CRS arguments:
##  +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
## +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Check raster attribute

ELEV.GP

## class       : RasterLayer 
## dimensions  : 310, 275, 85250  (nrow, ncol, ncell)
## resolution  : 5000, 5000  (x, y)
## extent      : -1250285, 124715, 988795, 2538795  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : F:\Spatial_Data_Processing_and_Analysis_R\Data\DATA_07\GP_ELEV.tif 
## names       : GP_ELEV 
## values      : 224.4348, 3828.2  (min, max)

• Nrow, Ncol: the number of rows and columns in the data (imagine a spreadsheet or a matrix).
• Ncells: the total number of pixels or cells that make up the raster.
• Resolution: the size of each pixel (deg decimal).
• Extent: the spatial extent of the raster. This value will be in the same coordinate units as the coordinate reference system of the raster.
• Coord ref: the coordinate reference system string for the raster. This raster is in Geographical Coordinate system with a datum of WGS 84.

Raster extent

extent(ELEV.GP)

## class       : Extent 
## xmin        : -1250285 
## xmax        : 124715 
## ymin        : 988795 
## ymax        : 2538795

Raster statistics

R allow you to calculate some basic statistics of raster data. We can use summary function to get some statistics of raster layer.

summary(ELEV.GP)

##            GP_ELEV
## Min.      224.4348
## 1st Qu.  1207.6884
## Median   1633.6725
## 3rd Qu.  2160.3500
## Max.     3828.2000
## NA's    42766.0000

Histogram of raster

R-base function hist is also useful to know frequency of raster values

HISTO.GP<-hist(ELEV.GP, 
     breaks=6,  
     col= "purple",                                              # color of histogram 
     main="Distribution of Elevation values with 5 breaks",     # tittle of the plot
     xlab= "Elevation (m)",                                     # x-axis name
     ylab = "Number of Pixels")                                 # y-axis name
box()

Plot raster

plot(ELEV.GP, main="Albers Equal Area Conic")

Plot Clip DEM

plot(ELEV.CO, main= "5 km x 5 Km DEM of CO State")
plot(CO.STATE, add=TRUE)

Histogram

Before reclassification, we create a histogram with 4 breaks, of ELEV.CO raster

HISTO.CO<-hist(ELEV.CO, 
     breaks=3,  
     col= "purple",                                             # color of histogram 
     main="Distribution of Elevation values with 4 breaks",     # tittle of the plot
     xlab= "Elevation (m)",                                     # x-axis name
     ylab = "Number of Pixels")                                 # y-axis name
box()

Now, you can check where are breaks and how many pixels in each category, and use these breaks to plot a raster layer

HISTO.CO$breaks

## [1] 1000 2000 3000 4000

HISTO.CO$counts

## [1] 5347 4083 1335

# plot using  the histogram breaks.
plot(ELEV.CO, 
     breaks = c(1000, 2000, 3000, 4000), 
     col = terrain.colors(5),
     main="DEM ")

We will use breaks of histogram to reclassify the DEM raster. First, you have to create classification matrix

• Class 1: 1000 - 2000 m (low elevation)
• Class 2: 2000 - 3000 m (medium elevation)
• Class 3: 3000 - 4000 (Inf) m (high elevation)

# Create a dataframe
reclass_df <- c(1000, 2000, 1,
                2000, 3000, 2,
                3000, Inf, 3)
reclass_df

## [1] 1000 2000    1 2000 3000    2 3000  Inf    3

Next, you have to reshape the data frame to matrix using matrix function

reclass_m <- matrix(reclass_df,
                ncol = 3,
                byrow = TRUE)
reclass_m

##      [,1] [,2] [,3]
## [1,] 1000 2000    1
## [2,] 2000 3000    2
## [3,] 3000  Inf    3

Now you have to apply reclassify() function of raster package using re-class_m object.

ELEV.CO.reclass <- reclassify(ELEV.CO,
                     reclass_m)

You can view the distribution of pixels assigned to each class using the bar-plot().

barplot(ELEV.CO.reclass,
        main = "Number of pixels in each class")
box()

Now we plot the reclassified raster using plot function with a nice looking legend

plot(ELEV.CO.reclass,
     legend = FALSE,
     col = c("green", "blue", "red"), axes = FALSE,
     main = "Reclassified DEM")

legend("topright",
       legend = c("Low", "Medium", "High"),
       fill = c("green", "blue", "red"),
       border = FALSE,
       bty = "n") # turn off legend border

Quantile

Sometimes quantile function is us full to the distribution of raster values

quantile(ELEV.CO)                     #  default 5 quantile

##       0%      25%      50%      75%     100% 
## 1032.333 1490.708 2006.680 2608.840 3828.200

quantile(ELEV.CO, prob=seq(0,1,0.1))  # 10 quantile

##       0%      10%      20%      30%      40%      50%      60%      70% 
## 1032.333 1280.194 1418.156 1567.555 1770.616 2006.680 2267.440 2480.909 
##      80%      90%     100% 
## 2756.125 3089.469 3828.200

value.vector <- values(ELEV.CO)
value.vector <- value.vector[value.vector != 0]
at = round(classIntervals(value.vector, n = 10, style = "quantile")$brks, 0) # 10 quantile
at

##  [1] 1032 1280 1418 1568 1771 2007 2267 2481 2756 3089 3828

rgb.palette <- colorRampPalette(c("lightblue","blue","green","yellow","red"),
space = "rgb")
polys<- list("sp.lines", as(CO.STATE, "SpatialLines"), col="grey", lty=1)  # CO state boundary
spplot(ELEV.CO, main = "Quantile of DEM (m)",
   at=at,sp.layout=list(polys),
   par.settings=list(axis.line=list(col="darkgrey",lwd=1)),
   colorkey=list(space="right",height=1, width=1.3,at=1:11,labels=list(cex=1.0,at=1:11,
   labels=c("1032","1280","1418","1568","1771","2007","2267","2481","2756","3089","3828"))),
   col.regions=rgb.palette(100))

Slope

The slope is an incline, or steepness, of a surface. The slope can be measured in degrees from horizontal (0-90), or percent slope (which is the rise divided by the run, multiplied by 100). A slope of 30 and 45 degrees equals 58 and 100 percent slope, respectively. As the slope angle approaches vertical (90 degrees), the percent slope approaches infinity. The slope for a cell in a raster is the steepest slope of a plane defined by the cell and its eight surrounding neighbors. (Source: ESRI online GIS Dictionary).

Slope(Source: http://resources.arcgis.com/EN/HELP/MAIN/10.2/index.html#/How_Focal_Statistics_works/009z000000r7000000/)

We will use terrain() function from raster package to create raster slope.

SLOPE <- terrain(ELEV.CO,           # raster DEM
    opt="slope",                        # Slope function 
    unit='degrees')                 # unit degree

plot(SLOPE,
       col = rainbow(16), main= "Slope (degree)")

Aspect

The compass direction that a slope faces is most commonly measured from degrees north. The Aspect command provides a 32-bit float raster that ranges from 0° and 360°. This represents the azimuth (the angular distance from the north to south point of the horizon that can be transversed by a vertical circle that intersects the horizon, which represents the direction of an object from the observer) that slopes are facing. In other words, 0° means that the slope is facing North, 90° indicates East, 180° indicates a Southward facing, and 270° is to the west. This however requires that the top of your input raster is north oriented, as most are. Values found in between these ranges assume a mixture of cardinal directions, ex. 250 indicates a hillside with a Southeast facing aspect. The aspect value -9999 is generally used as the ‘nodata’ value to indicate use of an undefined aspect in areas lacking variation in topography, with slope=0. You can also make use of the legend to indicate the appropriate aspect for each hillside in your study area. The legend permits you to identify North (~ 0 or 360 degrees), South (~ 180 degrees), East (~270 degrees) or West (~90 degrees) facing slopes).

Aspect(Source: http://resources.arcgis.com/EN/HELP/MAIN/10.2/index.html#/How_Focal_Statistics_works/009z000000r7000000/)

ASPECT <- terrain(ELEV.CO,        #  raster DEM
    opt="aspect",                       #  aspect function 
    unit='degrees')                   #  unit degree

plot(ASPECT,
       col = rainbow(16), main= "Aspect (degree)")

Hillshade

A hillshade is a grayscale 3D representation of the surface, with the sun’s relative position taken into account for shading the image. This function uses the altitude and azimuth properties to specify the sun’s position.

Hill Shade(Source: http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/hillshade-function.htm)

We generate a hill shade layer from slope and aspect layers (both in degree). Slope and aspect can be computed with function terrain. A hill shade layer is often used as a backdrop on top of which another, semi-transparent, layer is drawn.

hill <- hillShade(SLOPE, ASPECT, 45, 250)
plot(hill, col=grey(0:100/100), legend=FALSE, main='Hillshade Map of CO')
plot(ELEV.CO, col=rainbow(50, alpha=0.35), add=TRUE)

Terrain Ruggedness Index

TRI (Terrain Ruggedness Index) is the mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells (the amount of elevation difference between adjacent cells of a DEM). This provides a relative measure of elevation changes between a specified grid cell and neighbors. The TRI is useful for analyzing what environments might be suited to particular crops or species which may be sensitive to particular sloped environments, or for assessing the potential flow of soil during erosion events, and so on.

TRI<- terrain(ELEV.CO,          #  raster DEM
    opt="TRI",                  #  TRI function 
    unit='degrees')                 #  unit degree

plot(TRI,
       main= "Terrain Ruggedness Index (TRI)",
     col=rainbow(50))

Topographic Position Index

TPI (Topographic Position Index) is the difference between the value of a cell and the mean value of its 8 surrounding cells.Positive TPI values represent locations that are higher than the average of their surroundings, as defined by the neighborhood (ridges). Negative TPI values represent locations that are lower than their surroundings (valleys). TPI values near zero are either flat areas (where the slope is near zero) or areas of constant slope. Topographic position is an inherent scale-dependent phenomenon.

TPI<- terrain(ELEV.CO,          #  raster DEM
    opt="TPI",                  #  TPI function
    unit='degrees')                 #  unit degree

plot(TPI,
       main= "Topographic Position Index (TPI)",
     col=rainbow(50))

Roughness

Roughness is the difference between the maximum and the minimum value of a cell and its 8 surrounding cells.

ROUGH<- terrain(ELEV.CO,          # raster DEM
    opt="roughness",            # roughness function
    unit='degrees')                 # unit degree

plot(ROUGH,
       main= "Roughness",
     col=rainbow(50))

Curvature

Curvature is the change in the degree of slope for a given distance on that slope. It is a very helpful measure to understand the surface water flow through a landscape, and can be used in part to help design an irrigation scheme, for example. It is widely used by hydrological scientists. When interpreting data from a curve, any positive values indicate a convex (upward bulging) slope, while negative values correspond to a convex (downward facing) curve. The first is like a dinner bowl sanding correctly on a table, the latter is like a bowl that has been flipped over and is facing down. Any pixels observed that show positive curvature indicate the potential for liquid flow dispersal away from a central area, while and negative values indicate accumulation. In combination, convex and concave surfaces approximate actual topography.

A Plan curve is the measurement of the rate at which a horizontal curve changes. Positive and negative values respectively indicate divergence and convergence of a particular slope. Highlight a divergent slope, negative values a convergent slope.

A Profile curve is conversely a vertical measurement of the change of a slope. Positive values indicate convexity, while negative values represent concave profiles.

Hill Shade(Source: http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/hillshade-function.htm)

In there is no package available yet for Curvature analysis in R. I have a function here. We will call this function using source(“curveture_function.R”), But when I run this, it shows error.

#source("curveture_function.R")
# profile<-curv(DEM, calc = "profile")
# plan<-curv(DEM, calc = "plan")

However, in spatialEco package allow to calculate McNab’s and Bolstad’s variants of the surface curvature (concavity/convexity) index (McNab 1993; Bolstad & Lillesand 1992; McNab 1989). The index is based on features that confine the view from the center of a 3x3 window. In the Bolstad equation, edge correction is addressed be dividing by the radius distance to the outermost cell (36.2m). We will use curveture function of spatialEco package.

# McNab curvature
m.crv <- curvature(ELEV.CO, s=5, type="mcnab")
# Bolstad curvature
b.crv <- curvature(ELEV.CO, s=5, type="bolstad")

p1<-spplot(m.crv, main="McNab curvature")
p2<-spplot(b.crv, main="Bolstad curvature")
grid.arrange(p1, p2, ncol=2)

Flow Direction

The direction of the greatest drop in elevation (or the smallest rise if all neighbors are higher).

FLD<- terrain(ELEV.CO,          #  raster DEM
    opt="flowdir ",             # flow direction
    unit='degrees')                 # unit degree

plot(FLD,
       main= "Flow Direction",
     col=rainbow(50))

rm(list = ls())