(11)脚本、算法和函数

本文由SCY翻译自《Geocomputation with R》第十一章

前提条件

本章前提条件最小，因为它主要使用基础的R语言。sf包用于检查我们将开发的算法的结果，该算法用于计算多边形的面积。就先前的知识而言，本章假设您已经理解了R中的地理数据中介绍的地理类，以及如何使用它们来表示各种各样的输入文件格式(参见地理数据I/O)。

引言

引言章阐明了地理计算不仅仅是使用现有的工具，还包括以“可共享的R脚本和函数”的形式开发新工具。本章教授可复现代码的这些基础构件。它还介绍了低级几何算法，这种算法在GIS软件桥梁章中有应用。阅读本章应有助于您理解这类算法是如何工作的，并编写可多次、由多人、在多个数据集上使用的代码。然而，单凭本章自身无法使您成为一名熟练的程序员。编程是困难的，需要大量的实践。

要理解编程作为一种自成体系的智力活动，您必须转向计算机编程；您必须阅读和编写计算机程序——很多程序。

有充分的理由去学习编程。虽然本章本身并未教授编程——参见多种资源，它们教授R语言和其他语言的编程——但它确实提供了一些着眼于几何数据的起点，这些起点可能是发展编程技能的良好基础。

本章还演示并强调了可重复性的重要性。可重复性的优点不仅仅是允许其他人复制您的工作：相对于只运行一次的代码，可重复的代码在各方面通常都更优秀，包括计算效率、‘可扩展性’（代码在大数据集上运行的能力）以及更容易进行适应和维护。

脚本是可重复R代码的基础，这一主题在脚本节中有覆盖。算法是一系列用于修改输入并得出输出的步骤的配方，如几何算法节所述。为了简化共享和可重复性，算法可以被放置在函数中。这是函数节的主题。找出多边形的质心的例子将用于整合这些概念。几何操作章已经介绍了一个质心函数st_centroid()，但这个例子突出显示了看似简单的操作其实是相对复杂代码的结果，这证实了以下的观察：

关于空间数据问题最吸引人的一点是，对人来说看似极其简单的事情，在计算机上可能出奇地困难。

这个例子也反映了本章的次要目标，即不是“复制现有的内容，而是展示现有内容是如何工作的”。

脚本

如果说包中分发的函数是R代码的基础构件，那么脚本就是将它们粘合在一起的胶水。脚本应该按照逻辑顺序存储和执行，以创建可复现的工作流，这可以手动完成，也可以使用诸如targets这样的工作流自动化工具来完成。如果你是编程新手，当你第一次遇到脚本时，它们可能看起来让人望而却步，但实际上它们只是普通的文本文件。脚本通常保存为一个扩展名代表它们包含的语言的文件，例如，用Python编写的脚本使用.py扩展名，用Rust编写的脚本使用.rs扩展名。R脚本应该保存为.R扩展名，并且名字应该反映它们的功能。一个例子是11-hello.R，这是一个存储在书籍仓库的code文件夹中的脚本文件。11-hello.R是一个非常简单的脚本，只包含两行代码，其中一行是注释。

这些脚本通常是科研过程中非常有用的组成部分，特别是当你需要对大量栅格和矢量数据进行处理和后续分析时，它们可以帮助你实现高度自动化和可复现的工作流程。而使用例如targets这样的工作流工具，可以进一步提高代码的结构性和可维护性。

# Aim: provide a minimal R script
print("Hello geocompr")

这个脚本的内容并不特别令人兴奋，但它说明了一点：脚本不需要复杂。保存的脚本可以通过R命令行中的source()函数完整地调用和执行，如下面所示。这个命令的输出显示，注释会被忽略，但print()命令会被执行：

source("code/11-hello.R")
#> [1] "Hello geocompr"

你也可以从系统命令行shell，例如bash和PowerShell，调用R脚本，前提是RScript可执行文件已经配置好，例如如下：

Rscript code/11-hello.R

关于脚本文件中可以和不能包含什么，并没有严格的规定，也没有什么能阻止你保存错误的、不可复制的代码。不包含有效R代码的代码行应该被注释掉，通过在行的开始处添加#，以防止错误，如11-hello.R脚本中的第一行所示。然而，有一些值得遵循的惯例：

