专注系列化、高质量的R语言教程
（本号已支持快捷转载，无需白名单即可转载）
栅格数据能不能对齐主要看它们如下几何信息是否保持一致：
地理/投影坐标（projection）
原点（origin）
范围（extent）
行、列数（nnumber of rows and columns）
分辨率（resolution）
如果上述信息存在不一致，栅格数据就会对不齐，从而影响后续操作。关于坐标系可参见如下几篇推文：
本篇主要介绍后四种信息不一致的解决方法。目录如下：
1 怎么判断栅格对不齐
2 原点对齐
3 平移
4 翻转
5 数值替换
6 重采样
1 怎么判断栅格对不齐
raster工具包的compareRaster()函数可以用来对比两个及以上的栅格数据是否有相同的范围、行列数、投影、分辨率和原点。它的语法结构如下：
compareRaster(x, ...,
extent = TRUE,
rowcol = TRUE,
crs = TRUE,
res = FALSE,
orig = FALSE,
rotation = TRUE,
values = FALSE,
tolerance,
stopiffalse = TRUE,
showwarning = FALSE)参数含义如下：
x、...：对比的栅格对象；extent、rowcol、crs、res、orig、rotation、values：分别表示是否把范围、行列数、坐标系、原点、旋转度、像元值纳入对比；tolerance：误差允许范围，数值在0-0.5之间，表示占单个像元大小的比例；stopiffalse：比较指标不同时是否报错；showwarning：是否显示警告信息；在stopiffalse = FALSE时有效。
定义如下两个栅格作为示例数据，它们的行列数和分辨率都相同，但范围略微有点错位：
library(raster)
R1 <- raster(xmn = 0, xmx = 10,
ymn = 0, ymx = 10,
resolution = 1)
R2 <- raster(xmn = 0.2, xmx = 10.2,
ymn = 0.2, ymx = 10.2,
resolution = 1)
R1
## class : RasterLayer
## dimensions : 10, 10, 100 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 10, 0, 10 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
R2
## class : RasterLayer
## dimensions : 10, 10, 100 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0.2, 10.2, 0.2, 10.2 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs保持其他参数为默认值：
compareRaster(R1, R2)
## Error in compareRaster(R1, R2) : different extent设置stopiffalse和showwarning参数：
compareRaster(R1, R2, stopiffalse = F)
## [1] FALSE
compareRaster(R1, R2, stopiffalse = F,
showwarning = T)
## Warning in compareRaster(R1, R2, stopiffalse = F, showwarning = T): different
## extent
## [1] FALSE设置tolerance参数：
compareRaster(R1, R2, tolerance = 0.2)
## [1] TRUE2 原点对齐
使用随机数给像元赋值：
set.seed(1024)
values(R1) <- rpois(100, 5)
values(R2) <- rpois(100, 2)
r1 <- R1
r2 <- R2然后对其进行几何运算：
r1+r2
## Error in compareRaster(e1, e2, extent = FALSE, rowcol = FALSE, crs = TRUE, :
## different origin
运行结果报错并提示栅格“原点不同”（different origin）。
原点是指最接近(0,0)的点。可以使用origin()函数查看栅格的原点：
origin(r1)
## [1] 0 0
origin(r2)
## [1] 0.2 0.2采用origin(x) <- value的形式可以修改栅格的原点，从而实现对齐：
extent(r2)
## class : Extent
## xmin : 0.2
## xmax : 10.2
## ymin : 0.2
## ymax : 10.2
origin(r2) <- origin(r1)
origin(r2)
## [1] 5.551115e-17 0.000000e+00
extent(r2)
## class : Extent
## xmin : 5.551115e-17
## xmax : 10
## ymin : 0
## ymax : 10
r1 + r2
## class : RasterLayer
## dimensions : 10, 10, 100 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 10, 0, 10 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2, 14 (min, max)
修改后的原点可能仍然会有极其微小的误差，但不影响后续操作；
修改原点的同时会改变栅格的范围，相当于平移。
3 平移
在上面的例子中，两个栅格的行列数和分辨率系统，只是范围有些错位；在实际中可能是由一些偶然因素导致的。如果知道偏移量，可以通过“平移”将其复位，使用的函数是shift()。
显然，r1和r2之间的左右方向和上下方向的偏移量都为0.2。假设造成的原因是r1左偏和下偏了0.1，r2右偏和上偏了0.1，那么可以通过将它们分别向右向上、向左向下平移0.1来复位：
r1 <- R1
r2 <- R2
r1 <- shift(r1, dx = 0.1, dy = 0.1)
r2 <- shift(r2, dx = -0.1, dy = -0.1)
extent(r1)
## class : Extent
## xmin : 0.1
## xmax : 10.1
## ymin : 0.1
## ymax : 10.1
extent(r2)
## class : Extent
## xmin : 0.1
## xmax : 10.1
## ymin : 0.1
## ymax : 10.14 翻转
翻转操作和栅格对齐关系不大，但作为平移的相关操作在这里介绍一下。翻转分为水平翻转（direction = "y"，默认情况）和垂直翻转（direction = "x"）：
r1 <- R1
r11 <- flip(r1, direction = "y")
r12 <- flip(r1, direction = "x")对比翻转前后的图象：
plot(r1)
par(mfrow = c(1,2))
plot(r11)
plot(r12)
5 数值替换
在上面的情况中，偏移量小数点后的位数较少，可以通过平移“完美”地进行复位，但在实际情况中，偏移量小数点后往往“看不到尾”，因此很难精确地设置平移量，这时我们可以通过数值替换（replacement）来进行对齐。
r1 <- R1
r2 <- R2
r21 <- r1
values(r21) <- values(r2)
生成的
r21的属性值来自r2，相当于r2的副本，而几何信息来自r1，因此与r1是对齐的。
6 重采样
上述对齐操作的实质都是通过平移来实现，其前提是栅格的行列数和分辨率相同。若这些信息中至少有一个不同时，就无法通过平移将其对齐，而应使用重采样操作。该操作在如下推文中已有介绍：
825
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
