R实践——【rgplates】安装、介绍、入门

1. rgplates 安装

1.1 easy way

install.packages("rgplates")

1.2 备案方法

• 从GITHUB安装稳定版本

install.packages(
    "https://github.com/adamkocsis/rgplates/raw/main/_archive/source/rgplates_0.2.1.tar.gz", 
    repos=NULL, type="source")

• 从GITHUB安装开发版本
  如果您遇到任何地方都没有描述的软件包技术问题，您可能需要查看开发版本。

如果您希望安装开发版本，我建议您手动安装：

1. 将存储库克隆到您的本地硬盘。
2. 打开一个终端并导航到您克隆的目录。从那里应该可以看到 rgplates 目录。
3. 在终端中运行这一行

R CMD INSTALL rgplates

• 如果您看到提示未找到 R 的错误，则必须将其添加到您的 PATH 环境变量中。
• 如果rgplates依赖的R包没有安装，就得手动安装，否则会报错。

2. rgplates 介绍

该包开箱即用，依赖于 GPlates Web 服务 (GWS)，这是一种使用 URL 中提供的数据执行古地理重建的在线服务。rgplates (1) 将数据发送到 GWS，GWS (2) 计算古坐标，然后 (3) rgplates 读取返回的结果。 要使用这个在线服务，您必须连接到互联网。

3. rgplates 在线方法入门

3.1 加载rgplates

使用前，加载此包

library(rgplates)

加载后，rgplates 会自动加载 Simple Features for R (sf) 包，这是一个用于处理矢量空间数据的标准 R 包。

3.2 板块重建

所有古坐标重建都依靠 reconstruct() 函数完成。 它使用的默认模型是由 C. Scotese 开发的 PALEOMAP 项目的模型。

每个构造模型都依赖于在地球表面旋转的板块。板块位置可以重建到模型涵盖的任何年龄。可以使用reconstruct()函数查询板块的当前位置，将字符串“plates”作为第一个参数传递，并将目标age（年龄）设置为0 (以Ma为单位)。

pl0 = reconstruct("plates", age=0)
pl0

Simple feature collection with 501 features and 0 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 90
Geodetic CRS:  WGS 84
First 10 features:
                         geometry
