egTile: An example tile for debugging and didactic purposes
In henkelstone/NPEL.Classification: Classifaction and data-handling routines for NPEL Caribou Project

Description Usage Arguments Format Details Note Source See Also Examples

Given that it can be a challenge to build a raster dataset there is a mini dataset included with the package. Due to size limitations it does not contain very many layers and is limited to only 100 x 100 pixels. However, it is sufficient to illustrate the data format expected by the package, and it allows for new users to ‘play’ with the package. It also facilitates testing of code because it is so small and code executes relatively quickly.

1	egTile$variable

variable

is the layer to extract from the raster

A 100 x 100 grid containing a few useful Landsat and DEM derived layers

Given the limitations imposed by the need to have a portable library (read: small enough to easily distribute on over the internet), this raster is both small and has only a few layers. The layers included are those that fell out of our dimensional analysis as being the most significant of the Landsat and DEM data. The raw Landsat data is included so that users could generate other derived layers if they wished; code is provided in the example section. Following is a list of included raster variables:

base_1:7 the base Landsat layers
grnns, wetness, brtns the so-called ‘Tasseled-Cap’ variables; they reflect greenness, wetness, and brightness respectively
dem the basic elevations from the digital elevation model
slp slope, as derived from the DEM
asp aspect, as derived from the DEM; note that different algorithms use different values for NA, a reliable filter is that aspect is NA when slope == 0
hsd hillshade, a measure of how much light the surface receives – derived from the DEM

Note that to date, the aspect computed in the terrain package does not match that generated by ArcGIS. At the time of writing (15.Feb.2015) this difference has still not been resolved. Also aspect is tricky because ArcGIS and R have different default methods of dealing with 0 slope values, i.e. aspect is NA. As ever, use any variable with care, and the understanding that it may be flawed.

Also note that due to unresolved path issues, this dataset needs to be created when it is used. See the opening lines of the example code for how to do this.

Landsat data: available from the U.S. Geological Survey
DEM data: available from the U.S. Geological Survey National Elevation Dataset

raster package

library (raster)
library (NPEL.Classification) # Needs to be *before* raster on the search path!
gisPath = './'

egTile <- readTile(file.path(system.file("extdata", "egTile", package = "NPEL.Classification"),''),
                   layers=c('base','grnns','wetns','brtns','dem','slp','asp','hsd'))
summary (egTile)

# Create meaningful aliases for the Landsat layer names
blue <- egTile$'base.2'
green <- egTile$'base.3'
red <- egTile$'base.4'
NIR <- egTile$'base.5'
SWIR1 <- egTile$'base.6'
SWIR2 <- egTile$'base.7'

# Generate Landsat derived data
ndvi <- (NIR - red)/(NIR + red)
ndwi <- (green - NIR)/(green + NIR)
evi <- 2.5 * (NIR-red)/(1+NIR + 6*red + 7.5*blue)
mavi <- 0.5 + NIR - 0.5*sqrt( (1+2*NIR)^2 - 8*(NIR-red) )
nred <- NIR/red
ngrn <- NIR/green
grnns <- -0.1603*blue - 0.2819*green - 0.4934*red + 0.7940*NIR - 0.0002*SWIR1 - 0.1446*SWIR2
brtns <-  0.2043*blue + 0.4158*green + 0.5524*red + 0.5741*NIR + 0.3124*SWIR1 + 0.2303*SWIR2
wetns <-  0.0315*blue + 0.2021*green + 0.3102*red + 0.1594*NIR - 0.6806*SWIR1 - 0.6109*SWIR2

LandSAT <- stack (ndvi, ndwi, evi, mavi, nred, ngrn, grnns, brtns, wetns)
names (LandSAT) <- c('ndvi','ndwi','evi','mavi','nred','ngrn','grnns','brtns','wetns')
rm (blue, green, red, NIR, SWIR1, SWIR2, ndvi, ndwi, evi, mavi, nred, ngrn, grnns, brtns, wetns)
## Not run: writeRaster (LandSAT,filename=paste0(gisPath,'LandSAT'), overwrite=TRUE)

# Generate DEM derived data using built in functions
slope <-  terrain(egTile$dem, opt='slope', unit='degrees')
aspect <- terrain(egTile$dem, opt='aspect', unit='degrees')
aspect[slope < 1e-10] <- NA
aspect[abs(aspect-360) < 1e-10] <- 0
TPI <-    terrain(egTile$dem, opt='TPI')
TRI <-    terrain(egTile$dem, opt='TRI')
rough <-  terrain(egTile$dem, opt='roughness')
flow <-   terrain(egTile$dem, opt='flowdir')
hsd <-    hillShade(deg2rad(slope), deg2rad(aspect), angle=53.3, direction=199.8, normalize=TRUE)

# Generate DEM derived data manually to match ArcGIS output
 RI <- focal(egTile$dem, w=matrix(1, nrow=3, ncol=3),
             fun=function(x){sqrt( sum((x-x[5])^2,na.rm=TRUE)/8 )}, pad=TRUE, padValue=NA)
# SRR <- focal(egTile$dem, w=matrix(1, nrow=3, ncol=3),    # Cute but slow so do it the long way
#              fun=function(x){ m <- min(x); (mean(x)-m)/(max(x)-m) }, pad=TRUE, padValue=NA)
eMin <- focal(egTile$dem, w=matrix(1, nrow=3, ncol=3), fun=min, pad=TRUE, padValue=NA)
eMax <- focal(egTile$dem, w=matrix(1, nrow=3, ncol=3), fun=max, pad=TRUE, padValue=NA)
eAvg <- focal(egTile$dem, w=matrix(1/9, nrow=3, ncol=3), pad=TRUE, padValue=NA)
SRR <- (eAvg-eMin) / (eMax-eMin)
SRR[abs(eMax-eMin) < 0.0001] <- NA
rm (eMin, eMax, eAvg)

DEM <- stack (slope, aspect, TPI, TRI, rough, SRR, RI, flow, hsd)
names (DEM) <- c('slp','asp','TPI','TRI','rough','srr','ri','flow','hsd')
rm (slope, aspect, TPI, TRI, rough, SRR, RI, flow, hsd)
## Not run: writeRaster (DEM,filename=paste0(gisPath,'DEM'), overwrite=TRUE)
detach('package:raster',unload=TRUE)