fixes #948

DESCRIPTION

@@ -2,7 +2,7 @@ Package: terra
 Type: Package
 Title: Spatial Data Analysis
 Version: 1.6-49
-Date: 2022-12-12
+Date: 2022-12-15
 Depends: R (>= 3.5.0)
 Suggests: parallel, tinytest, ncdf4, sf (>= 0.9-8), deldir, XML, leaflet
 LinkingTo: Rcpp

NEWS.md

@@ -1,5 +1,8 @@
-# version 1.6-48
+# version 1.6-49
 
+## enhancements
+
+- `shade` is now vectorized for arguments `angle` and `direction` to facilitate generating multiple hillshades that can be combined for a better result [#948](https://github.com/rspatial/terra/issues/948) by Jürgen Niedballa
 
 ## bug fixes 
 

R/plot.R

@@ -300,23 +300,10 @@ setMethod("barplot", "SpatRaster",
 
 
 
-
 shade <- function(slope, aspect, angle=45, direction=0, normalize=FALSE, filename="", overwrite=FALSE, ...) {
-
-	x <- c(slope[[1]], aspect[[1]])
-
-	direction <- direction[1] * pi/180
-	zenith <- (90 - angle[1]) * pi/180
-
-	if (normalize) {
-		fun <- function(slp, asp) {
-			shade <- cos(slp) * cos(zenith) + sin(slp) * sin(zenith) * cos(direction-asp)
-			shade[shade < 0] <- 0
-			shade * 255
-		}
-	} else {
-		fun <- function(slp, asp) { cos(slp) * cos(zenith) + sin(slp) * sin(zenith) * cos(direction-asp) }
-	}
-	lapp(x, fun=fun, filename=filename, overwrite=overwrite, wopt=list(...))
+	stopifnot(inherits(slope, "SpatRaster"))
+	opt <- terra:::spatOptions(filename, overwrite=overwrite, ...)
+	slope@ptr <- slope@ptr$hillshade(aspect@ptr, angle, direction, normalize[1], opt)
+	messages(slope, "shade")
 }
 

man/shade.Rd

@@ -18,8 +18,8 @@ shade(slope, aspect, angle=45, direction=0, normalize=FALSE,
 \arguments{
   \item{slope}{SpatRasterwith slope values (in radians) }
   \item{aspect}{SpatRaster with aspect values (in radians) }
-  \item{angle}{ The the elevation angle of the light source (sun), in degrees}
-  \item{direction}{ The direction (azimuth) angle of the light source (sun), in degrees}
+  \item{angle}{ The the elevation angle(s) of the light source (sun), in degrees}
+  \item{direction}{ The direction (azimuth) angle(s) of the light source (sun), in degrees}
   \item{normalize}{Logical. If \code{TRUE}, values below zero are set to zero and the results are multiplied with 255}
   \item{filename}{character. Output filename}
   \item{overwrite}{logical. If \code{TRUE}, \code{filename} is overwritten}  
@@ -42,6 +42,14 @@ aspect <- terrain(alt, "aspect", unit="radians")
 hill <- shade(slope, aspect, 40, 270)
 plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
 plot(alt, col=rainbow(25, alpha=0.35), add=TRUE)
+
+
+# A better hill shade may be achieved by combining 
+# different angles and directions. For example
+
+h <- shade(slope, aspect, angle = c(45, 45, 45, 80), direction = c(225, 270, 315, 135))
+h <- Reduce(mean, h)				
+
 }
 
 

src/RcppModule.cpp

@@ -891,6 +891,7 @@ RCPP_MODULE(spat){
 		.method("scale", &SpatRaster::scale)
 		.method("shift", &SpatRaster::shift)
 		.method("terrain", &SpatRaster::terrain)
+		.method("hillshade", &SpatRaster::hillshade)
 		.method("summary", &SpatRaster::summary)
 		.method("summary_numb", &SpatRaster::summary_numb)
 		.method("transpose", &SpatRaster::transpose)

src/distRaster.cpp

@@ -24,6 +24,7 @@
 #include "math_utils.h"
 #include "vecmath.h"
 #include "file_utils.h"
+#include "string_utils.h"
 #include "crs.h"
 #include "sort.h"
 
@@ -3623,3 +3624,102 @@ SpatRaster SpatRaster::terrain(std::vector<std::string> v, unsigned neighbors, b
 	return out;
 }
 
+
+
+
+SpatRaster SpatRaster::hillshade(SpatRaster aspect, std::vector<double> angle, std::vector<double> direction, bool normalize, SpatOptions &opt) {
+
+	SpatRaster out = geometry(1);
+	if ((nlyr() != 1) || (aspect.nlyr() != 1)) {
+		out.setError("slope and aspect should have one layer");
+		return out;
+	}
+	if ((angle.size() == 0) || (direction.size() == 0)) {
+		out.setError("you must provide a value for aspect and direction");
+		return out;		
+	}
+
+	std::vector<std::string> nms;
+
+	if ((angle.size() > 1) || (direction.size() > 1)) {
+		recycle(angle, direction);
+		recycle(direction, angle);	
+		//nms = opt.names;
+		SpatOptions ops(opt);
+		size_t nl = angle.size();
+		out.source.resize(nl);
+		if (ops.names.size() == nl) {
+			nms = opt.names;
+		} else {
+			nms.reserve(nl);
+			for (unsigned i=0; i<nl; i++) {
+				std::string nmi = "hs_" + double_to_string(angle[i]) + "_" + double_to_string(direction[i]);
+				nms.push_back(nmi);
+			}
+		}
+
+		for (unsigned i=0; i<nl; i++) {
+			ops.names = {nms[i]};
+			SpatRaster r = hillshade(aspect, {angle[i]}, {direction[i]}, normalize, ops);
+			out.source[i] = r.source[0];
+		}
+		if (opt.get_filename() != "") {
+			out = out.writeRaster(opt);
+		}
+		return out;
+	} else {
+		if ((opt.names.size() == 1) && (opt.names[0] != "")) {
+			out.setNames(  {opt.names[0] } );
+		} else {
+			out.setNames( {"hillshade"} );
+		}
+	}
+
+	double dir = direction[0] * M_PI / 180.0;
+	double zen = (90.0 - angle[0]) * M_PI/180.0;
+	double coszen = cos(zen);
+	double sinzen = sin(zen);
+
+	if (!readStart()) {
+		out.setError(getError());
+		return(out);
+	}
+	if (!aspect.readStart()) {
+		out.setError(getError());
+		return(out);
+	}
+
+  	if (!out.writeStart(opt, filenames())) {
+		readStop();
+		return out;
+	}
+
+	for (size_t i = 0; i < out.bs.n; i++) {
+		std::vector<double> slp;
+		std::vector<double> asp;
+		readBlock(slp, out.bs, i);
+		aspect.readBlock(asp, out.bs, i);
+		if (normalize) {
+			for (size_t i=0; i<slp.size(); i++) {
+				slp[i] = cos(slp[i]) * coszen + sin(slp[i]) * sinzen * cos(dir-asp[i]);
+				if (slp[i] < 0) {
+					slp[i] = 0;
+				} else {
+					slp[i] *= 255;
+				}
+			}
+		} else {
+			for (size_t i=0; i<slp.size(); i++) {
+				slp[i] = cos(slp[i]) * coszen + sin(slp[i]) * sinzen * cos(dir-asp[i]);
+			}
+		}
+		if (!out.writeBlock(slp, i)) return out;
+	}
+	out.writeStop();
+	readStop();
+	aspect.readStop();
+	return out;
+}
+
+
+

src/spatRaster.h

@@ -740,6 +740,7 @@ class SpatRaster {
 
 		SpatRaster scale(std::vector<double> center, bool docenter, std::vector<double> scale, bool doscale, SpatOptions &opt);
 		SpatRaster terrain(std::vector<std::string> v, unsigned neighbors, bool degrees, unsigned seed, SpatOptions &opt);
+		SpatRaster hillshade(SpatRaster aspect, std::vector<double> angle, std::vector<double> direction, bool normalize, SpatOptions &opt);
 
 		SpatRaster selRange(SpatRaster x, int z, int recycleby, SpatOptions &opt);
 		SpatRaster selectHighest(size_t n, bool low, SpatOptions &opt);