• 按顺序编写脚本：就像电影的剧本一样，脚本应该有一个明确的顺序，比如’设置’、‘数据处理’和’保存结果’（大致相当于电影中的’开头’、‘中间’和’结尾’）。
• 向脚本添加注释，以便其他人（以及你未来的自己）能够理解它。
  至少，注释应该说明脚本的目的（见图@ref(fig:codecheck)）并（对于长脚本）将其划分为几个部分。
  这可以在RStudio\index{RStudio}中通过快捷键Ctrl+Shift+R来完成，它会创建’可折叠’的代码段标题。
• 最重要的是，脚本应该是可复制的：能在任何计算机上运行的自包含脚本比只能在你的计算机上、在好天气下运行的脚本更有用。
  这涉及到在开始时附加所需的包，从持久性来源（如可靠网站）读取数据，并确保已经采取了之前的步骤。^[
  之前的步骤可以用注释或者if语句来指出，比如if (!exists("x")) source("x.R")（如果对象x不存在，则会运行脚本文件x.R）。
  ]

在R脚本中强制可复制性是困难的，但有一些工具可以帮助。默认情况下，RStudio \index{RStudio} 会’代码检查’ R脚本，并用红色波浪线标出有问题的代码，如下图所示：

Code checking in RStudio. This example, from the script 11-centroid-alg.R, highlights an unclosed curly bracket on line 19.

📌A useful tool for reproducibility is the reprex package.
Its main function reprex() tests lines of R code to check if they are reproducible, and provides markdown output to facilitate communication on sites such as GitHub.
See the web page reprex.tidyverse.org for details.

本节的内容适用于任何类型的R脚本。对于地理计算脚本的特殊考虑是，它们往往具有外部依赖性，例如在第@ref(read-write)章中重度使用的用于处理地理数据的核心R包所需的GDAL依赖，用于数据导入和导出。运行更专业的地理算法可能需要GIS软件依赖，如第十章所概述。处理地理数据的脚本还常常要求输入数据集以特定格式提供。这种依赖应在脚本的注释中或其所属项目的其他地方进行说明，如脚本 11-centroid-alg.R 所示。‘防御性’ 编程技巧和良好的错误信息可以通过检查依赖性和与用户沟通（如果某些需求未得到满足）来节省时间。如果语句，用R中的if ()实现，可用于发送消息或运行代码行，只有在满足某些条件时才会这样做。例如，以下代码行会在缺少某个文件时向用户发送消息：

if (!file.exists("required_geo_data.gpkg")) {
  message("No file, required_geo_data.gpkg is missing!")
} 
#> No file, required_geo_data.gpkg is missing!

11-centroid-alg.R脚本所进行的工作在下面的可复制示例中有所展示，该示例创建了一个名为poly_mat的预备对象，代表一个边长为9单位的正方形。这个例子显示了source()可以与URLs一同工作，假设你有互联网连接。如果没有，相同的脚本可以通过source("code/11-centroid-alg.R")来调用，前提是你之前已经下载了github.com/geocompx/geocompr仓库，并且你正在geocompr文件夹中运行R。

poly_mat = cbind(
  x = c(0, 9, 9, 0, 0),
  y = c(0, 0, 9, 9, 0)
)
# Short URL to code/11-centroid-alg.R in the geocompr repo
source("https://t.ly/0nzj")

#> [1] "The area is: 81"
#> [1] "The coordinates of the centroid are: 4.5, 4.5"

几何算法

我们可以这么理解算法，它相当于烘焙食谱。它们是一整套指导方针，当对输入进行操作时，会得到有用/美味的结果。输入在烘焙的情况下是如面粉和糖的成分，在算法的情况下是数据和输入参数。而烘焙食谱可能导致美味的蛋糕，成功的算法应有带来环境/社会/其他利益的计算结果。在深入可复制示例之前，下面简短的历史将展示算法是如何与脚本（在脚本节中涉及）和函数（可以用来概括算法，使其更便携和易于使用，我们将在函数节中看到）相关联的。

"算法"这个词源自9世纪在巴格达出版的早期数学教科书 Hisab al-jabr w’al-muqabala。这本书被译成拉丁文，并变得非常流行，以至于作者的姓氏，al-Khwārizmī，“被永久地当作一个科学术语来记忆：Al-Khwarizmi 变成了 Alchoarismi、Algorismi，最终变成了 algorithm” [@bellos_alex_2011]。在计算时代，算法指的是解决问题的一系列步骤，从而得出预定义的输出。输入必须在适当的数据结构中正式定义[@wise_gis_2001]。算法通常从流程图或伪代码开始，展示过程的目标，然后在代码中实现。为了简化可用性，常见的算法通常被封装在函数内，这些函数可能会隐藏一些或所有已采取的步骤（除非你查看函数的源代码，参见函数节）。

