(7)地理数据重投影

发表于 更新于 分类于 ， 阅读次数： Waline： 本文字数： 11k 阅读时长 ≈ 40 分钟

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

本章演示了如何设置并将地理数据从一个CRS转换到另一个CRS，并进一步突出了由于忽略CRS而可能出现的特定问题，当您的数据以经纬度坐标存储时，应该特别注意。

前提条件

  • 本章要求下列包：
1
2
3
4
5
library(sf)
library(terra)
library(dplyr)
library(spData)
library(spDataLarge)

引言

前面章节介绍了坐标参照系统（CRSs），主要集中在两个主要类型：地理（‘lon/lat’，其单位为经度和纬度的度数）和投影（通常以米为单位，从基准点开始）坐标系统。本章在此基础上进一步拓展。它演示了如何设置并将地理数据从一个CRS转换到另一个CRS，并进一步突出了由于忽略CRS而可能出现的特定问题，您应该特别注意，特别是当您的数据以经纬度坐标存储时。

在许多项目中，不需要担心不同的CRS，更不用说在不同的CRS之间转换了。了解您的数据是在投影坐标系统还是地理坐标系统中，以及这对几何操作的影响非常重要。然而，如果您知道您的数据的CRS以及几何操作的影响（在下一节中介绍），CRS应该在幕后正常工作，人们通常在事情出错时突然需要了解CRS。拥有明确定义的CRS，并且所有项目数据都在其中，再加上理解如何以及为什么使用不同的CRS，可以确保事情不会出错。此外，学习坐标系统将加深您对地理数据集及其有效使用方法的了解。

本章教授CRS的基础知识，展示了使用不同CRS的后果（包括可能出错的情况），以及如何将数据集从一个坐标系统“重新投影”到另一个坐标系统。在下一节中，我们将介绍R中的CRS，展示如何获取和设置与空间对象关联的CRS。通过创建缓冲区的实际示例演示了了解数据所在CRS的重要性。本章解决了何时重新投影以及使用哪个CRS的问题。本章介绍了重新投影矢量和栅格对象，并在章节中修改了地图投影。

坐标参考系统

大多数现代地理工具都需要CRS转换，包括核心R空间包和桌面GIS软件如QGIS，它们与 PROJ 接口，这是一个开源C++库，用于“将坐标从一个坐标参考系统（CRS）转换到另一个”。CRS可以用多种方式描述，包括以下几种。

  1. 简单但可能含糊不清的陈述，如“它在经/纬坐标中”
  2. 正式化但现已过时的“proj4字符串”（也称为“proj-string”），如 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
  3. 使用标识的’authority:code’文本字符串，如 EPSG:4326

每个都指的是同一件事：构成全球定位系统（GPS）坐标和许多其他数据集基础的“WGS84”坐标系统。但哪一个才是正确的？

简短的回答是，识别CRSs的第三种方式是首选：EPSG:4326 由本书涵盖的 sf（以及扩展的 stars）和 terra 包以及许多其他用于处理地理数据的软件项目所理解，包括 QGIS 和 PROJ。EPSG:4326 在未来可靠的。此外，尽管它是机器可读的，但与proj-string表示法不同，“EPSG:4326” 简短、易记且在网上非常容易找到（例如，搜索EPSG:4326会在网站 epsg.io 上找到专 dedicated 页）。更简洁的标识符4326由sf理解，但我们建议使用更明确的 AUTHORITY:CODE 表示法来消除歧义并提供上下文。

较长的回答是，这三个描述都不够，并且需要更多的细节来进行明确的CRS处理和转换：由于CRS的复杂性，不可能在如此短的文本字符串中捕获有关它们的所有相关信息。出于这个原因，开放地理空间联盟（OGC，还开发了sf 包实现的简单要素规范）开发了一种用于描述CRS的开放标准格式，称为WKT（众所周知的文本）。这在100多页的文档中详细介绍了“定义了坐标参考系统的抽象模型的文本字符串实现的结构和内容，该抽象模型描述在ISO 19111:2019”中。WGS84 CRS的WKT表示法具有标识符EPSG:4326，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
st_crs("EPSG:4326")
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#>     ENSEMBLE["World Geodetic System 1984 ensemble",
#>         MEMBER["World Geodetic System 1984 (Transit)"],
#>         MEMBER["World Geodetic System 1984 (G730)"],
#>         MEMBER["World Geodetic System 1984 (G873)"],
#>         MEMBER["World Geodetic System 1984 (G1150)"],
#>         MEMBER["World Geodetic System 1984 (G1674)"],
#>         MEMBER["World Geodetic System 1984 (G1762)"],
#>         MEMBER["World Geodetic System 1984 (G2139)"],
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ENSEMBLEACCURACY[2.0]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["Horizontal component of 3D system."],
#>         AREA["World."],
#>         BBOX[-90,-180,90,180]],
#>     ID["EPSG",4326]]

命令的输出显示了CRS标识符（也称为空间参考标识符或SRID）的工作方式：它只是一个查找，提供与CRS的更完整WKT表示形式关联的唯一标识符。这提出了一个问题：如果标识符和CRS的较长WKT表示之间存在不匹配怎么办？关于这一点，@opengeospatialconsortium_wellknown_2019很明确，详细的WKT表示优先于identifier：

