R/terrain_unix.R

In hrvg/RiverML: Machine Learning for River Science

#' Modified `terrain()` function
#' 
#' This modified version of `raster::terrain()` can compute curvature. See details.
#' 
#' Author: Robert J. Hijmans
#' Date : February 2011
#' Version 1.0
#' Licence GPL v3
#' Modified by Herve Guillon, November 2017
#' Added options for computation of planform and profile curvature
#' @param x raster
#' @param opt type of metrics
#' @param unit radians or tangent
#' @param neighbors default to queen case
#' @param filename passed to WriteRaster
#' @return a raster with the computed metrics
#' @importFrom raster blockSize brick canProcessInMemory couldBeLonLat getValues hasValues isLonLat nlayers pbClose pbCreate pbStep pointDistance projection raster res setValues trim values writeRaster writeStart writeStop writeValues yFromRow
#' @useDynLib RiverML _do_terrain_
#' @export
#' @keywords tam-dm
terrain_ <- function(x, opt='slope', unit='radians', neighbors=8, filename='', ...) {

	

	if (nlayers(x) > 1) {

		warning('first layer of x is used')

		x <- subset(x, 1)

	}

	stopifnot(hasValues(x))

	

	stopifnot(is.character(filename))

	filename <- trim(filename)

 

	stopifnot(is.character(opt))

	opt <- unique(trim(tolower(opt)))

	i <- which(! opt %in% c('tri', 'tpi', 'roughness','slope', 'aspect', 'flowdir', 'curvplan', 'curvprof'))

	if (length(i) > 0) {

		stop('invalid value in "opt", choose from:\n "tri", "tpi", "roughness", "slope", "aspect", "flowdir", "curvplan", "curvprof"')

	}

	stopifnot(length(opt) > 0 ) 

	

 

	nopt <- rep(0, 10)

	if ('tri' %in% opt) {

		nopt[1] <- 1

	} 

	if ('tpi' %in% opt) {

		nopt[2] <- 1

	} 

	if ('roughness' %in% opt) {

		nopt[3] <- 1

	}

	if ('slope' %in% opt) {

		if (neighbors == 4) {

			nopt[4] <- 1

		} else {

			nopt[6] <- 1

		}

	}

	if ('aspect' %in% opt) {

		if (neighbors == 4) {

			nopt[5] <- 1

		} else {

			nopt[7] <- 1

		}

	} 

	if ('flowdir' %in% opt) {

		nopt[8] <- 1

	}

	if ('curvplan' %in% opt) {

		nopt[9] <- 1

	}

	if ('curvprof' %in% opt) {

		nopt[10] <- 1

	}
	

	nopt <- as.integer(nopt)

	nl <- sum(nopt)

	

	if (nl == 1) {

		out <- raster(x)

	} else {

		out <- brick(x, values=FALSE, nl=nl)

	}

	names(out) <- c('tri', 'tpi', 'roughness','slope', 'aspect', 'slope', 'aspect', 'flowdir', 'curvplan', 'curvprof')[as.logical(nopt)]

 

	rs <- as.double(res(out))

	un <- as.integer(1)

	lonlat <- FALSE

	if ('slope' %in% opt | 'aspect' %in% opt | 'flowdir' %in% opt | 'curvplan' %in% opt | 'curvprof' %in% opt) {

		stopifnot(is.character(unit))

		unit <- trim(tolower(unit))

		stopifnot(unit %in% c('degrees', 'radians', 'tangent'))

		if (unit=='degrees') {

			un <- as.integer(0)

		} else if (unit=='tangent') {

			un <- as.integer(2)		

		}

		stopifnot(neighbors %in% c(4, 8))

		stopifnot(! is.na(projection(x)) )

		lonlat <- isLonLat(out)

		if (!lonlat & couldBeLonLat(out)) {

			warning('assuming CRS is longitude/latitude')

			lonlat <- TRUE

		}

 

		if (lonlat) {		

			rs[2] <- pointDistance(cbind(0,0), cbind(0, rs[2]), longlat=TRUE)

		}

	}

	lonlat <- as.integer(lonlat)

	
	if (canProcessInMemory(out)) {

		if (lonlat) {

			y <- yFromRow(x, 1:nrow(x))

		} else {

			y <- 0

		}

		v <- .Call('_do_terrain_', as.double(values(x)), as.integer(dim(out)), rs, un, nopt, lonlat, y, NAOK=TRUE)

		out <- setValues(out, v)

		if (filename  != '') {

			out <- writeRaster(out, filename, ...)

		}

		

	} else {

 

		out <- writeStart(out, filename, ...)

		tr <- blockSize(out, minblocks=3, minrows=3)

		pb <- pbCreate(tr$n, label='terrain', ...)

		

		nc <- ncol(out)

		buf <- 1:nc

		v <- getValues(x, row=1, nrows=tr$nrows[1]+1)

		y <- 0

		if (lonlat) { y <- yFromRow(out, 1:(tr$nrows[1]+1)) }
		
		v <- .Call('_do_terrain_', as.double(v), as.integer(c(tr$nrows[1]+1, nc)), rs, un, nopt, lonlat, y, NAOK=TRUE)

		out <- writeValues(out, matrix(v, ncol=nl), 1)

		pbStep(pb, 1)

		for (i in 2:(tr$n-1)) {

			v <- getValues(x, row=tr$row[i]-1, nrows=tr$nrows[i]+2)

			if (lonlat) { y <- yFromRow(out, (tr$row[i]-1) : (tr$row[i]+tr$nrows[i])) }

			v <- .Call('_do_terrain_', as.double(v), as.integer(c(tr$nrows[i]+2, nc)), rs, un, nopt, lonlat, y, NAOK=TRUE)

			v <- matrix(v, ncol=nl)[-buf,]

			out <- writeValues(out, v, tr$row[i])

			pbStep(pb, i)		

		}

		i <- tr$n

		v <- getValues(x, row=tr$row[i]-1, nrows=tr$nrows[i]+1)

		if (lonlat) { y <- yFromRow(out, (tr$row[i]-1) : (tr$row[i]+tr$nrows[i]-1)) }

		v <- .Call('_do_terrain_', as.double(v), as.integer(c(tr$nrows[i]+1, nc)), rs, un, nopt, lonlat, y, NAOK=TRUE)

		v <- matrix(v, ncol=nl)[-buf,]

		out <- writeValues(out, v, tr$row[i])

		pbStep(pb, tr$n)

		

		out <- writeStop(out)

		pbClose(pb)

	}

	return(out)

}

 

 

# x <- terrain(utm, out='tri')