writeTile: Compute and output a landscape

In henkelstone/NPEL.Classification: Classifaction and data-handling routines for NPEL Caribou Project

Description

This function renders a landscape, which means taking a model and a set of raster layers (stack or brick), then dropping all the pixels through the model to get outputs (see details), and writing the result to a file. Since many raster.* objects are too big to (reasonably) hold in memory, the landscape is broken into blocks and processed in series similar to calc.

Usage

writeTile(model, inRdata, outFilename, layers = c("class"), threshold = 0,
  labels.all = NULL, ...)

Arguments

model	the model to use for predicting output values, often created using generateModels.
inRdata	a raster stack (or brick) of the input data; can be loaded using readTile
outFilename	the output filename
layers	(optional) one or more of: ‘class’ the chosen class—also use for continuous data; ‘prob’ the probability of this class; ‘threshold’ same as class where probability is greater than or equal to threshold (below), otherwise 0; ‘all’ the probabilities for all classes—one layer per class... see Details.
threshold	(optional) this needs to be provided if ‘threshold’ output is selected: 0 < Th < 1.
labels.all	(optional) if 'all' is selected, specify labels for these layers; the default label is "prob.<class name>"
...	variable(s) to pass to the model prediction function.

Details

Given the variety of inputs and outputs that could be used, this function (attempts) to streamline the process for the user by providing a consistent interface, dealing with at least some of the differences internally. Various conditions that influence what type and class of data are returned are: the model being used, whether the data is categorical or continuous, and of course, user preference.

Model packages/types:

• Random Forest models are probabilistic—that is, they output the class that has the largest number of occurrences in each individual tree. Hence they can output a wide variety of different statistics.

• Nearest Neighbour algorithms, by their very nature do not generate class probabilities; they output information about the nearest existing data-point (in phase space). Hence, the option to output probabilities is disabled for the two primitive nearest-neighbour packages (fnn, class).

• The GBM package reports probabilities, so if desired, this can be used as a measure of certainty.

Input data type: for continuous data it is meaningful to talk about the probability of a class occurring, especially in contrast to the probability of another class occurring Of course with continuous data, this doesn't make sense; a model generated on continuous data (that is a regression model) can only report it's best estimate of the continuous variable at each point.

User preference: given those particulars, the user can choose what output they might want:

• For continuous data there isn't much choice... use the 'class' option.

• For categorical data it depends on whether the package/algorithm allows for probabilities. Assuming that these are generated the following output options are offered:

  • ‘class’ will output the ‘winning’ class, that is, the class with the highest probability of occurrence.

  • ‘prob’ the probability with which the ‘winning’ class won.

  • ‘threshold’ if we don't trust the classes with ‘winning’ probabilities below a certain threshold, then this can be selected and classes that don't meet this requirement will be set to 0 in the output map.

  • ‘all’ output all the probability for all the classes as a multi-layer map. Note that this can be an exceptionally large datafile— as in the order of many gigabytes. Plan accordingly.

Value

a new raster.brick object pointing to ‘outFilename’ containing the output landscape.

Warning

Runtimes may be long! On a 2015 iMac (Intel I5 quad-core 3.3 GHz, 8GB RAM, and SSD) it takes 12-15 hours to do a full rendering of approximately 15k x 15k pixels (=225 Mpix). Also, output file sizes can be large(ish) — ~2+ GB.

See Also

writeTiles for an automated way to render all models in a model block

Examples

egTile <- readTile(file.path(system.file("extdata", "egTile", package = "NPEL.Classification"),''),
                   layers=c('base','grnns','wetns','brtns','dem','slp','asp','hsd'))
models <- generateModels (siteData, suppModels, x=c('brtns','grnns','wetns','dem','slp','asp','hsd'), y='ecoType')
fName <- paste0(dirname(tempfile()),'/Tmp_rfsrc.tif')
writeTile (models$rfsrc, egTile, fName, layers='class')
writeTiles (list(models$rfsrc), egTile, base='/Tmp_', path=dirname(fName), layers=c('class','prob','threshold','all'), threshold=0.5)
unlink (fName)

## Not run: 
# Example with the option of cropping to a specific area
data ('Output/siteData.dat')
egTile <- readTile(file.path(system.file("extdata", "egTile", package = "NPEL.Classification"),''),
                   layers=c('base','grnns','wetns','brtns','dem','slp','asp','hsd'))

# Choose a model with specified groupings to load
type <- c('Full','Reduced')[2]                          # Select either full or reduced grouping
load (paste0('Output/',type,' Models/identity.dat'))
load (paste0('Output/',type,' Models/domSpecies.dat'))
load (paste0('Output/',type,' Models/domGroup.dat'))
load (paste0('Output/',type,' Models/MaxGranularity.dat'))
rm (type)

inRdata <- egTile
# Choose a block and crop if desired
IB_1 <- maptools::readShapePoly('Input/Intensive Blocks/IB_1')
IB_2 <- maptools::readShapePoly('Input/Intensive Blocks/IB_2')
IB_3 <- maptools::readShapePoly('Input/Intensive Blocks/IB_3')
IB_4 <- maptools::readShapePoly('Input/Intensive Blocks/IB_4')
IB_5 <- maptools::readShapePoly('Input/Intensive Blocks/IB_5')
IB_6 <- maptools::readShapePoly('Input/Intensive Blocks/IB_6')
inRdata <- crop (inRdata,extent(IB_3))

# Output a tile--be sure the filename matches the parameters... there is no other check!!
writeTile(modelRun$gbm,
          egTile,
          paste0(gisPath,'Renderings/identity_full_gbm_IB3.tif'),
          layers=c('class','prob'))
writeTiles(modelRun,
           egTile,
           'Renderings/identity_reduced_', gisPath, '_IB3.tif',
           layers=c('class'))

# To mask of the water areas, i.e. to set them to some unique water value
data(water)
plot (egTile)
plot (mask(egTile,water,maskvalue=1,updatevalue=NA)) 
## End(Not run)

henkelstone/NPEL.Classification documentation built on May 17, 2019, 3:42 p.m.