地理算法，比如我们在GIS章节中遇到的，是接收地理数据并通常返回地理结果的算法（同样的东西还有其他术语，比如GIS算法和几何算法）。这可能听起来简单，但这是一个深刻的主题，有一个专门的学术领域，计算几何，致力于它们的研究[@berg_computational_2008]，以及关于该主题的大量书籍。@orourke_computational_1998，例如，使用可复制和免费提供的C代码，以一系列逐渐复杂的几何算法介绍该主题。

几何算法的一个例子是找出多边形的质心的算法。有许多方法来计算质心，其中一些仅适用于特定类型的空间数据。为了本节的目的，我们选择一种容易可视化的方法：将多边形分解成许多三角形，并找出每个三角形的质心，这种方法由 @kaiser_algorithms_1993 讨论，还简短地在 @orourke_computational_1998 中提到。在编写任何代码之前，进一步将这种方法分解成离散任务会有所帮助（后来称为步骤1到步骤4，这些也可以以示意图或伪代码的形式呈现）：

1. 将多边形分成相邻的三角形。
2. 找出每个三角形的质心。
3. 找出每个三角形的面积。
4. 找出三角形质心的面积加权平均值。

这些步骤可能听起来很简单，但将单词转换成工作代码需要一些工作和大量的反复试验，即使输入有限制：算法仅适用于凸多边形，它们不包含大于180°的内角，不允许星形（decido和sfdct软件包可以使用外部库对非凸多边形进行三角剖分，如algorithm 小册子中所示，该小册子托管在 geocompx.org 上）。

表示多边形最简单的数据结构是一个矩阵，其中x和y坐标在每行表示一个顶点，以追踪多边形边界的顺序，其中第一行和最后一行是相同的[@wise_gis_2001]。在这种情况下，我们将用基础的R创建一个有五个顶点的多边形，该多边形建立在GIS Algorithms [@xiao_gis_2016参见github.com/gisalgs 以获取 Python 代码]的一个示例上，如下图所示。

# generate a simple matrix representation of a polygon:
x_coords = c(10, 20, 12, 0, 0, 10)
y_coords = c(0, 15, 20, 10, 0, 0)
poly_mat = cbind(x_coords, y_coords)

现在我们有了一个示例数据集，我们准备进行上面概述的第1步。下面的代码展示了通过创建一个单一的三角形（T1）来实现这一点，该代码演示了该方法；它还通过基于公式$1/3(a + b + c)$ 来计算其质心，从而演示了第2步，其中$a$到$c$是表示三角形顶点的坐标：

# create a point representing the origin:
Origin = poly_mat[1, ]
# create 'triangle matrix':
T1 = rbind(Origin, poly_mat[2:3, ], Origin) 
C1 = (T1[1,] + T1[2,] + T1[3,]) / 3

Illustration of polygon centroid calculation problem.

第3步是找出每个三角形的面积，以便计算一个加权平均值，该值会考虑到大三角形的不成比例影响。计算三角形面积的公式如下[@kaiser_algorithms_1993]：

$$
\frac{A_x ( B_y − C_y ) + B_x ( C_y − A_y ) + C_x ( A_y − B_y )}{ 2 }
$$

式中，$A$到$C$是三角形的三个点，而$x$和$y$分别指代 x 和 y 的维度。将这个公式翻译成能够处理矩阵表示形式的三角形 T1 数据的 R 代码如下（函数 abs() 确保了结果为正）：

# calculate the area of the triangle represented by matrix T1:
abs(T1[1, 1] * (T1[2, 2] - T1[3, 2]) +
    T1[2, 1] * (T1[3, 2] - T1[1, 2]) +
    T1[3, 1] * (T1[1, 2] - T1[2, 2])) / 2
#> [1] 85

该代码块输出了正确的结果。^[可以用以下公式（该公式假设底边是水平的）来验证结果：面积是底边宽度与高度乘积的一半，$A = B * H / 2$。在这种情况下 $10 * 10 / 2 = 50$。]问题在于代码比较笨拙，如果我们想要在另一个三角形矩阵上运行它，就必须重新键入。为了使代码更具通用性，我们将在函数节中看到如何将其转换为一个函数。

