raster | R语言中的空间栅格对象及其基本处理方法（Ⅲ）：切片/掩膜、图层叠加...

本篇介绍如下处理方法：

• 根据矢量数据的范围对栅格数据进行裁剪；

• 通过对栅格数据进行运算得到新的栅格对象，如通过人口和GDP的栅格数据得到人均GDP的栅格数据。

示例的栅格数据是2015年中国人口和GDP的公里网格数据，各位读者可以使用下方提供的链接自行下载，也可以在公众号后台发送关键词“示例数据”获取。

library(raster)
# 2015年中国人口空间分布公里网格数据集
# 数据来源：https://www.resdc.cn/DOI/DOI.aspx?DOIid=32
pop <- raster("60-1.pop2015/tpop2015/w001001.adf")

# 2015中国GDP空间分布公里网格数据集
# 数据来源：https://www.resdc.cn/DOI/doi.aspx?DOIid=33
gdp <- raster("61-1.gdp2015/gdp2015/w001001.adf")

切片/掩膜

当需要使用矢量数据对栅格数据进行裁剪时，可以使用切片或掩膜操作。

切片函数crop()使用矢量数据的矩形范围对栅格数据进行裁剪。函数语法结构如下：

crop(x, y, filename = "", snap = 'near',
     datatype = NULL, ...)

• x：待裁剪的栅格对象；

• y：裁剪范围或包含裁剪范围的空间数据对象（矢量、栅格数据均可）；

• filename：裁剪后的数据储存地址，可选；

• snap：数据对齐方式，可选项有near、in、out。

掩膜函数mask()保持矢量数据实际范围内的属性数据不变，而将范围外的属性数据都转为其他值（默认转换为空值NA）。函数语法结构如下：

mask(x, mask, filename = "", inverse = FALSE, 
     maskvalue = NA, updatevalue = NA,
     updateNA = FALSE, ...)

• x：栅格对象；

• mask：作为掩膜的空间数据对象（矢量、栅格数据均可）；

• inverse：是否执行反向掩膜操作；

• maskvalue：掩膜内的栅格单元的属性值若为该值则将其转换为updatevalue参数对应的数值；

• updatevalue：掩膜外的所有栅格单元的属性值（空值除外）会被转换为该值，默认为空值NA；

• updateNA：若为TRUE，则掩膜外的空值也会被转换为updatevalue参数对应的数值。

示例矢量数据：

library(sf)
# 北京市县级行政区划矢量数据
bj <- read_sf("60-2.ChinaShape/县.shp") %>%
  dplyr::filter(省 == "北京市")

切片与掩膜操作要求两个对象的投影信息一致，否则会报错。如下：

bjpop <- crop(pop, bj)
## Error in .local(x, y, ...) : extents do not overlap

比较切片与掩膜操作的效果：

# 转换投影
bj_proj <- st_transform(bj, crs(pop))
# 切片
bj_crop <- crop(pop, bj_proj)
# 掩膜
bj_mask <- mask(pop, bj_proj)
# 可视化
plot(bj_crop, main = "crop")
plot(bj_mask, main = "mask")

• 切片会改变原有栅格对象的范围，但其得到的是一个矩形范围；

• 掩膜不会改变原有栅格对象的范围，仅是将矢量范围外的数据变为空值，因此右图看起来北京的范围很小。

对切片后的栅格对象再进行掩膜操作，可以得到预期裁剪效果：

bj_pop <- mask(crop(pop, bj_proj), bj_proj)
plot(bj_pop)

调用mask()函数的inverse参数可进行反向掩膜操作：

bj_pop2 <- mask(crop(pop, bj_proj), bj_proj,
                inverse = T)
plot(bj_pop2)

图层叠加

本部分的示例数据是北京市范围内的人口和GDP栅格数据：

bj_pop <- mask(crop(pop, bj_proj), bj_proj)
bj_gdp <- mask(crop(gdp, bj_proj), bj_proj)

如果想通过叠加人口和GDP的栅格数据得到人均GDP的栅格数据，有如下两种方法：

• 第一种方法在本系列的第一篇推文已经进行了介绍，即直接使用四则运算：

bj_pergdp <- bj_gdp/bj_pop

• 第二种方法是使用overlay()函数，这种方法也容易理解。语法结构如下：

overlay(x, y, ..., fun, filename = "",
        recycle = TRUE, forcefun = FALSE)

• x、y、...：参与叠加的栅格对象；

• fun：进行属性数据运算的函数。

bj_pergdp2 <- overlay(bj_gdp, bj_pop, 
                      fun = function(x, y) x/y)

以上两种方法得到的结果是一样的：

# 结果可视化
par(mfrow = c(1,2))
plot(bj_pergdp)
plot(bj_pergdp2)

当只对单个栅格对象进行运算时，可以使用calc()函数；cover()函数会将它的第一个参数中的空值使用第二个参数对象位置的属性值进行替代：

bj_calc <- calc(bj_pop, fun = function(x) return(30000))
bj_cover <- cover(bj_pop, bj_calc)
# 结果可视化
par(mfrow = c(1,2))
plot(bj_calc, main = "calc")
plot(bj_cover, main = "cover")

下期预告

图层叠加操作要求这些栅格对象之间的原点、分辨率和范围必须一致，本系列的下篇推文将介绍在上述指标不一致时如何对其进行调整。