raster | R语言中的空间栅格对象及其基本处理方法（Ⅱ）：投影、属性提取-CSDN博客

本文链接：https://blog.csdn.net/weixin_54000907/article/details/119091651

library(raster)

投影

设置投影

在使用raster()函数创建栅格对象时，调用crs参数可以定义其投影。投影参数具有默认值：

r <- raster()
crs(r)
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

类似于+proj=longlat +datum=WGS84 +no_defs的文本就是描述投影信息的WKT（Well-known text）；关于WKT可查询网站http://geotiff.maptools.org/proj_list/。

crs()函数既可以读取栅格对象的投影信息，也可以定义其投影：

# 读取投影
crs(x, asText=FALSE, ...)
# 定义投影
crs(x, ...) <- value

crs(r) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84"
crs(r)
## CRS arguments:
##  +proj=lcc +lat_0=0 +lon_0=-100 +lat_1=48 +lat_2=33 +x_0=0 +y_0=0
## +datum=WGS84 +units=m +no_defs

projection()函数的功能与crs()函数一致。

转换投影

投影转换是将栅格对象的投影由一种转换到另一种，使用的函数是projectRaster()。语法结构如下：

projectRaster(from, to, res, crs,
              method="bilinear", alignOnly=FALSE,
              over=FALSE, filename="", ...)

from：需转换投影的栅格对象。

projectRaster()函数有两种方式进行投影：

指定包含目标投影方式的栅格对象作为模板，即参数to；
指定目标投影方式的WKT文本，即参数crs。

r_2 <- raster("D:/R-3.6.3/library/raster/external/test.grd")
crs(r_2)
## CRS arguments:
##  +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
## +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs

r_3 <- projectRaster(r, r_2)
crs(r_3)
## CRS arguments:
##  +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
## +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs

属性提取

属性提取操作是将一个空间对象（记为x）的属性信息映射到另外一个空间对象（记为y）上，使用的函数是extract()。

extract(x, y, ...)

矢量对象y既可以是sp对象，也可以是sf对象。

根据矢量对象y的数据类型有如下几种情况：

y为点要素

返回每个点要素所在栅格的属性值，或者缓冲区的栅格对象的属性汇总值，此时extract()函数的语法结构如下：

extract(x, y, method = 'simple',
        buffer = NULL, small = FALSE,
        cellnumbers = FALSE, fun = NULL,
        na.rm = TRUE, layer, nl,
        df = FALSE, factors = FALSE, ...)

method：两个可选项；simple表示直接返回所在栅格的属性值，bilinear表示返回最近的4个栅格的插值结果；
buffer：数值型参数，为每个点要素指定缓冲区范围；当使用缓冲区时，其他参数的调用方式同线要素和面要素的情况。

y为线要素

由于线要素可能与栅格对象的多个要素有交点，返回的是相交栅格的属性汇总值，此时extract()函数的语法结构如下：

extract(x, y, fun = NULL, na.rm = FALSE,
        cellnumbers = FALSE, df = FALSE,
        layer, nl, factors = FALSE,
        along = FALSE, sp=FALSE, ...)

fun：汇总函数名；如mean表示取相交栅格属性的平均值；
cellnumbers：是否返回相交栅格的个数；
df：结果是否以数据框的形式返回。

y为面要素

返回中心在面要素内的栅格的属性汇总值：

extract(x, y, fun = NULL, na.rm = FALSE,
        weights = FALSE, normalizeWeights = TRUE,
        cellnumbers = FALSE, small = TRUE,
        df = FALSE, layer, nl, factors = FALSE,
        sp=FALSE, ...)

weights：属性汇总时是否以相交的面积作为权重，默认为FALSE；
normalizeWeights：当使用面积加权时权重是否标准化，默认为TRUE；
small：对于面积较小的面要素的处理方式。当没有任何一个栅格的中心在该面要素内，small = TRUE且面要素落在一个栅格内，则返回该栅格的属性值；否则返回空值。

library(sf)
library(raster)
# 2015年中国人口空间分布公里网格数据集
# 数据来源：https://www.resdc.cn/DOI/DOI.aspx?DOIid=32
pop <- raster("60-1.pop2015/tpop2015/w001001.adf")
beijing <- read_sf("60-2.ChinaShape/县.shp") %>%
        dplyr::filter(省 == "北京市")

beijing$pop <- extract(pop, beijing, fun = sum,
                       df = T, na.rm = T)[,2]
beijing$pop
##  [1]  826160.1 1322401.1 3841823.7 2255492.6  646506.4 3682926.9  469271.1
##  [8] 1058469.6 1265216.9 1054529.7 1936018.9 1594486.1  429912.2  460520.9
## [15]  510374.8  337011.8

plot(beijing["pop"], breaks = "quantile",
     nbreaks = 5, pal = hcl.colors(5, "Blue-Red"))