📌如果引用的标识符中给出的任何属性或值与WKT描述中明确给出的属性或值冲突，则WKT值应优先。

以AUTHORITY:CODE形式引用CRS标识符的惯例，也被用于其他languages编写的地理软件，允许广泛引用正式定义的坐标系统。^[可以使用几种其他方法来引用唯一的CRS，QGIS和其他标识符类型（例如EPSG:4326标识符的更详细变体，urn:ogc:def:crs:EPSG::4326）接受五种标识符类型（EPSG代码，PostGIS SRID，INTERNAL SRID，proj-string和WKT字符串）。]CRS标识符中最常用的权威机构是EPSG，这是欧洲石油调查小组的首字母缩写，该小组发布了CRS的标准化列表（EPSG由油气机构Geomatics Committee of the International Association of Oil & Gas Producers在2005年taken over）。CRS标识符中可以使用其他权威机构。
例如，ESRI:54030指的是ESRI对罗宾逊投影的实现，其具有以下WKT字符串（仅显示前8行）：

1
2
3
4
5
6
7
8
9
10
st_crs("ESRI:54030")
#> Coordinate Reference System:
#>   User input: ESRI:54030 
#>   wkt:
#> PROJCRS["World_Robinson",
#>     BASEGEOGCRS["WGS 84",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]]],
#> ...

WKT字符串详尽、精确且准确，允许无歧义地存储和转换CRS。它们包含了任何给定CRS的所有相关信息，包括其基准面和椭球体、本初子午线、投影和单位。^[在WKTCRS定义出现之前，proj-string是指定坐标操作和存储CRS的标准方式。这些基于键=值形式（例如，+proj=longlat +datum=WGS84+no_defs）构建的字符串表示，在大多数情况下已经被或者将来应该被WKT表示所取代。]

最近的PROJ版本（6+）仍允许使用proj-string来定义坐标操作，但一些proj-string键（如+nadgrids、+towgs84、+k、+init=epsg:）要么不再受支持，要么不鼓励使用。
此外，只有三个基准面（即WGS84、NAD83和NAD27）可以在proj-string中直接设置。关于CRS定义和PROJ库演变的更长解释可以在@bivand_progress_2021，@pebesma_spatial_2022的第2章，以及Floris Vanderhaeghe的博客文章中找到。此外，如PROJ文档所概述的，WKT CRS格式有不同的版本，包括WKT1和WKT2的两个变体，其中后者（WKT2，2018规范）对应于ISO 19111:2019。

查询和设置坐标系

让我们看看CRS是如何存储在R的空间对象中的，以及如何查询和设置它们。首先，我们将研究如何在矢量地理数据对象中获取和设置CRS，从以下示例开始：

1
2
vector_filepath = system.file("shapes/world.gpkg", package = "spData")
new_vector = read_sf(vector_filepath)

我们的新对象new_vector是一个sf类的数据框，代表了全球各个国家（有关详细信息，请参阅帮助页面?spData::world）。可以使用sf函数st_crs()来检索CRS。

1
2
3
4
5
st_crs(new_vector) # get CRS
#> Coordinate Reference System:
#>   User input: WGS 84 
#>   wkt:
#>   ...

输出是一个包含两个主要组件的列表：

  1. User input（在此情况下为WGS 84，即EPSG:4326的同义词，它在这种情况下是从输入文件中获取的），对应于上述的CRS标识符
  2. wkt，包含有关CRS的所有相关信息的完整WKT字符串

input元素是灵活的，根据输入文件或用户输入，可以包含AUTHORITY:CODE表示（例如EPSG:4326）、CRS的名称（例如WGS 84）或甚至proj-string定义。wkt元素存储WKT表示，用于保存对象到文件或进行任何坐标操作。以上，我们可以看到new_vector对象具有WGS84椭球体，使用格林威治基准经线，以及纬度和经度轴的顺序。在这种情况下，我们还有一些附加元素，如USAGE解释了使用该CRS的适合区域，以及ID指向CRS的标识符：EPSG:4326。

st_crs函数还具有一个有用的功能——我们可以检索有关所用CRS的一些附加信息。例如，尝试运行：

  • st_crs(new_vector)$IsGeographic来检查CRS是否是地理的
  • st_crs(new_vector)$units_gdal找出CRS的单位
  • st_crs(new_vector)$srid提取其’SRID’标识符（如果可用）
  • st_crs(new_vector)$proj4string提取proj-string表示

在坐标参考系统（CRS）缺失或设置错误的CRS的情况下，可以使用st_set_crs()函数（在这种情况下，由于文件在读入时已经正确设置了CRS，WKT字符串保持不变）：

1
new_vector = st_set_crs(new_vector, "EPSG:4326") # set CRS

在栅格地理数据对象中获取和设置坐标参考系统（CRS）的工作方式与矢量数据对象类似。terra包中的crs()函数从SpatRaster对象访问CRS信息（注意使用cat()函数以便漂亮地打印）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
my_rast = rast(raster_filepath)
cat(crs(my_rast)) # get CRS
#> GEOGCRS["WGS 84",
#>     ENSEMBLE["World Geodetic System 1984 ensemble",
#>         MEMBER["World Geodetic System 1984 (Transit)"],
#>         MEMBER["World Geodetic System 1984 (G730)"],
#>         MEMBER["World Geodetic System 1984 (G873)"],
#>         MEMBER["World Geodetic System 1984 (G1150)"],
#>         MEMBER["World Geodetic System 1984 (G1674)"],
#>         MEMBER["World Geodetic System 1984 (G1762)"],
#>         MEMBER["World Geodetic System 1984 (G2139)"],
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ENSEMBLEACCURACY[2.0]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["Horizontal component of 3D system."],
#>         AREA["World."],
#>         BBOX[-90,-180,90,180]],
#>     ID["EPSG",4326]]

