R语言rasterImage r语言raster分块写入

这里小编推出一个系列来介绍R语言中管理空间栅格数据的工具包：raster。

实际上，在R中，管理栅格数据也有两个工具包，即raster和terra。与sf和sp不同的是，这两个工具包出自同一作者之手，后者可能会逐渐取代前者作为R语言中管理栅格数据的基础工具包。

加载工具包：

library(raster)

栅格对象

从外界读取栅格数据

使用的函数是与工具包同名的函数raster()，语法结构如下：

raster(filename, layer = 0)

• filename：数据的存储位置；
• layer：图层。

system.file("external/test.grd", package="raster")
## "D:/R-3.6.3/library/raster/external/test.grd"

r <- raster("D:/R-3.6.3/library/raster/external/test.grd")

r的数据结构如下：

• 可以看出，raster工具包的栅格对象与sp工具包中的矢量对象的数据结构是类似的。

r
## class      : RasterLayer 
## dimensions : 115, 80, 9200  (nrow, ncol, ncell)
## resolution : 40, 40  (x, y)
## extent     : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
## crs        : +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 
## source     : D:/R-3.6.3/library/raster/external/test.grd 
## names      : test 
## values     : 138.7071, 1736.058  (min, max)

• class：数据结构类型；
• dimensions：维度，即行数、列数、空间单元数（nrow、ncol、ncell）；
• resolution：分辨率，即单个空间单元的宽（x）和高（y）;
• extent：数据范围，即x和y方向的最小值和最大值（xmin、xmax、ymin、 ymax）；
• crs：投影信息；
• source：导入数据的位置；
• names：数据名称；
• values：空间单元属性值的范围（min、max）。

以上参数之间存在数量关系：

undefined

创建栅格数据

创建栅格数据的函数也是raster()，默认情况下创建的是全球范围1*1经纬度的栅格对象：

x1 <- raster()
x1
## class      : RasterLayer 
## dimensions : 180, 360, 64800  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs

完整的语法结构如下：

raster(nrows = 180, ncols = 360,
       xmn = -180, xmx = 180, ymn = -90, ymx = 90,
       crs, ext, resolution, vals = NULL)

由于分辨率、维度和数据范围之间存在数量关系，因此任意指定其中两者就可以确定栅格数据的形态。

指定维度与数据范围

x2 <- raster(nrows = 18, ncols = 36, 
             xmn = -900, xmx = 900,
             ymn = -900, ymx = 900)
x2
## class      : RasterLayer 
## dimensions : 18, 36, 648  (nrow, ncol, ncell)
## resolution : 50, 100  (x, y)
## extent     : -900, 900, -900, 900  (xmin, xmax, ymin, ymax)
## crs        : NA

指定数据范围与分辨率

x3 <- raster(xmn = -900, xmx = 900,
             ymn = -900, ymx = 900,
             resolution = c(50, 100))
x3
## class      : RasterLayer 
## dimensions : 18, 36, 648  (nrow, ncol, ncell)
## resolution : 50, 100  (x, y)
## extent     : -900, 900, -900, 900  (xmin, xmax, ymin, ymax)
## crs        : NA

指定维度与分辨率

x4 <- raster(nrows = 18, ncols = 36, 
             resolution = c(50, 100))
x4
## class      : RasterLayer 
## dimensions : 18, 36, 648  (nrow, ncol, ncell)
## resolution : 50, 100  (x, y)
## extent     : -900, 900, -900, 900  (xmin, xmax, ymin, ymax)
## crs        : NA

在指定数据范围时，如果想和已知的栅格数据的范围保持一致，可以调用ext参数：

x5 <- raster(nrows = 18, ncols = 36, 
             ext = extent(x2))
x5
## class      : RasterLayer 
## dimensions : 18, 36, 648  (nrow, ncol, ncell)
## resolution : 50, 100  (x, y)
## extent     : -900, 900, -900, 900  (xmin, xmax, ymin, ymax)
## crs        : NA

• extent()函数用于提取栅格对象的数据范围；
• ext = extent(x2)表示数据范围与x2一致。

多图层栅格数据

栅格数据可以包含多个图层。stack()函数可以将多个栅格数据合并成一个栅格数据的多个图层：

r1 <- raster(ncols = 3, nrows = 3,
             vals = c(1:9))
r2 <- raster(ncols = 3, nrows = 3,
             vals = rep(c(0,1,2), 3))
r <- stack(r1, r2)
r
## class      : RasterStack 
## dimensions : 3, 3, 9, 2  (nrow, ncol, ncell, nlayers)
## resolution : 120, 60  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : layer.1, layer.2 
## min values :       1,       0 
## max values :       9,       2

plot(r)

undefined

代数运算

栅格数据之间可以进行代数运算。

r1 <- raster(ncols = 3, nrows = 3,
             vals = c(1:9))
r2 <- raster(ncols = 3, nrows = 3,
             vals = rep(c(0,1,2), 3))

getValues()函数可以查询栅格对象各单元的属性值：

getValues(r1)
## [1] 1 2 3 4 5 6 7 8 9
getValues(r2)
## [1] 0 1 2 0 1 2 0 1 2

四则运算

加减乘除、开方等运算都可以适用于栅格数据：

(a) 单个栅格数据的运算

getValues(r1+10)
## [1] 11 12 13 14 15 16 17 18 19
getValues(sqrt(r1))
## [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427
## [9] 3.000000

(b) 多个栅格数据间的运算

getValues(r1 + r2)
## [1]  1  3  5  4  6  8  7  9 11
getValues(r1 * r2)
## [1]  0  2  6  0  5 12  0  8 18

逻辑运算

>、==、!等逻辑运算符也可适用于栅格数据：

plot(r2 == 1)

• r2 == 1表示筛选出r2中属性值为1的单元。

汇总运算

平均值、极值等汇总函数也适用于栅格数据：

getValues(mean(r1, r2))
## [1] 0.5 1.5 2.5 2.0 3.0 4.0 3.5 4.5 5.5

可视化

基础绘图系统中的plot()函数经过再封装可以用于栅格数据的可视化，具体可在raster工具包中查看plot()函数的用法，而graphics包中的原有参数仍然可以使用，但部分默认值发生了变化。

filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)
plot(r)

此前，本号也介绍了一些graphics工具包中与地理分析相关的绘图函数，如contour()、persp()等，这些函数也可以用于栅格对象的可视化：

par(mfrow = c(2,2), plt = c(0.1,0.95,0.15,0.9))
contour(r)
persp(r)
hist(r)
density(r)