专注系列化、高质量的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] TRUE
2 原点对齐
使用随机数给像元赋值:
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.1
4 翻转
翻转操作和栅格对齐关系不大,但作为平移的相关操作在这里介绍一下。翻转分为水平翻转(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)
![27e2bf43c0a985da27b1dfec3b844f5a.png](https://i-blog.csdnimg.cn/blog_migrate/a735df21040afd9954ea793f419fd3d0.png)
![2c035c06ec1557d1ae18024bdb335061.png](https://i-blog.csdnimg.cn/blog_migrate/2d1d0eb54575c5b138202d1a3567829b.png)
5 数值替换
在上面的情况中,偏移量小数点后的位数较少,可以通过平移“完美”地进行复位,但在实际情况中,偏移量小数点后往往“看不到尾”,因此很难精确地设置平移量,这时我们可以通过数值替换(replacement)来进行对齐。
r1 <- R1
r2 <- R2
r21 <- r1
values(r21) <- values(r2)
生成的
r21
的属性值来自r2
,相当于r2
的副本,而几何信息来自r1
,因此与r1
是对齐的。
6 重采样
上述对齐操作的实质都是通过平移来实现,其前提是栅格的行列数和分辨率相同。若这些信息中至少有一个不同时,就无法通过平移将其对齐,而应使用重采样操作。该操作在如下推文中已有介绍: