R语言实践——chronosphere入门（古地理图+古地理坐标重建）

目录

chronosphere简介

功能结构

1. RasterArray

 2. SpatialArray

 访问RasterArray和SpatialArray的属性和功能

构建子集和重组

1. 挑选图层（layer）——双方括号

 2. 单方括号

从raster和sp中继承的功能

 1. RasterArray

 2. SpatialArray

绘图

调色板

RasterLayer

 3. RasterArray

添加一个统一的图例

 ​编辑

RasterArray处理NA图层

列数和页数

---

chronosphere简介

chronosphere目的是为了更好地处理rasters，并提供了GPlates的功能（包括rgplates）

功能结构

1. RasterArray

直译就是栅格组。它就是一组RasterLayer，可以像操作R语言的数组一样对其进行切割、组合等操作。
RasterArray类有两部分：@stack和@index。stack包含了栅格数据，该部分就是一个RasterStack对象。index描述RasterArray的结构。

以下内容参照示例数据dems

library(chronosphere)

data("dems")
dems

 第一部分展示RasterLayer共有的属性：维度（dimensions），分辨率（resolution）以及坐标系统（CRS，即extent和coord. ref.）。

第二部分展示RasterArray结构（Array properties）。

使用proxy（dems）可以展示RasterArray的第二部分数据，即index内容。

proxy(dems)

 从已存在的RasterArray对象中提取数据并生成新的RasterArray对象：

library(chronosphere)

data("dems")
dems

proxy(dems)

# 准备stack栅格数据
stackOfLayers = dems@stack

#准备index数据
ind = 1:10
names(ind) = letters[1:10]

#生成新的RasterArray对象
nra = RasterArray(index=ind, stack=stackOfLayers)
nra

 index对象决定了RasterArray对象的结构。

还可以使用combine（）函数将多个RasterLayer或RasterArray对象构建成RasterArray：

#新建两个相同结构的栅格对象r1和r2
r1 = raster()
values(r1) = 1
r2 = raster()
values(r2) = 2

#combine构建RasterArray对象
comb = combine(r1, r2)
comb

 还可以通过cbind（）和rbind（）函数构建类似矩阵的RasterArray对象：

cbind(dems, dems)

 2. SpatialArray

该类同样具有@index和@stack，其与RasterArray不同之处在于，它的@stack中有一个SpatialStack对象，由Spatial*对象组成的列表。

SpatialStack对象可以是SpatialPoints、SpatialPointDataFrame、SpatialLines、SpatialLineDataFrame、SpatialPolygons或SapatialPolygonsDataFrame任意对象，只要求它们的坐标系统相同即可。

以下示例参照coasts

data(coasts)
coasts

 这是一个二维SpatialArray对象，有5行2列：分别是大陆边缘和海岸线。proxy（）函数可以很容易地访问并分析该对象。

proxy(coasts)

 同样也可以用combine（）将数个Spatial*对象生成为新的SpatialArray对象：

margin1 = coasts[1, 1]
margin2 = coasts[2, 1]
combination = combine(margin1, margin2)
combination

当然也可以通过更底层的操作完成，首先要新建一个SpatialStach对象：

coast0 = coasts[1, 2]
coast5 = coasts[1, 2]
spStack = stack(margin0, margin5, coast0, coast5)

 然后用合适的index和spStack合并起来：

ind = matrix(1:4, ncol = 2)
colnames(ind) = c("margin", "coast")
rownames(ind) = c("0", "5")
spArray = SpatialArray(index=ind, stack=spStack)
spArray

 访问RasterArray和SpatialArray的属性和功能

访问并更改RasterArray属性的方法与Raster*对象更相似，而不是数组。它们都连接到@index，RasterArray和SpatialArray对象共享。因此两者的一些维度更变方法是通用的。

RasterArray对象的元素数量通过length（）查看：

length(dems)

 行列数的查询方法：

nrow(coasts)
ncol(coasts)

以上三个方法都汇总为dim（）。此方法不同于向量的dim（）。