输出结果为代表CRS的WKT字符串。相同的函数crs()也可以用来给栅格对象设置CRS。

1
crs(my_rast) = "EPSG:26912" # set CRS

在这里，我们可以使用标识符（在大多数情况下推荐）或完整的WKT字符串表示。设置crs的替代方法包括proj-string字符串或从其他现有对象中提取的CRS，尽管这些方法可能在将来的适用性上较差。

重要的是，st_crs()和crs()函数不会改变坐标的值或几何图形。它们的作用仅仅是设置关于对象CRS的元数据信息。

在某些情况下，地理对象的CRS是未知的，就像下面的代码块中创建的london数据集一样，引入的伦敦示例中就有这样的情况：

1
2
3
4
london = data.frame(lon = -0.1, lat = 51.5) |> 
  st_as_sf(coords = c("lon", "lat"))
st_is_longlat(london)
#> [1] NA

输出的NA显示sf不知道CRS是什么，并且不愿意猜测（NA字面意思是’不可用’）。除非人为指定了CRS或者从具有CRS元数据的源中加载，否则sf不会对坐标系统做出任何明确的假设，除了说“我不知道”。鉴于可用CRS的多样性，这种行为是有道理的，但与一些方法不同，例如GeoJSON文件格式规范，它简化地假设所有坐标都具有经纬度坐标系：EPSG:4326。没有指定CRS的数据集可能会引起问题：所有地理坐标都有坐标系统，只有知道它正在处理的CRS类型，软件才能在绘图和几何操作方面做出正确的决策。因此，再次强调，始终检查数据集的CRS并设置它（如果它丢失）是非常重要的。

1
2
3
london_geo = st_set_crs(london, "EPSG:4326")
st_is_longlat(london_geo)
#> [1] TRUE

投影和未投影数据上的几何操作

自 sf 版本 1.0.0 以来，得益于引入的S2 球面几何引擎，R处理具有经纬度CRS的地理矢量数据集的能力有了很大的提高。如图所示，sf根据CRS的类型和S2是否被禁用（默认启用）来选择使用 GEOS 或 S2。对于投影数据和没有 CRS 的数据，总是使用GEOS；对于地理数据，S2 默认启用，但可以通过 sf::sf_use_s2(FALSE) 禁用。


The behavior of the geometry operations in the sf package depending on the input data’s CRS.

为了展示 CRS 的重要性，我们将在前一节中的london对象周围创建 100 公里的缓冲区。我们还将故意创建一个“距离”为1度的有错误的缓冲区，这大约相当于100公里（赤道上1度大约为111公里）。在深入代码之前，值得简短地跳到下图来直观了解您应该能够通过下面的代码块重现的输出。

第一阶段是围绕上面创建的 london 和 london_geo 对象创建三个缓冲区，边界距离分别为 1 度和 100 公里（或 100,000 米，可以用科学记数法表示为 1e5）距离伦敦市中心。

1
2
3
london_buff_no_crs = st_buffer(london, dist = 1)   # incorrect: no CRS
london_buff_s2 = st_buffer(london_geo, dist = 100000) # silent use of s2
london_buff_s2_100_cells = st_buffer(london_geo, dist = 100000, max_cells = 100)

在上述的第一行中，sf假设输入是投影的，并生成了以度为单位的缓冲区的结果，这是有问题的，我们将看到。在第二行中，sf使用球面几何引擎S2（在章节@ref(spatial-class) 中引入）默默地计算缓冲区的范围，使用 max_cells = 1000 的默认值——在第三行中设置为 100——这些后果将很快显现（详见 ?s2::s2_buffer_cells）。为了强调 sf 对未投影（地理）坐标系统的S2几何引擎的影响，我们将在下面的代码块中暂时禁用它，使用命令sf_use_s2()（默认为 TRUE）。与 london_buff_no_crs一样，新的london_geo对象是一个地理怪物：它的单位是度，这在绝大多数情况下都没有意义：

1
2
3
4
5
6
7
8
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
london_buff_lonlat = st_buffer(london_geo, dist = 1) # incorrect result
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
sf::sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on

上述的警告信息暗示了在经纬度数据上执行平面几何操作的问题。当通过命令 sf::sf_use_s2(FALSE)关闭球面几何操作时，缓冲区（以及其他几何操作）可能会产生无价值的输出，因为它们使用的是纬度和经度的单位，这与以米为单位的适当距离单位相比是一个糟糕的替代品。

📌两条经线（称为子午线）之间的距离在赤道处约为111公里（执行geosphere::distGeo(c(0,0),c(1,0))找到精确距离）。这在极地缩小为零。例如，在伦敦的纬度上，子午线之间的距离不到70公里（挑战：执行验证此事的代码）。

相比之下，纬线彼此之间的距离不受纬度的影响：它们总是相距约111公里，包括在赤道和靠近极地的地方。