1  POLYGON ((-155.2068 19.5404...
2  POLYGON ((-77.6715 24.7979,...
3  POLYGON ((-77.6715 24.7979,...
4  POLYGON ((-71.9396 20.9241,...
5  POLYGON ((-78.3856 27.3354,...
6  POLYGON ((-180 90, 180 90, ...
7  POLYGON ((164.2599 89.0123,...
8  POLYGON ((-106.0379 21.2018...
9  POLYGON ((-112.1621 28.3617...
10 POLYGON ((-108.7002 23.7579...

gplates 包依赖 sf 来处理矢量空间数据——比如板块。你可以用它做任何你通常可以用和 sf 对象做的事情：操纵它，改变它的地图投影，导出或使用它进行计算。这些数据是标准的等矩形投影，配有经度和纬度数据。

您可以使用plot()函数绘制板块分布图。要关注空间数据而不是要素的属性(本质上是技术性的)，可以绘制该对象的几何图形。

plot(pl0$geometry)

设置年龄参数允许您访问过去时间点的板块构造状态。年龄接受以百万年为单位的日期。重建三叠纪/侏罗纪边界附近的板块位置(大约200Ma)，您必须将年龄参数设置为200:

pl200 = reconstruct("plates", age=200)
pl200

Simple feature collection with 142 features and 0 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -77.0625 xmax: 180 ymax: 90
Geodetic CRS:  WGS 84
First 10 features:
                         geometry
1  POLYGON ((17.3375 17.8127, ...
2  POLYGON ((-19.2369 72.1231,...
3  POLYGON ((-19.3706 71.9604,...
4  POLYGON ((-19.3706 71.9604,...
5  POLYGON ((-2.4637 32.5774, ...
6  POLYGON ((-19.7326 72.67, -...
7  POLYGON ((-19.7326 72.67, -...
8  POLYGON ((-19.6809 73.003, ...
9  POLYGON ((-18.0301 72.7282,...
10 POLYGON ((-18.6014 71.6729,...

您可以用类似的方式绘制结果：

plot(pl200$geometry)

请注意，与今天 (0Ma) 的结果不同，过去的板块配置不会返回海洋板块。（大部分模型不返回大洋板块，这是PALEOMAP模型的一个特点）

同样，pl0和pl200对象都是sf类对象。你可以自定义他们的绘图，就像你对任何其他sf对象一样——例如为多边形设置一个填充颜色，而不绘制它们的边界。

plot(pl200$geometry, col="gray", border=NA)

3.3 独立的地点坐标

这很好，但单独绘制板块的用处也仅有此了。 reconstruct() 的真正用途是能够计算给定年龄的当前位置的古坐标。

3.3.1 单个现存坐标点

我们来考虑一下伦敦的位置！无论是在网址里，还是用谷歌地图的用户界面，都可以查到市中心的坐标在北纬51.52(纬度)，西经0.38(经度)左右。为了找出这个城市在三叠纪/侏罗纪界线的位置，你需要:(1)在R中输入这些坐标；以及(2)将它们提供给reconstruct()函数。

请注意，我们在这里处理的是全球范围的近似坐标。为了获得更精确的结果，请始终确保点的CRS与地图（包括椭圆）的CRS匹配！

因为在通常使用的等矩形投影中，绘图的x轴为经度，y轴为纬度，所以我们通常先按经度的顺序记录坐标，然后按纬度的顺序记录坐标(东经度和北纬度记录为正值)。

为了使这种结构绝对清晰，最好将坐标输入为两列矩阵，第一列是经度，第二列是纬度：

london = c(-0.38, 51.52)
londonMat = matrix(london, ncol=2, byrow=TRUE)
colnames = c("long", "lat")
londonMat

      [,1]  [,2]
[1,] -0.38 51.52

因为地图的坐标参考系统(CRS)是经度-纬度，所以您可以使用这些坐标直接在当前地图上用点()表示城市的位置x，在本例中是用红色加号：

plot(pl0$geometry)
points(londonMat, col="red", pch=3)

3.3.2 单个点的古坐标

找到这些地点的古坐标就像计算板块的坐标一样简单。您必须使用之前定义的矩阵 (londonMat) 作为 reconstruct() 的第一个参数 (之前给出了“plates”) ，并提供以百万年为单位的目标年龄：

londonMat200 = reconstruct(londonMat, age=200)
londonMat200

     paleolong paleolat
[1,]    2.3267  38.1654

此计算的结果是一个类似的矩阵：包括 paleolongitude 和 paleolatitude 列。如果您将坐标作为普通矩阵提供，则坐标会被推断为经度和纬度。

你可以用同样的方式来可视化这一点，就像你可视化该地点的当前位置一样：

plot(pl200$geometry, col="gray", border=NA)
points(londonMat200, col="red", pch=3)

3.3.3 多个点的古坐标

如果要重建多个位置，只需在矩阵中提供更多的行。例如，如果您还想计算澳大利亚悉尼(南纬33.85度，东经151.11度)和加拿大蒙特利尔(北纬45.52度，西经73.61度)的位置，您需要以类似的方式添加这些位置。

sydney = c(151.17, -33.85)
montreal = c(-73.61, 45.52)
cities = rbind(london, sydney, montreal)
colnames(cities) = c("long", "lat")

现在我们有了经度和纬度的矩阵，我们需要做的就是使用它作为 reconstruct() 函数的第一个参数：

cities200 = reconstruct(cities, age=200)

请注意，实体的顺序保持不变。我们可以用类似于单个城市的方法来绘制这些图。

plot(pl200$geometry, col="gray", border=NA)
points(cities200, col="red", pch=3)

3.4 现今的海岸线

某些模型（包括 PaleoMAP 模型）允许您访问当今海岸线的重建位置，这可能有助于您更好地定位重建板块。如果将字符串“coastlines”作为 reconstruct() 函数的第一个参数传递，则可以获得这些数据：

coast200 = reconstruct("coastlines", age=200)

这将返回一个 sf 类对象，类似于plates。您可以使用 plot() 可视化这些并将它们绘制在其他地图之上，但您必须设置 add=TRUE，否则将创建一个新的绘图。

plot(pl200$geometry, col="gray", border=NA)
points(cities200, col="red", pch=3)
plot(coast200$geometry, add=TRUE)

3.5 其他的重建模型

如果你查看reconstruct（）引用的细节，你会发现除了“PALEOMAP”之外，模型参数还可以设置为其他字符串。这些表示可以通过GPLates Web服务访问的其他模型。例如，如果您想使用Seton等人2012模型在200Ma下执行相同的计算来重建板的位置，您所要做的就是设置model=“SETON2012”：

pl200seton = reconstruct("plates", age=200, model="SETON2012")

您可以将这个结果与PALEOMAP模型进行比较，方法是将这个结果绘制在具有一定透明度的结果之上(例如，在HTML RGBA中为半透明的红色:" #FF000077 ")：

plot(pl200$geometry, col="gray", border=NA)
plot(pl200seton$geometry, col="#FF000077", add=TRUE, border=NA)

你可以看到，这两个重建有很大的不同，随着时间的推移，这一点变得更加明显。此外，通过透明度，你们可以看到板块是如何在收敛区重叠的。

3.6 在线方法的缺点

尽管设置起来很容易，但在线方法有几个局限性。为了更完全地控制旋转，您可能想尝试离线方法。

如果想坚持使用在线方法，您需要了解以下限制因素：

3.6.1 小数位和数据质量

在线方法的数据通过 URL 传递到 GWS，这对重建施加了一些技术限制。第一个限制是年龄必须四舍五入为整数，这可能会影响结果的质量。

第二个潜在的问题是，可以传输的数据量有些有限。这对于数百个点来说很好，但是对于重建数千个坐标来说，您肯定需要使用离线方法。

离线方法可以毫无问题地处理17.24Ma这样的年龄，以及任意数量的数据。

3.6.2 预建模型

虽然在线方法很容易设置，但它使用预定义的模型来实际执行重建。如果您想使用 GWS 尚未实现的自定义或其他模型，则需要使用离线方法。

3.6.3 网络连接

由于计算不是在您的计算机上执行的，而是由托管 GWS 的服务器执行的，因此如果您没有互联网连接，则无法使用此方法。此外，如果您考虑再现性，很容易看出您无法控制该服务在未来是否仍然可用。

有一些docker容器允许你实例化GWS，这可能有助于解决这个问题——但是如果你真的想确保你的代码保持可复制性——你可能想尝试离线方法。

3.6.4 其他选项

您还可以自定义一些其他选项，这些选项不会在GWS中公开供您设置，这可能会影响您的结果。

4. rgplates离线方法入门

虽然在线重建方法可以帮助完成快速计算，但 rgplates 包的真正强大之处在于我们所说的离线重建方法。

4.1 介绍

离线方法不依赖于与GPlates网络服务的连接，而是允许你在你的计算机上本地执行板块构造旋转。为此，除了rgplates包之外，您的计算机上还需要两样东西：

• GPlates桌面应用程序，用于重建计算
• 用于重建的实际构造模型文件。

4.2 GPlates桌面应用程序

GPlates是目前我们用来建造、操作和应用板块构造重建的标准程序。GPlates有一个开发得非常好的图形用户界面(GUI ),您可以从本教程的任何内容中独立使用它。如果你想了解更多关于GPlates的知识，你可以在这里找到很多资料。

为了使用GPlates的离线方法，您不必离开R环境或使用gplates的GUI。你只需要确保你的电脑上安装了应用程序，rgplates会处理剩下的事情。

4.2.1 GPlates安装

除了安装rgplates包之外，您还必须在计算机上安装gplates桌面应用程序。在看到您的操作系统需要哪个文件后，您必须向下滚动到页面底部，选择您希望下载的文件，然后按照说明进行操作。

这已经在所有主要的操作系统(Windows、Mac和GNU/Linux)上进行了测试，并且已经在gplates版本2.2和2.3上可靠地工作。rgplates将能够找到GPlates在Mac和GNU/Linux上的开箱即用安装。只要你把GPlates安装到默认的安装目录，在Windows上也不会有任何问题(不要在安装过程中改变建议的内容——除非你自找麻烦！)

4.3 构造模型

要使用离线方法，您的计算机上必须有一个构造模型。实际上，并没有板块构造模型这样的东西。从技术角度来看，每个模型至少包含两个数据项：

• the Partitioning polygons
• and the Reconstruction tree (aka. rotation file)

分割多边形是矢量空间数据(多边形),代表现今构造板块的分布。GPlates通过旋转在地球表面移动这些多边形，如重建树中所述。这些通常是一起开发的(我们通常称之为模型)，但是它们分别存储在不同的文件中。换句话说，GPlates只是执行旋转/重建的程序，但是您还需要有文件

在线方法依赖于GPlates Web服务服务器上的文件。例如，古地图模型的分割多边形是可见的，如果你只要求“板块”，即今天的年龄(即0Ma)。

partPol = reconstruct("plates", age=0)
plot(partPol$geometry)

为了使用离线方法，我们需要

• 在我们的计算机上同时拥有分区多边形和重建树
• 用R表示它们
• 将它们传递给 reconstruct() 函数

4.4 PALEOMAP模型文件

为了简化第一次试验，我们将使用与 rgplates 包一起分发的 PaleoMAP 项目的模型文件。

这些文件位于 zip 存档中包的安装目录中。我们要做的是：

• 创建一个新的临时目录来存储我们的模型文件
• 在那里解压缩模型文件，然后
• 将它们加载到R中。

注意：有更简单的方法来获取这些文件并以其他方式将它们表示到 R 中（例如使用chronosphere，我们将在教程 4 中进行研究。）但出于教学目的，我建议首先使用以下方式执行此操作如下。

4.4.1 创建临时文件夹

临时目录是存储临时数据、文件的最佳位置，它们是在计算过程中创建的，在 R 会话之外不需要（当然你不需要在这里存储这些，你可以使用计算机上的任何目录） .

临时目录很容易创建，你只需要使用 tempdir() 函数，它会为你创建这样一个目录，并告诉你在哪里可以找到它：

td = tempdir()

[1] "C:\\Users\\bailo\\AppData\\Local\\Temp\\RtmpgRghHL"

每个操作系统都会以不同的方式处理这些问题，所以如果您的操作系统与上面的结果不匹配，请不要感到惊讶——实际上，当您运行上面的函数时，每个会话中的情况都可能不同。

4.4.2 解压模型文件

模型的文件很大，比 CRAN 指南推荐的大，所以必须压缩。我们将不得不使用 unzip() 函数解压缩它们——但首先我们必须找到它们的位置。

可以使用 system.file() 函数找到 R 包的安装路径。

rgPath = system.file(package="rgplates")

[1] "D:/ALL_Softwares/R-4.2.0/library/rgplates"

您可以通过使用 list.files() 列出该目录的内容来确认这是正确的路径

list.files(rgPath)

 [1] "CITATION"     "DESCRIPTION"  "extdata"      "help"         "html"         "INDEX"       
 [7] "MD5"          "Meta"         "NAMESPACE"    "NEWS.md"      "R"            "rgplates.pdf"

模型数据位于 extdata 目录中，在存档 paleomap_v3.zip 中。我们必须找到这个存档并将其解压缩到临时目录（exdir 参数）

unzip(file.path(rgPath,"extdata/paleomap_v3.zip"), exdir=td)

您可以通过查看临时目录来确认提取确实发生了：

list.files(file.path(td))

 [1] "20e8cc60218b4655a3881a246b5b14f2.png"            
 [2] "3d5d7aba18ec4ea09ad3afe6d6c3ab4f.png"            
 [3] "617dc5e784fb413bb4f181837cd7bf7b.png"            
 [4] "9b9c99ee368f4989bfb89c9755d96e87.png"            
 [5] "c1e2330ed82e491fbb8afd581c4af79e.png"            
 [6] "d44dcec927d64eb99923a4da09e245cd.png"            
 [7] "ded46f127eb44d17a1b85ecdbbcd66f5.png"            
 [8] "ef2d31134c49404684bf9638119de99c.png"            
 [9] "ffad242897fd47839324f4fa08c96d63.png"            
[10] "PALEOMAP_PlateModel.rot"                         
[11] "PALEOMAP_PlatePolygons.gpml"                     
[12] "PALEOMAP_PoliticalBoundaries.gpml"               
[13] "rs-graphics-9d334d49-dded-4ee1-be86-5b71e69d6008"

有两个文件对我们很重要：一个以 .gpml 结尾（GPlates 标记语言格式的分区多边形）和一个以 .rot 结尾（重建树/旋转文件）。我们知道这些文件的绝对路径，可以在以下位置找到它们：

pathToPolygons = file.path(td, "PALEOMAP_PlatePolygons.gpml")
pathToRotations = file.path(td, "PALEOMAP_PlateModel.rot")

> pathToPolygons
[1] "C:\\Users\\bailo\\AppData\\Local\\Temp\\RtmpgRghHL/PALEOMAP_PlatePolygons.gpml"
> pathToRotations
[1] "C:\\Users\\bailo\\AppData\\Local\\Temp\\RtmpgRghHL/PALEOMAP_PlateModel.rot"

4.5 板块模型类

现在我们有了文件并且知道如何找到它们，我们已经以某种方式在 R 中表示它们，因此我们可以表明我们希望将它们与 GPlates 一起使用来进行板块构造重建。

我们需要创建一个 platemodel 类对象来执行此操作。我们需要提供的是这些文件的路径，分别用于分区多边形和旋转文件：

pm = platemodel(polygons = pathToPolygons,
                rotation = pathToRotations)

GPlates plate tectonic model.
static polygons:  "PALEOMAP_PlatePolygons.gpml" 
rotation:         "PALEOMAP_PlateModel.rot" 

这是一个platemodel类对象，它实际上只是一个包含文件实际位置的包装器。
现在是时候提一下，这就是离线方法真正的灵活性所在。您可以在此处使用任何第三方模型文件，只要您提供文件的适当绝对路径，下面的方法应该有效。请参阅下面的详细信息

4.6 重建板块

重建本身的工作方式与在线方法相同，唯一的区别是我们将使用 platemodel 类对象 pm，而不是使用字符串作为模型参数（默认为“PALEOMAP”）。

以下是使用 PaleoMAP 模型为 200Ma 重建的板块 - 现在使用离线方法计算：

pl0ff200 = reconstruct("plates", age=200, model=pm, path.gplates="D:/ALL_Softwares/GPlates/gplates-2.2.0.exe")

请注意，这是一个类似于我们之前在在线方法中看到的 sf 类对象，尽管它具有更多的属性。不过，您可以使用 plot() 以相同的方式将其可视化。

plot(pl0ff200$geometry, border=NA, col="gray")

4.7 重建地点

重建点的位置与在线方法类似。同样，唯一的区别是您必须提供 platemodel 对象作为模型参数。这是之前计算和可视化一些城市过去位置的示例，现在使用离线方法：

首先，我们获取地点数据：

london = c(-0.38, 51.52)
sydney = c(151.17, -33.85)
montreal = c(-73.61, 45.52)

cities = rbind(london, sydney, montreal)
colnames(cities) = c("long", "lat")

           long    lat
london    -0.38  51.52
sydney   151.17 -33.85
montreal -73.61  45.52

其次，我们用我们的模型进行重建

cities200 = reconstruct(cities, age=200, model=pm, path.gplates="D:/ALL_Softwares/GPlates/gplates-2.2.0.exe")

               long       lat
london     2.326686  38.16536
sydney    79.175489 -60.67032
montreal -20.890271  26.08778

然后我们绘图：

plot(pl0ff200$geometry, col="gray", border=NA)
points(cities200, col="red", pch=3, add=TRUE)

4.8 获取模型文件的其他途径

模型文件可从EarthByte资源页面获得。例如，我们在这里使用的古地图模型的文件可以在https://www.earthbyte.org/paleomap-paleoatlas-for-gplates/.获得。

在页面底部，您会找到一个链接：Link to rasters, reconstruction files and tutorial （指向栅格、重建文件和教程的链接），您可以通过该链接下载 Scotese_Paleoatlas_v3.zip 文件。在此 zip 文件中，在 PALEOMAP Global Plate Model 目录中，您可以看到我们之前拥有的相同文件。

只要提供文件的绝对路径，就可以使用这些文件或任何类似的文件来创建platemodel对象。

4.9 chronosphere

正如你所看到的，处理这些文件是相当乏味的，而我们，rgplates的开发者，很少以这种方式使用它们。在实际应用中，我们推荐使用离线的rgplates方法和chronosphere，你可以在下一个教程中了解更多。

5. 使用chronosphere获取的数据进行重建

最初，rgplates中的所有功能都是在chronosphere项目中开发的。然而，自项目开始以来，这些模块的功能发生了分化:rgplates包含用R执行板块构造重建的功能，而chronosphere是一个数据分发和版本化系统。这两者与构造模型数据相关联:我们将下载这些数据并将其导入chronosphere，并与rgplates一起使用。

5.1 加载包

library(rgplates)
library(chronosphere)

为了表明我们调用了哪些函数，我们将使用::运算符。这个操作符允许我们指出我们想从哪个包中使用函数。

5.2 下载数据

稍后将提供有关如何使用 chronosphere-portal 的详细教程。现在，我们只需要知道如何访问 platemodel 类对象。

从 chronosphere 下载数据很容易：您只需要一行代码，以及您选择的数据的适当 ID。要下载我们之前使用的构造模型（板块多边形和重建树）（PaleoMAP 模型），我们只需要这个命令：

pm = chronosphere::fetch(dat="paleomap", var="model", ver="v3-GPlates")

GPlates plate tectonic model.
PALEOMAP plate tectonic reconstruction files (C. Scotese) 
 - v3 (GPlates resource) 
static polygons:  "PALEOMAP_PlatePolygons.gpml" 
rotation:         "PALEOMAP_PlateModel.rot" 

chronosphere::fetch() 的参数是您要下载的项目的数据集 (dat)、变量 (var) 和版本 (ver) ID。

5.3 实际重建

我们可以立即将其与离线重建方法一起使用。

pl400 = rgplates::reconstruct("plates", age=400, model=pm,
                              path.gplates="D:/ALL_Softwares/GPlates/gplates-2.2.0.exe")

plot(pl400$geometry, col="gray", border=NA)

5.4 结束语

通过三行代码，我们下载了构造信息，计算了旋转并绘制了重建的几何图形。

这就是我们最初打算如何使用reconstruct()函数的离线方法，这也是我们建议作为一般工作流采用的方法。

请注意，我们打算将尽可能多的模型添加到chronosphere中，chronosphere正在进行彻底检查，以使其尽可能灵活和有用。然而，保留所有当前可用的数据并保持所有功能并不容易——但预计将在2023年上半年完成。