第4步要求在所有三角形上（如上面所示），而不仅仅是一个三角形上，进行第2步和第3步。这就需要迭代以创建表示多边形的所有三角形，如下图所示。这里使用 lapply()和 vapply()来迭代每个三角形，因为它们在基础 R 中提供了一个简洁的解决方案：^[有关文档，请参见 ?lapply，更多关于迭代的信息，请参见第@ref(location)章。]

i = 2:(nrow(poly_mat) - 2)
T_all = lapply(i, function(x) {
  rbind(Origin, poly_mat[x:(x + 1), ], Origin)
})

C_list = lapply(T_all,  function(x) (x[1, ] + x[2, ] + x[3, ]) / 3)
C = do.call(rbind, C_list)

A = vapply(T_all, function(x) {
  abs(x[1, 1] * (x[2, 2] - x[3, 2]) +
        x[2, 1] * (x[3, 2] - x[1, 2]) +
        x[3, 1] * (x[1, 2] - x[2, 2]) ) / 2
  }, FUN.VALUE = double(1))

Illustration of iterative centroid algorithm with triangles. The X represents the area-weighted centroid in iterations 2 and 3.

现在，我们有条件完成第4步，使用sum(A)来计算总面积，以及使用 weighted.mean(C[, 1], A) 和 weighted.mean(C[, 2], A) 来计算多边形的质心坐标（给警觉的读者一个练习：验证这些命令是否有效）。为了展示算法与脚本之间的联系，本节的内容已经被压缩成 11-centroid-alg.R。我们在脚本节末尾看到，这个脚本如何计算一个正方形的质心。编写脚本 的优点是，它适用于新的 poly_mat 对象（参见下面的练习，以 st_centroid() 为参考来验证这些结果）：

source("code/11-centroid-alg.R")
#> [1] "The area is: 245"
#> [1] "The coordinates of the centroid are: 8.83, 9.22"

上面的示例表明，使用基础 R 从基础原理出发，确实可以开发出低级地理操作。它还表明，如果已经存在一个经过验证的解决方案，那么重新发明轮子可能是不值得的：如果我们的目标仅仅是找到多边形的质心\index{centroid}，那么将poly_mat表示为一个sf对象并使用现有的sf::st_centroid() 函数可能会更快。然而，从1^st^ 原理编写算法的巨大好处是，你将理解整个过程的每一个步骤，这是使用其他人代码时无法保证的。另一个需要考虑的因素是性能，与低级语言如C++相比，R在进行数字计算方面可能较慢（参见地理计算软件节），并且难以优化。如果目的是开发新方法，那么计算效率不应该是优先考虑的。这一点体现在“过早优化是编程中一切（或至少大部分）邪恶的根源”[@knuth_computer_1974]这一说法中。