不要将关于地理（经度/纬度）坐标参照系的警告解释为“坐标参照系不应该被设置”：它几乎总是应该被设置的！更好的理解是建议将数据重新投影到投影坐标参照系上。此建议并不总是需要遵循：在某些情况下，执行空间和几何操作几乎没有差别（例如，空间子集）。但是对于涉及距离的操作，例如缓冲区，确保良好结果的唯一方法（不使用球形几何引擎）是创建数据的投影副本，并在其上运行操作。在下面的代码块中完成了这一操作：

1
2
london_proj = data.frame(x = 530000, y = 180000) |> 
  st_as_sf(coords = c("x", "y"), crs = "EPSG:27700")

结果是一个与london相同但重新投影到合适坐标参照系（在这种情况下为英国国家网格，其EPSG代码为27700）的新对象，该坐标参照系的单位为米。我们可以使用st_crs()来验证CRS已更改，如下所示（部分输出已被替换为...）：

1
2
3
4
5
6
7
8
9
10
st_crs(london_proj)
#> Coordinate Reference System:
#>   User input: EPSG:27700 
#>   wkt:
#> PROJCRS["OSGB36 / British National Grid",
#>     BASEGEOGCRS["OSGB36",
#>         DATUM["Ordnance Survey of Great Britain 1936",
#>             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
#>                 LENGTHUNIT["metre",1]]],
#> ...

