【R语言】使用r语言进行单个或多个栅格数据间的计算

写在前面

本文介绍如何使用r语言进行单个或多个栅格数据间的计算，类似于ArcGIS中的栅格计算器。

单个栅格数据：

library(raster)
calc(x, fun, filename='', na.rm, forcefun=FALSE, forceapply=FALSE, ..)

元素：

x：Raster* object
fun：function
filename：character. Output filename (optional)
Examples
# NOT RUN {
r <- raster(ncols=36, nrows=18)
values(r) <- 1:ncell(r)

# multiply values with 10
fun <- function(x) { x * 10 }
rc1 <- calc(r, fun)

# set values below 100 to NA. 
fun <- function(x) { x[x<100] <- NA; return(x) }
rc2 <- calc(r, fun)

# set NA values to -9999
fun <- function(x) { x[is.na(x)] <- -9999; return(x)} 
rc3 <- calc(rc2, fun)

# using a RasterStack as input
s <- stack(r, r*2, sqrt(r))
# return a RasterLayer
rs1 <- calc(s, sum)

# return a RasterBrick
rs2 <- calc(s, fun=function(x){x * 10})
# recycling by layer
rs3 <- calc(s, fun=function(x){x * c(1, 5, 10)})

# use overlay when you want to refer to individual layer in the function
# but it can be done with calc: 
rs4 <- calc(s, fun=function(x){x[1]+x[2]*x[3]})

## 
# Some regression examples
## 

# create data
r <- raster(nrow=10, ncol=10)
s1 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))
s2 <- lapply(1:12, function(i) setValues(r, rnorm(ncell(r), i, 3)))
s1 <- stack(s1)
s2 <- stack(s2)

# regression of values in one brick (or stack) with another
s <- stack(s1, s2)
# s1 and s2 have 12 layers; coefficients[2] is the slope
fun <- function(x) { lm(x[1:12] ~ x[13:24])$coefficients[2] }
x1 <- calc(s, fun)

# regression of values in one brick (or stack) with 'time'
time <- 1:nlayers(s)
fun <- function(x) { lm(x ~ time)$coefficients[2] }
x2 <- calc(s, fun)

# get multiple layers, e.g. the slope _and_ intercept
fun <- function(x) { lm(x ~ time)$coefficients }
x3 <- calc(s, fun)


### A much (> 100 times) faster approach is to directly use 
### linear algebra and pre-compute some constants

## add 1 for a model with an intercept
X <- cbind(1, time)

## pre-computing constant part of least squares
invXtX <- solve(t(X) %*% X) %*% t(X)

## much reduced regression model; [2] is to get the slope
quickfun <- function(y) (invXtX %*% y)[2]
x4 <- calc(s, quickfun) 
# }

参考：https://www.rdocumentation.org/packages/raster/versions/3.1-5/topics/calc

---

多个栅格数据：

overlay(x, y, .., fun, datatype, format, overwrite, progress)
1) x and y are Raster* objects
 2) x is a Raster* object, y is missing (this is equivalent to calc)

Examples
r <- raster(ncol=10, nrow=10)
r1 <- init(r, fun=runif)
r2 <- init(r, fun=runif)
r3 <- overlay(r1, r2, fun=function(x,y){return(x+y)})

# long version for multiplication
r4 <- overlay(r1, r2, fun=function(x,y){return(x*y)} )

#use a stack
s <- stack(r1, r2)
r5 <- overlay(s, fun=function(x,y){return(x*y)} )

# use a single RasterLayer (same as calc function)
r6 <- overlay(r1, fun=function(x){return(sqrt(x))} )

# multiplication with more than two layers (make sure the number of RasterLayers matches the arguments of 'fun'
r7 <- overlay(r1, r2, r3, r4, fun=function(a,b,c,d){return(a*b+c*d)} )  
# equivalent function, efficient if values can be loaded in memory
r8 <- r1 * r2 + r3 * r4

参考：https://www.rdocumentation.org/packages/raster/versions/1.1.7/topics/overlay