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"))
![](https://i-blog.csdnimg.cn/blog_migrate/3bd10be1168f2c5c4dfc92026c0a0487.png)