此CRS描述的显著组成部分包括EPSG代码（EPSG:27700）和详细的wkt字符串（仅显示前5行）。^[有关最相关投影参数和相关概念的简短描述，请参阅Jochen Albrecht在http://www.geography.hunter.cuny.edu/~jochen/GTECH361/lectures/ 主持的第四讲，并在https://proj.org/usage/projections.html 查看信息。关于投影的其他优秀资源是spatialreference.org和progonos.com/furuti/MapProj。]CRS的单位描述在LENGTHUNIT字段中，单位是米（而不是度），这告诉我们这是一个投影的CRS：st_is_longlat(london_proj)现在返回FALSE，并且对london_proj的几何操作将无需警告即可工作。对london_proj的缓冲区操作将使用GEOS，结果将以适当的距离单位返回。以下代码行在投影数据周围创建了恰好100公里的缓冲区：

1
london_buff_projected = st_buffer(london_proj, 100000)

上述三个具有指定CRS的london_buff*对象（london_buff_s2，london_buff_lonlat和london_buff_projected）在前面的代码块中创建的几何图形在下图中有所展示。

1
2
3
4
5
6
#> Warning: The `returnclass` argument of `ne_download()` sp as of rnaturalearth 1.0.0.
#> ℹ Please use `sf` objects with {rnaturalearth}, support for Spatial objects
#>   (sp) will be removed in a future release of the package.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.


Buffers around London showing results created with the S2 spherical geometry engine on lon/lat data (left), projected data (middle) and lon/lat data without using spherical geometry (right). The left plot illustrates the result of buffering unprojected data with sf, which calls Google’s S2 spherical geometry engine by default with max cells set to 1000 (thin line). The thick ‘blocky’ line illustrates the result of the same operation with max cells set to 100.

从上图中可以清楚地看到，基于s2和适当投影的CRS的缓冲区并没有被“压扁”，意味着缓冲区边界的每一部分都与伦敦等距离。当s2未被使用时，从经纬度CRS生成的结果，无论是因为输入缺少CRS还是因为sf_use_s2()被关闭，都会严重失真，结果在南北轴向延伸，突出了在经纬度输入上使用假设投影数据的算法（如GEOS所做）的危险性。然而，使用S2生成的结果也同样失真，尽管不太显著。图中（左）的两个缓冲区边界都是锯齿状的，尽管这可能只在用s2参数max_cells设置为100创建的粗边界表示的缓冲区时才显现或相关。

教训是，通过S2从经纬度数据获得的结果将与使用投影数据获得的结果不同。随着max_cells值的增加，S2派生的缓冲区和GEOS派生的投影数据上的缓冲区之间的差异减小，此参数的“正确”值可能取决于许多因素，缺省值1000是合理的。在选择max_cells值时，应在计算速度和结果分辨率之间取得平衡。在曲线边界有优势的情况下，在缓冲（或执行其他几何操作）之前转换到投影CRS可能是合适的。

从上述示例中可以清楚地看到CRS的重要性（主要是它们是投影的还是地理的）以及sf默认使用S2在经纬度数据上进行缓冲的影响。后续章节将更深入地探讨何时需要投影CRS以及重投影矢量和栅格对象的细节。

何时重投影？

上一节展示了如何手动设置坐标参考系统（CRS），通过st_set_crs(london, "EPSG:4326")。然而，在现实世界的应用中，CRS通常在数据读入时自动设置。在许多项目中，主要与CRS有关的任务是将对象从一个CRS转换为另一个CRS。但是，何时应该转换数据？转换为哪个CRS？这些问题没有明确的答案，CRS的选择总是涉及权衡。然而，本节提供了一些一般性原则，可以帮助您做决策。

首先值得考虑的是何时进行转换。在某些情况下，转换为地理CRS是至关重要的，例如使用leaflet包在线发布数据时。另一个情况是当必须比较或组合两个具有不同CRS的对象时，例如当我们试图找到具有不同CRS的两个对象之间的距离时：

1
2
st_distance(london_geo, london_proj)
# > Error: st_crs(x) == st_crs(y) is not TRUE

要使london和london_proj对象在地理上可比较，必须将其中一个转换为另一个的CRS。但是应该使用哪个CRS呢？答案取决于具体情况：许多项目，特别是涉及网络映射的项目，需要以EPSG:4326的输出，此时最好转换投影对象。然而，如果项目需要平面几何操作而不是球面几何操作引擎（例如，为了创建具有平滑边缘的缓冲区），那么将具有地理CRS的数据转换为具有投影CRS的等效对象可能是值得的，例如英国国家网格（EPSG:27700）。

选择哪个CRS使用?

选择哪个CRS的问题很棘手，很少有“正确”的答案：“不存在全能的投影，当远离指定框架的中心时，所有投影都会发生失真”。此外，你不应该只依附于一个投影来完成每项任务。有可能在分析的某些部分使用一种投影，在另一部分使用另一种投影，甚至在可视化中使用其他投影。始终尝试选择最能实现你目标的CRS！在选择**地理坐标参照系（CRS）**时，答案通常是WGS84。它不仅用于网络地图绘制，而且由于GPS数据集以及数以千计的栅格和矢量数据集默认都提供在这个CRS下。WGS84是世界上最常用的CRS，因此值得记住其EPSG代码：4326。这个“神奇的数字”可以用来将具有不寻常的投影CRS的对象转换成广泛理解的形式。

那么在需要**投影坐标参照系（CRS）**时怎么办呢？在某些情况下，这不是我们可以自由决定的问题：“通常投影的选择是由公共测绘机构做出的”这意味着当使用本地数据源时，最好使用数据提供的CRS来工作，以确保兼容性，即使官方的CRS不是最精确的。伦敦的例子很容易回答，因为（a）英国国家网格（及其相关的EPSG代码27700）众所周知，（b）原始数据集（london）已经有了该CRS。

常用的默认值是通用横轴墨卡托（UTM），一组将地球划分为60个经度楔形和20个纬度段的CRS。UTM CRS所使用的横轴墨卡托投影是共形的，但随着距离UTM区域中心的距离增加，面积和距离的失真程度逐渐加剧。因此，GIS软件Manifold的文档建议将使用UTM区域的项目的经度范围限制在中央子午线的6度之内（来源：manifold.net）。所以，我们建议仅在您关注保持相对较小区域的角度时使用UTM！

地球上几乎每个地方都有一个UTM代码，例如“60H”，指的是新西兰北部R语言发明的地方。UTM的EPSG代码对于北半球的位置依次从32601到32660，对于南半球的位置则从32701到32760。

要展示该系统是如何工作的，让我们创建一个函数lonlat2UTM()，用来计算地球上任何点所对应的EPSG代码，如下所示：

1
2
3
4
5
6
7
8
lonlat2UTM = function(lonlat) {
  utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1
  if (lonlat[2] > 0) {
    utm + 32600
  } else{
    utm + 32700
  }
}

以下命令使用这个函数来识别奥克兰和伦敦的UTM区域和相应的EPSG代码。这些代码对于进一步的地理空间数据分析和操作非常重要。

1
2
3
4
lonlat2UTM(c(174.7, -36.9))
#> [1] 32760
lonlat2UTM(st_coordinates(london))
#> [1] 32630

目前，我们也有一些工具可以帮助我们选择合适的CRS（坐标参考系统），其中包括crssuggest包。该包的主要函数 suggest_crs()接受一个带有地理CRS的空间对象，并返回可能用于给定区域的投影CRS列表。^[该包还允许在没有任何CRS信息附加的情况下找出数据的真实CRS。]另一个有用的工具是一个网页https://jjimenezshaw.github.io/crs-explorer/，它根据所选的位置和类型列出CRS。

重要提示：虽然这些工具在许多情况下都很有帮助，但在应用推荐的CRS之前，你需要了解其属性。

当不立即清楚合适的CRS时，CRS的选择应取决于后续地图和分析中最重要的要保留的属性。所有的CRS要么是等面积的，要么是等距的，要么是保形的（形状保持不变），要么是这些特性的某种折中组合。可以为感兴趣的区域创建带有本地参数的自定义CRS，并且当单个CRS不适用于所有任务时，项目中可以使用多个CRS。如果没有合适的CRS，'测地计算’可以提供一个备选方案（参见proj.org/geodesic.html）。无论使用哪个投影CRS，对于覆盖数百公里的几何体，结果可能不准确。

在决定定制CRS时，我们建议如下：^[非常感谢一位匿名审稿人的评论，他们构成了这些建议的基础。]

  • 对于定制的本地投影（将原点的纬度和经度设置为研究区域的中心），可以使用Lambert方位等面积（LAEA）投影，这是一种在所有位置的等面积投影，但在数千公里以外会扭曲形状。
  • 对于特定准确的直线距离，可以使用方位等距（AEQD）投影，该投影的起点与本地投影的中心点之间的距离。
  • 对于覆盖数千公里的区域，可以使用Lambert保角圆锥（LCC）投影，其锥体设置可在割线之间保持距离和面积属性合理。
  • 对于极地区域，可以使用立体（STERE）投影，但要注意不要依赖于离中心数千公里的面积和距离计算。

一种可能的方法是为研究区域的中心点创建方位等距（AEQD）投影，以自动选择特定于本地数据集的投影坐标参考系统（CRS）。这涉及基于数据集的中心点，使用米为单位创建自定义CRS（没有EPSG代码）。请注意，应谨慎使用此方法：没有其他数据集与创建的自定义CRS兼容，当用于覆盖数百公里的广泛数据集时，结果可能不准确。

本节概述的原则同样适用于矢量和栅格数据集。然而，CRS转换的某些特点是每个地理数据模型所独有的。我们将介绍矢量数据转换的特殊性和栅格转换的特殊性。接下来，最后一节将展示如何创建自定义地图投影。

矢量几何重投影

空间类章演示了矢量几何是如何由点组成的，以及点如何形成更复杂的对象，如线和多边形。因此，重新投影矢量包括转换这些点的坐标，它们构成了线和多边形的顶点。

何时重投影节包含了一个示例，其中至少一个sf对象必须转换为具有不同CRS的等效对象，以计算两个对象之间的距离。

1
london2 = st_transform(london_geo, "EPSG:27700")

现在已经使用 sf 函数 st_transform() 创建了 london 的转换版本，就可以找到两个伦敦表示之间的距离了。^[替代 st_transform()的是lwgeom中的st_transform_proj()，它可以进行绕过GDAL的转换，并可以支持GDAL不支持的投影。在撰写本文时（2022年），我们没有找到任何仅由st_transform_proj() 支持而不由 st_transform() 支持的投影。]令人惊讶的是，london 和 london2 相距超过2公里！^[两点之间的位置差异并不是由于转换操作的不完美（实际上非常精确），而是由于手动创建的london和london_proj的坐标精度低所造成的。同样令人惊讶的是，结果以米为单位的矩阵形式提供。这是因为 st_distance()可以提供许多特征之间的距离，而且CRS的单位是米。使用 as.numeric() 将结果强制转换为常规数字。]

1
2
3
4
st_distance(london2, london_proj)
#> Units: [m]
#>      [,1]
#> [1,] 2018

下面展示了使用 cycle_hire_osm（来自 spData 的一个 sf 对象，代表伦敦可以租借自行车的“停靠站”）查询和重新投影CRS的函数。可以用函数 st_crs() 查询sf对象的CRS。输出以包含坐标系统信息的多行文本打印：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
st_crs(cycle_hire_osm)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCS["WGS 84",
#>     DATUM["WGS_1984",
#>         SPHEROID["WGS 84",6378137,298.257223563,
#>             AUTHORITY["EPSG","7030"]],
#>         AUTHORITY["EPSG","6326"]],
#>     PRIMEM["Greenwich",0,
#>         AUTHORITY["EPSG","8901"]],
#>     UNIT["degree",0.0174532925199433,
#>         AUTHORITY["EPSG","9122"]],
#>     AUTHORITY["EPSG","4326"]]

主要的CRS组件 User input 和 wkt被打印为单一实体，st_crs()的输出实际上是一个名为crs的类的命名列表，具有两个元素，分别命名为 input 和 wkt，如以下代码块的输出所示：

1
2
3
4
5
crs_lnd = st_crs(london_geo)
class(crs_lnd)
#> [1] "crs"
names(crs_lnd)
#> [1] "input" "wkt"

可以使用 $ 运算符检索附加元素，包括 Name、proj4string 和 epsg（有关详细信息，请参阅?st_crs和GDAL网站 上的 CRS 和转换教程）：

1
2
3
4
5
6
crs_lnd$Name
#> [1] "WGS 84"
crs_lnd$proj4string
#> [1] "+proj=longlat +datum=WGS84 +no_defs"
crs_lnd$epsg
#> [1] 4326

WKT 表示形式存储在 crs_lnd 对象的 $wkt 元素中，是最终的真实来源。这意味着前面代码块的输出是来自 PROJ 提供的 wkt 表示的查询，而不是对象及其 CRS 的固有属性。

当对象的CRS被转换时，CRS的wkt和User Input元素都会发生改变。在下面的代码块中，我们用投影 CRS 创建了一个新版本的 cycle_hire_osm（为了简洁，仅显示 CRS 输出的前4行）：

1
2
3
4
5
6
7
cycle_hire_osm_projected = st_transform(cycle_hire_osm, "EPSG:27700")
st_crs(cycle_hire_osm_projected)
#> Coordinate Reference System:
#>   User input: EPSG:27700 
#>   wkt:
#> PROJCRS["OSGB36 / British National Grid",
#> ...

生成的对象具有新的CRS，其EPSG代码为27700。但是如何了解有关此EPSG代码或任何代码的更多详细信息呢？一个选项是在线搜索，

1
2
3
4
5
6
7
crs_lnd_new = st_crs("EPSG:27700")
crs_lnd_new$Name
#> [1] "OSGB36 / British National Grid"
crs_lnd_new$proj4string
#> [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
crs_lnd_new$epsg
#> [1] 27700

结果显示，EPSG代码27700代表了英国国家网格，这个结果可以通过在线搜索“EPSG 27700”找到。

📌在控制台中打印一个空间对象会自动返回其坐标参考系统。要显式地访问和修改它，可以使用st_crs函数，例如st_crs(cycle_hire_osm)。

栅格几何重投影

在前一部分描述的投影概念同样适用于栅格数据。然而，矢量和栅格的重新投影之间存在重要差异：转换矢量对象涉及更改每个顶点的坐标，但这不适用于栅格数据。栅格由相同大小的矩形单元格组成（以地图单位表示，如度或米），因此单独转换像素的坐标通常是不切实际的。栅格重新投影涉及创建一个新的栅格对象，与原始对象相比，通常具有不同数量的列和行。随后必须重新估计属性，从而允许用适当的值“填充”新的像素。换句话说，栅格重新投影可以看作是两个单独的空间操作：栅格范围到另一个CRS的矢量重新投影，以及通过重新采样计算新的像素值。因此，在大多数情况下，当使用栅格和矢量数据时，最好避免重新投影栅格，而是重新投影向量。

📌常规栅格的重新投影也被称为扭曲。
此外，还有一种称为“变换”的类似操作。与重新采样所有值不同，它保留所有值不变，但重新计算每个栅格单元格的新坐标，改变网格几何形状。例如，它可以将输入栅格（一个常规网格）转换为曲线网格。在R中，可以使用stars包执行变换操作。

栅格重新投影过程是使用terra包中的project()函数完成的。与前一节中演示的st_transform()函数类似，project()采用地理对象（在本例中为栅格数据集）和第二个参数中的某些CRS表示。顺便说一下–第二个参数也可以是具有不同CRS的现有栅格对象。

让我们看两个光栅转换的例子：使用分类和连续数据。土地覆盖数据通常由分类地图表示。nlcd.tif文件为美国犹他州的一个小区域提供了信息，该信息来自National Land Cover Database 2011，在NAD83/UTM zone 12 N CRS中，如下面的代码块输出所示（仅显示输出的第一行）：

1
2
3
4
cat_raster = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
crs(cat_raster)
#> PROJCRS["NAD83 / UTM zone 12N",
#> ...

在这个区域内，区分了8个土地覆盖类型（NLCD2011土地覆盖类型的完整列表可以在mrlc.gov上找到）：

1
2
3
4
5
6
7
8
9
10
unique(cat_raster)
#>       levels
#> 1      Water
#> 2  Developed
#> 3     Barren
#> 4     Forest
#> 5  Shrubland
#> 6 Herbaceous
#> 7 Cultivated
#> 8   Wetlands

在重新投影分类光栅时，估计值必须与原始值相同。这可以使用最近邻方法（near）完成，该方法将每个新单元格的值设置为输入光栅的最近单元格（中心）的值。一个例子是将cat_raster重新投影到WGS84，这是一种非常适合网络映射的地理坐标参照系。第一步是获得该CRS的PROJ定义，例如可以使用http://spatialreference.org网页来完成。最后一步是使用project()函数重新投影光栅，在分类数据的情况下，它使用最近邻方法（near）。

1
cat_raster_wgs84 = project(cat_raster, "EPSG:4326", method = "near")

重新投影后的光栅与原始光栅在许多属性上有所不同，包括列数和行数（因此单元格数）、分辨率（从米转换为度）和范围，如表@ref(tab:catraster)所示（请注意，类别数量从8增加到9，这是由于增加了NA值，而不是因为创建了新类别 —— 土地覆盖类别得以保留）。

Table: Key attributes in the original (‘cat_raster’) and projected (‘cat_raster_wgs84’) categorical raster datasets.

CRS nrow ncol ncell resolution unique_categories
NAD83 1359 1073 1458207 31.5275 8
WGS84 1246 1244 1550024 0.0003 9

对数值光栅（具有 numeric 或在此例中的 integer 值）的重新投影过程几乎相同。下面使用 spDataLarge 中的 srtm.tif 文件展示这一过程，该文件来自 Shuttle Radar Topography Mission（SRTM），以WGS84坐标参考系表示海平面以上的米高（海拔）：

1
2
3
con_raster = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
crs(con_raster)
#> [1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

我们将重新投影这个数据集到一个投影的坐标参考系中，但不是用适用于分类数据的最近邻方法。相反，我们将使用双线性方法，该方法基于原始栅格中的四个最近单元格计算输出单元格的值。^[
还可以在重采样节中提到的其他方法。]投影数据集中的值是这四个单元格的值的距离加权平均值：输入单元格距离输出单元格的中心越近，其权重就越大。以下命令创建代表WGS 84 / UTM zone 12N的文本字符串，并使用bilinear方法将栅格重新投影到此CRS中。

1
2
3
con_raster_ea = project(con_raster, "EPSG:32612", method = "bilinear")
crs(con_raster_ea)
#> [1] "PROJCRS[\"WGS 84 / UTM zone 12N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 12N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-111,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA).\"],\n        BBOX[0,-114,84,-108]],\n    ID[\"EPSG\",32612]]"

数值变量的栅格重新投影也会导致数值和空间属性的微小变化，例如单元格数量、分辨率和范围。这些变化在表 @ref(tab:rastercrs) 中得到展示。^[另一个未在表中表示的细微变化是新投影栅格数据集中值的类别是 numeric。这是因为bilinear方法适用于连续数据，结果很少强制转换为整数值。这可能会在保存栅格数据集时影响文件大小。]

Table: (#tab:rastercrs)Key attributes in the original (‘con_raster’) and projected (‘con_raster_ea’) continuous raster datasets.

CRS nrow ncol ncell resolution mean
WGS84 457 465 212505 0.0008 1843
UTM zone 12N 515 422 217330 83.5334 1842

📌当然，2D地球投影的局限性同样适用于矢量数据和栅格数据。我们最多只能遵循三个空间属性（距离、面积、方向）中的两个。因此，手头的任务决定了选择哪个投影。例如，如果我们对密度感兴趣（每个网格单元的点数或每个网格单元的居民数），我们应该使用等面积投影。

自定义坐标系

使用AUTHORITY:CODE标识符（例如EPSG:4326）捕获的现有CRS非常适合许多应用。然而，在某些情况下，可能希望使用替代投影或创建自定义CRS。*选择哪个CRS使用?*节提到了使用自定义CRS的原因，并提供了几种可能的方法。在这里，我们展示如何在R中应用这些想法。

其中一种方法是采用现有的CRS的WKT定义，修改其中的一些元素，然后使用新定义进行重新投影。对于空间向量，可以使用st_crs()$wkt和st_transform()，对于空间栅格，可以使用crs()和project()，如下面的示例所示，该示例将zion对象转换为自定义的方位等距（AEQD）CRS。

1
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))

使用自定义AEQD CRS需要知道数据集中心点的坐标（以度为单位，地理CRS）。在我们的情况下，可以通过计算zion区域的质心并将其转换为WGS84来提取这些信息。

1
2
3
4
zion_centr = st_centroid(zion)
zion_centr_wgs84 = st_transform(zion_centr, "EPSG:4326")
st_as_text(st_geometry(zion_centr_wgs84))
#> [1] "POINT (-113 37.3)"

接下来，我们可以使用新获得的值来更新下面所示的方位等距（AEQD）CRS的WKT定义。注意，我们只修改了下面的两个值 —— 将“Central_Meridian”改为我们的质心的经度，并将“Latitude_Of_Origin”改为质心的纬度。

1
2
3
4
5
6
7
8
9
10
my_wkt = 'PROJCS["Custom_AEQD",
 GEOGCS["GCS_WGS_1984",
  DATUM["WGS_1984",
   SPHEROID["WGS_1984",6378137.0,298.257223563]],
  PRIMEM["Greenwich",0.0],
  UNIT["Degree",0.0174532925199433]],
 PROJECTION["Azimuthal_Equidistant"],
 PARAMETER["Central_Meridian",-113.0263],
 PARAMETER["Latitude_Of_Origin",37.29818],
 UNIT["Meter",1.0]]'