dim(dems)
dim(coasts)

 

 RasterArray可以很便捷地通过name来编排。names（），colnames（），rownames（）和dimnames（）。

names（）常用于类似向量的RasterArray对象，colnames（）和rownames（）常用于类似矩阵的SpatialArray对象。dimnames（）是前三种方法的汇总。

names(dems)
colnames(coasts)
rownames(coasts)

 以上的name方法也可用于重命名：

coasts2 = coasts
colnames(coasts2) = c("m", "c")
coasts2

 如同操作普通数组一样，也可以用dimnames（）查询或重命名：

dimnames(coasts2)[[1]] = paste(rownames(coasts), "Ma")
coasts2

 RasterArray和SpatialArray的元素除了names外，每个layer都有name，可以通过layers（）：

layers(dems)

RasterLayer栅格总数或堆栈（stack）的栅格总数可以通过ncell（）和nvalues（）：

ncell(dems)
nvalues(dems)

构建子集和重组

 RasterArray和SpatialArray对象的首要目的是便捷地获取条目（items）。或是对stack专注于layer（双方括号），或是RasterArray专注于元素（单方括号）。

1. 挑选图层（layer）——双方括号

输出的是一个RasterLayer或RasterStack对象。挑选的方法有通过下标、通过名称、或通过位置的逻辑值。

可以通过名称选择RasterArray的图层。plot（）或mapplot（）将选择的图层可视化：

one = dems[["dem_45"]]
mapplot(one, col="earth")

 双方括号会返回一个RasterArray。用多个元素创建子集会返回一个RasterStack。

dems[[c(1,2)]]

 2. 单方括号

如果RasterArray是一维数组，单方括号会检索并替换目标

dems["30"]

 返回30Ma图层。默认情况下，RasterArray容器会清除掉，但将drop设为FALSE即可保留。

dem30 = dems["30", drop=FALSE]
class(dem30)

获取一维RasterArray的元素：

dems[4:12]

如果要从多维的RasterArray或SpatialArray中获取一个图层：

one = coasts["15", "coast"]
mapplot(one)

 

 或者获取整行/整列：

coasts["15", ]
coasts[, "coast"]

从raster和sp中继承的功能

 1. RasterArray

RasterArray的RasterLayer可以裁剪为单行数据：

# 设置裁剪区域
ext = extent(c(xmin=106.58, xmax=157.84, ymin=-45.23, ymax=1.14))
# 对dems进行裁剪
au = crop(dems, ext)
mapplot(au[1], col="earth")

 还可以重新调整分辨率：

# 调整分辨率
template = raster(res=5)
coarse = resample(dems, template)
# 可视化
mapplot(coarse["45"], col="earth")

 2. SpatialArray

对于SpatialArray也一样。可以对整个集合的投影进行修改：

# 将coasts的矩形图转为椭圆图
mollCoasts = spTransform(coasts, CRS("+proj=moll"))
# 设置底图
edge = spTransform(mapedge(), CRS("+proj=moll"))
# 绘制底图
plot(edge, col=ocean(8)[5])
# 绘制大陆架
plot(mollCoasts["20", "margin"], col=ocean(8)[7], border=NA, add=T)
# 绘制海岸线
plot(mollCoasts["20", "coast"], col=terra(8)[5], add=T, border=NA)

绘图

mapplot（）

调色板

showPal（）

 

RasterLayer

mapplot（）可作用于RasterLayer，mapplot函数作用于Raster*对象时会包括一个默认调色板并省略图例、坐标轴和边框。

 

 使用特定的调色板：

mapplot(dems[1], col="earth", main="Using the earth palette")

 3. RasterArray

mapplot（）可以将RasterArray对象绘制为多图

data(dems)
mapplot(dems)

添加一个统一的图例

 

RasterArray处理NA图层

如果有NA图层，那么在该图层位置上会显示“Plot not available”。

data(dems)
# 将一图层设置为NA
dems[5] = NA
is.na(dems)

 

列数和页数

默认绘图为3列，也可用ncol参数设置：

 

 有时，在一张图上会绘制太多个地图。可用multi参数设置为TRUE将地图绘制在多个图上：