算法开发是困难的。这一点应该从开发一个仅仅是一种相对低效的解决方案、对实际问题（实际中凸多边形并不常见）具有有限应用的基础 R 的质心算法中显而易见。这种经验应该会让人更加重视像GEO（它是 sf::st_centroid() 的基础）和CGAL（计算几何算法库）这样的低级地理库，这些库不仅运行快，还适用于多种类型的输入几何。
这些库开源的一个重要优点是，其源代码便于学习、理解，并且（对于具备技能和信心的人）进行修改。^[事实上，CGAL函数 CGAL::centroid() 是由7个子函数组成的，如https://doc.cgal.org/latest/Kernel_23/group__centroid__grp.html 所描述，使其能够适用于多种类型的输入数据，而我们创建的解决方案仅适用于一种非常特定的输入数据类型。GEOS 函数 Centroid::getCentroid() 的底层源代码可以在 https://github.com/libgeos/geos/search?q=getCentroid 上找到。]

函数

与算法类似，函数接收一个输入并返回一个输出。然而，函数是特定编程语言中实现的，而不是“配方”本身。在 R 中，函数本身就是对象，可以以模块化的方式创建和组合。例如，我们可以创建一个执行我们质心生成算法第二步的函数，如下所示：

t_centroid = function(x) {
  (x[1, ] + x[2, ] + x[3, ]) / 3
}

上面的示例演示了函数的两个关键组件：1）函数 主体，即定义函数如何处理输入的大括号内的代码；和 2） 形式参数，即函数使用的参数列表——在这种情况下是 x（第三个关键组件，环境，超出了本节的范围）。默认情况下，函数返回最后一个被计算的对象（在 t_centroid() 的情况下是质心的坐标）。^[你也可以通过在函数主体中添加 return(output) 明确设置函数的输出，其中 output 是要返回的结果。]

该函数现在适用于你传递给它的任何输入，如下面的命令所示，该命令计算了上一节中示例多边形中第1个三角形的面积：

t_centroid(T1)
#> x_coords y_coords 
#>     14.0     11.7

我们还可以创建一个函数\index{function}来计算三角形的面积，我们将其命名为 t_area()：

t_area = function(x) {
  abs(
    x[1, 1] * (x[2, 2] - x[3, 2]) +
    x[2, 1] * (x[3, 2] - x[1, 2]) +
    x[3, 1] * (x[1, 2] - x[2, 2])
  ) / 2
}

注意，在创建函数之后，可以用一行代码计算三角形的面积，避免了冗长代码的重复：函数是一种 泛化 代码的机制。新创建的函数t_area()接受任何对象x，假设其具有与我们一直在使用的“三角形矩阵”数据结构相同的维度，并返回其面积，如下面在T1上的示例：

t_area(T1)
#> [1] 85

我们可以通过使用它来找到一个新的三角形矩阵的面积来测试函数的泛化性，该三角形矩阵的高度为1，底边为3：

t_new = cbind(x = c(0, 3, 3, 0),
              y = c(0, 0, 1, 0))
t_area(t_new)
#>   x 
#> 1.5

函数的一个有用特性是它们是模块化的。只要你知道输出会是什么，一个函数就可以用作另一个函数的构建块。因此，t_centroid() 和 t_area()可以用作更大的函数的子组件，以执行脚本11-centroid-alg.R的工作：计算任何凸多边形的面积。下面的代码块创建了 poly_centroid() 函数，以模仿针对凸多边形的 sf::st_centroid() 的行为：^[注意，我们创建的函数在 lapply()和 vapply() 函数调用中被迭代地调用。]

poly_centroid = function(poly_mat) {
  Origin = poly_mat[1, ] # create a point representing the origin
  i = 2:(nrow(poly_mat) - 2)
  T_all = lapply(i, function(x) {rbind(Origin, poly_mat[x:(x + 1), ], Origin)})
  C_list = lapply(T_all, t_centroid)
  C = do.call(rbind, C_list)
  A = vapply(T_all, t_area, FUN.VALUE = double(1))
  c(weighted.mean(C[, 1], A), weighted.mean(C[, 2], A))
}

poly_centroid(poly_mat)
#> [1] 8.83 9.22

像 poly_centroid() 这样的函数可以进一步扩展以提供不同类型的输出。例如，要将结果作为类 sfg 的对象返回，可以使用 ‘包装器’ 函数来修改 poly_centroid() 的输出，然后返回结果：

poly_centroid_sfg = function(x) {
  centroid_coords = poly_centroid(x)
  sf::st_point(centroid_coords)
}

我们可以通过以下方式验证输出与 sf::st_centroid() 的输出相同：

poly_sfc = sf::st_polygon(list(poly_mat))
identical(poly_centroid_sfg(poly_mat), sf::st_centroid(poly_sfc))
#> [1] FALSE

编程

在本章中，我们的速度很快，从脚本到函数，再到算法这个棘手的主题。我们不仅在抽象层面上讨论了它们，还创建了针对特定问题的工作示例：

• 介绍并演示了在’多边形矩阵’上运行的脚本11-centroid-alg.R
• 描述了允许此脚本工作的各个步骤，被称为算法\，一个计算性的配方
• 为了通用化这个算法，我们将其转换成模块化的函数，最终组合成上一节中的函数poly_centroid()

每一个可能看似简单。然而，熟练的编程是复杂的，涉及将每个元素——脚本、算法和函数——组合成一个系统，具有效率和风格。最终的结果应该是健壮且用户友好的工具，供其他人使用。如我们预期本书的大多数读者将是编程新手，能够遵循并复现前面几节的结果是一个重大成就。编程需要许多小时的专门学习和实践才能熟练。

面对希望以有效方式实现新算法的开发者的挑战，可以通过考虑创建一个简单但并非用于生产的函数所付出的努力量来考虑：在当前状态下，poly_centroid()在大多数（非凸）多边形上失败！这提出了一个问题：如何通用化这个函数？两个选项是（1）找到三角化非凸多边形的方法（这是本章支持的在线算法文章所涵盖的主题）和（2）探索其他不依赖于三角网格的质心算法。

一个更广泛的问题是：当已经有高性能的算法被实现并打包成诸如st_centroid()这样的函数时，是否值得编程解决方案？在这个特定情况下，简单的答案是’否’。在更广泛的背景下，考虑到学习编程的好处，答案是’视情况而定’。在编程中，很容易浪费时间尝试实现一种方法，只有到最后才发现有人已经做了艰苦的工作。你可以把这一章看作是通往几何算法编程魔法的垫脚石。然而，它也可以被看作是何时尝试编程一个通用解决方案，以及何时使用现有的更高级解决方案的一堂课。肯定会有编写新函数是最佳方案的时候，但也会有使用已经存在的函数是最佳方案的时候。

“不要重复发明轮子”一样适用于编程，甚至更适用于生活中的其他方面。项目开始时进行一点研究和思考可以帮助决定编程时间最好如何使用。以下三个原则也可以帮助最大化编写代码时的努力，无论它是一个简单的脚本还是由数百个函数组成的包：

1. DRY（不要重复自己）：尽量减少代码重复，并以更少的代码行解决特定问题。
   这一原则在《R for Data Science》[@grolemund_r_2016]的函数章节中有解释。
2. KISS（保持简单愚蠢）：这一原则建议首先尝试简单解决方案，优先于复杂解决方案，需要时使用依赖项，目标是保持脚本简洁。
   这一原则是名言“事物应尽量简单，但不能更简单”的计算类比。
3. 模块性：如果你的代码被划分成明确定义的片段，它将更容易维护。一个函数应该只做一件事，但要做得非常好。如果你的函数变得太长，考虑将其拆分成多个小函数，每个都可以用于其他目的，支持DRY和KISS原则。

我们不能保证这一章会立即让你能够为你的工作创建完美的函数。然而，我们确信其内容将帮助你决定何时是适当的尝试时机（当没有其他现有的函数解决问题，当编程任务在你的能力范围内，当解决方案的好处可能大于开发它的时间成本）。通过使用上述原则，结合通过完成上面的实例获得的实践经验，你将建立你的脚本编写、包编写和编程技能。编程的第一步可能会很慢（下面的练习不应该匆忙完成），但长期的回报可能会很大。

练习

E1. 阅读脚本 11-centroid-alg.R，它位于书籍GitHub仓库的 code 文件夹中。

• 它遵循了在章节1.3中介绍的哪些最佳实践？
• 在你的电脑上通过一个IDE\index{IDE}（比如 RStudio\index{RStudio}）创建这个脚本的一个版本（最好是通过逐行键入脚本，用你自己的编码风格和你自己的注释，而不是复制粘贴——这将帮助你学会如何键入脚本）。以一个方形多边形为例（例如，通过 poly_mat = cbind(x = c(0, 9, 9, 0, 0), y = c(0, 0, 9, 9, 0)) 创建），逐行执行脚本。
• 可以对脚本进行哪些更改以使其更具可重复性？
• 如何改进文档？

E2. 在几何算法部分，我们计算了由 poly_mat 表示的多边形的面积和地理质心分别为245和8.8, 9.2。

• 参考脚本 11-centroid-alg.R 在你自己的电脑上重现结果（奖励：键入命令 - 尽量避免复制粘贴）。
• 结果正确吗？通过将 poly_mat 转换为名为 poly_sfc 的 sfc 对象（使用 st_polygon()（提示：此函数接受 list() 类对象）），然后使用 st_area() 和 st_centroid() 来验证结果。

E3. 之前提到我们创建的算法只适用于凸包。定义凸包（见几何运算章节），并测试该算法在一个非凸包的多边形上的效果。

• 奖励1：思考为什么该方法仅适用于凸包，并注意要对算法做哪些更改才能使其适用于其他类型的多边形。
• 奖励2：基于 11-centroid-alg.R 的内容，仅使用基础R函数编写一个算法，找出以矩阵形式表示的线串的总长度。

E4. 在函数部分，我们创建了不同版本的 poly_centroid() 函数，它们生成了 sfg 类的输出（poly_centroid_sfg()）和类型稳定的 matrix 输出（poly_centroid_type_stable()）。
进一步扩展该函数，创建一个版本（例如，名为 poly_centroid_sf()），该版本是类型稳定的（仅接受 sf 类的输入）并且返回 sf 对象（提示：你可能需要用命令 sf::st_coordinates(x) 将对象 x 转换为矩阵）。

• 通过运行 poly_centroid_sf(sf::st_sf(sf::st_sfc(poly_sfc))) 来验证它的工作原理
• 当你尝试运行 poly_centroid_sf(poly_mat) 时，你得到什么错误信息？