这种方法的最后一步是将我们的原始对象（zion）转换为我们的新自定义CRS（zion_aeqd）。

1
zion_aeqd = st_transform(zion, my_wkt)

自定义投影也可以交互式地制作，例如，使用Projection Wizard网页应用程序。这个网站允许您选择数据的空间范围和失真属性，并返回可能的投影列表。该列表还包含您可以复制和用于重新投影的投影的WKT定义。有关使用WKT字符串创建自定义CRS定义的详细信息，请参阅@opengeospatialconsortium_wellknown_2019。

PROJ字符串也可以用来创建自定义投影，接受投影本身固有的限制，特别是在涵盖大地理区域的几何体上。已经开发了许多投影，并可以通过PROJ字符串的+proj=元素进行设置，仅在PROJ网站上就详细描述了几十个项目。

当绘制世界地图并保持面积关系时，Mollweide投影是一种流行且通常合理的选择，如下图所示。要使用此投影，我们需要在st_transform函数中使用proj-string元素"+proj=moll"来指定它。

1
world_mollweide = st_transform(world, crs = "+proj=moll")

<
Mollweide projection of the world.

当绘制世界地图时，人们通常希望最小化所有空间属性（面积、方向、距离）的失真。实现这一目标的最受欢迎的投影之一是Winkel tripel，如下图所示。^[该投影被国家地理学会等机构使用。]这个结果是通过以下命令创建的：

1
world_wintri = st_transform(world, crs = "+proj=wintri")


Winkel tripel projection of the world.

此外，在大多数CRS定义中可以修改proj-string参数，例如可以使用+lon_0和+lat_0参数调整投影的中心点。下面的代码将坐标转换为以纽约市的经度和纬度为中心的兰伯特等面积投影。

1
2
world_laea2 = st_transform(world,
                           crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=-74 +lat_0=40")


Lambert azimuthal equal-area projection of the world centered on New York City.

关于CRS修改的更多信息可以在使用PROJ文档中找到。

练习

E1. 通过将 nz 对象转换为 WGS84 坐标参考系统（CRS）来创建一个名为 nz_wgs 的新对象。

  • 为两者创建一个 crs 类的对象，并用它来查询它们的CRS。
  • 根据每个对象的边界框，每个CRS使用什么单位？
  • 从 nz_wgs 中移除CRS并绘制结果：这张新西兰的地图有什么问题，为什么？

E2. 将 world 数据集转换为横向墨卡托投影（"+proj=tmerc"）并绘制结果。发生了什么变化，为什么？尝试将其转换回 WGS 84 并绘制新对象。为什么新对象与原对象不同？

E3. 使用最近邻插值方法将连续栅格（con_raster）转换为 NAD83 / UTM区域12N。发生了什么变化？它如何影响结果？

E4. 使用双线性插值方法将分类栅格（cat_raster）转换为 WGS 84。发生了什么变化？它如何影响结果？

-------------已经到底啦-------------