(8)地理数据 I/O

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

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

本章关于读取和写入地理数据。地理数据输入对地理计算至关重要，没有数据就不可能有真实世界的应用。数据输出也非常关键,它能够让其他人使用由你的工作产生的有价值的新数据集或改进后的数据集。总而言之,输入/输出这些进程可以称为数据I/O。

前提条件

本章要求下面的包

1
2
3
4
library(sf)
library(terra)
library(dplyr)
library(spData)

引言

本章关于读取和写入地理数据。地理数据输入对地理计算至关重要，没有数据就不可能有真实世界的应用。数据输出也非常关键，它能够让其他人使用由你的工作产生的有价值的新数据集或改进后的数据集。总而言之,输入/输出这些进程可以称为数据I/O。

地理数据I/O通常只需要在项目开始和结束时编写几行代码。它经常被简单地视为一个步骤的过程而被忽略。然而,在项目开始时(例如使用过时或有缺陷的数据集)的错误会在稍后导致大问题,因此值得在确定哪些数据集是可用的、它们可以在哪里找到以及如何检索它们上投入相当多的时间。后面章节描述了各种地理门户网站，这些网站共同包含了数十TB的数据,以及如何使用它们。为了进一步简化数据访问,已经开发了许多用于下载地理数据的软件包。存在许多地理文件格式,各有其优缺点。关于实际高效读取和写入这些文件格式的过程直到第数据输入节和数据输出节才涵盖。最后一节输出可视化演示了保存可视化输出(地图)的方法，为高级制图章节关于可视化做准备。

获取开放数据

在互联网上可以获取到大量且不断增加的地理数据，其中许多数据是免费提供并可供访问和使用的（需适当给予提供者适当的署名）。[^08-read-write-plot-1]从某种程度上说，现在存在着过多的数据，意味着通常有多个地方可以获取相同的数据集。有些数据集的质量较差。在这种情况下，了解从何处获取数据变得至关重要，因此第一部分介绍了一些最重要的数据来源。各种“地理门户”（提供地理空间数据的网络服务，如Data.gov）是一个很好的起点，提供了广泛的数据，但通常仅适用于特定位置（如在更新的维基百科页面上所示）。

[^08-read-write-plot-1]：例如，访问https://freegisdata.rtwilson.com/获取免费可用的地理数据网站列表。

一些全球性的地理门户解决了这个问题。例如，GEOSS门户和科普尼库斯开放获取中心包含许多具有全球覆盖范围的栅格数据集。从由美国国家航空航天局（NASA）运营的SEDAC门户和欧洲联盟的INSPIRE地理门户可以访问到丰富的矢量数据集，涵盖全球和区域。

大多数地理门户提供图形界面，允许基于空间和时间范围等特征查询数据集，美国地质调查局的EarthExplorer就是一个典型的例子。通过浏览器交互式地探索数据集是了解可用图层的有效方式。然而，从可重现性和效率的角度来看，下载数据最好使用代码进行。可以通过多种技术从命令行启动下载，主要是通过URL和API（例如查看哨兵API）。托管在静态URL上的文件可以使用download.file()进行下载，正如下面的代码块所示，该代码块访问了来自doi.pangaea.de的PeRL：多年冻土区域池塘和湖泊数据库：

1
2
3
4
5
download.file(url = "https://hs.pangaea.de/Maps/PeRL/PeRL_permafrost_landscapes.zip",
              destfile = "PeRL_permafrost_landscapes.zip", 
              mode = "wb")
unzip("PeRL_permafrost_landscapes.zip")
canada_perma_land = read_sf("PeRL_permafrost_landscapes/canada_perma_land.shp")

地理数据包

已经开发了许多用于访问地理数据的R包，其中一些在表中呈现出来。这些包提供了一个或多个空间库或地理门户的接口，旨在使命令行中的数据访问更加快速。

(\#tab:datapackages)Selected R packages for geographic data retrieval.
Package Description
FedData Datasets maintained by the US Federal government, including elevation and land cover.
geodata Download and import imports administrative, elevation, WorldClim data.
osmdata Download and import small OpenStreetMap datasets.
osmextract Download and import large OpenStreetMap datasets.
rnaturalearth Access to Natural Earth vector and raster data.
rnoaa Imports National Oceanic and Atmospheric Administration (NOAA) climate data.

需要强调的是，表中仅代表了少量可用的地理数据包。例如，存在大量的R包用于获取各种社会人口数据，如tidycensus和tigris（美国），cancensus（加拿大），eurostat和giscoR（欧盟），或者idbr（国际数据库）–阅读分析美国人口普查数据以找到如何分析此类数据的一些示例。同样，还有一些R包可供访问各个地区和国家的空间数据，如bcdata（不列颠哥伦比亚省），geobr（巴西），RCzechia（捷克）或rgugik（波兰）。另一个值得注意的包是GSODR，它在R中提供全球每日天气摘要数据（请参阅该包的README以获取天气数据源的概述）。

每个数据包都有自己的访问数据的语法。这种多样性在随后的代码块中得到了展示，这些代码块展示了如何使用表@ref(tab:datapackages)中的三个包获取数据。^08-read-write-plot-2 国家边界通常很有用，可以使用rnaturalearth包中的ne_countries()函数来访问这些边界，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(rnaturalearth)
#> Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
usa = ne_countries(country = "United States of America") # United States borders
#> 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.
class(usa)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"
# alternative way of accessing the data, with geodata
# geodata::gadm("USA", level = 0, path = tempdir())

默认情况下，rnaturalearth返回Spatial*类的对象。可以使用st_as_sf()将结果转换为sf对象，如下所示：

1
usa_sf = st_as_sf(usa)

第二个示例使用geodata包下载了一系列具有十分钟空间分辨率（在赤道处约为18.5公里）的全球月降水总和栅格数据。结果是一个SpatRaster类的多层对象。

1
2
3
library(geodata)
worldclim_prec = worldclim_global("prec", res = 10, path = tempdir())
class(worldclim_prec)

第三个示例使用osmdata包从OpenStreetMap（OSM）数据库中查找公园。如下面的代码块所示，查询从函数opq()（OpenStreetMap查询的缩写）开始，其第一个参数是边界框或表示边界框的文本字符串（在此示例中为利兹市）。结果传递给一个函数，用于选择我们感兴趣的OSM元素（在本例中为公园），表示为键-值对。接下来，它们被传递给osmdata_sf()函数，该函数完成了下载数据并将其转换为sf对象列表的工作（有关详细信息，请参阅vignette('osmdata')）：

1
2
3
4
library(osmdata)
parks = opq(bbox = "leeds uk") |> 
  add_osm_feature(key = "leisure", value = "park") |> 
  osmdata_sf()

osmdata包的一个限制是它受到速率限制，这意味着它无法下载大型的OSM数据集（例如一个大城市的所有OSM数据）。为了克服这个限制，开发了osmextract包，它可以用来下载和导入包含预定义区域压缩版本的二进制.pbf文件，其中包含了OSM数据库的数据。

OpenStreetMap是一个庞大的全球众包数据数据库，每天都在不断增长，并且有一个更广泛的工具生态系统，可以轻松访问数据，从用于快速开发和测试OSM查询的Overpass turbo网络服务，到将数据导入PostGIS数据库的osm2pgsql。尽管从OSM派生的数据集的质量各不相同，但数据源和更广泛的OSM生态系统具有许多优势：它们提供全球范围内的免费数据集，并且得益于众多志愿者的帮助不断改进。使用OSM鼓励“公民科学”和对数字公共领域的贡献（您可以从www.openstreetmap.org开始编辑代表您熟悉的世界某个部分的数据）。

有时，软件包附带了内置数据集。可以通过四种方式访问这些数据集：通过附加软件包（如果软件包使用spData的“惰性加载”），使用data(dataset, package = mypackage)，通过mypackage::dataset引用数据集，或者使用system.file(filepath, package = mypackage)来访问原始数据文件。以下代码块使用world数据集（通过附加其父包library(spData)已加载）说明了后两种选项：[^08-read-write-plot-3]

[^08-read-write-plot-3]：有关使用R软件包进行数据导入的更多信息，请参见@gillespie_efficient_2016的第5.5和第5.6节。

1
2
world2 = spData::world
world3 = read_sf(system.file("shapes/world.gpkg", package = "spData"))

最后一个示例system.file("shapes/world.gpkg", package = "spData")返回了一个路径，该路径指向spData软件包中的"shapes/"文件夹内的world.gpkg文件。

获取空间信息的另一种方法是进行地理编码——将位置描述（通常是地址）转换为其坐标。这通常通过向在线服务发送查询并获得位置作为结果来完成。存在许多此类服务，它们在地理编码方法、使用限制、成本或API密钥要求方面存在差异。R有几个用于地理编码的包；然而，tidygeocoder似乎允许使用一致的接口连接到最多数量的地理编码服务。tidygeocoder的主要函数是geocode，它接受一个包含地址的数据框，并将坐标添加为"lat"和"long"。该函数还允许使用method参数选择地理编码服务，并具有许多其他参数。

下面的示例搜索位于伦敦苏豪区一幢建筑上的约翰·斯诺蓝色牌匾的坐标。

1
2
3
4
library(tidygeocoder)
geo_df = data.frame(address = "54 Frith St, London W1D 4SJ, UK")
geo_df = geocode(geo_df, address, method = "osm")
geo_df

得到的数据框可以使用st_as_sf()转换为sf对象。

1
geo_sf = st_as_sf(geo_df, coords = c("long", "lat"), crs = "EPSG:4326")

这个包还允许执行相反的过程，称为反向地理编码，用于基于一对坐标获取一组信息（名称、地址等）。

地理网络服务

为了标准化访问空间数据的Web API，开放地理信息联盟（OGC）制定了一系列Web服务规范（统称为OGC Web Services，简称OWS）。这些规范包括Web要素服务（Web Feature Service，WFS）、Web地图服务（Web Map Service，WMS）、Web地图切片服务（Web Map Tile Service，WMTS）、Web覆盖服务（Web Coverage Service，WCS），甚至还有Web处理服务（Web Processing Service，WPS）。像PostGIS等地图服务器已经采用了这些协议，从而实现了查询的标准化。与其他Web API一样，OWS API使用“base URL”、“endpoint”和在?后的“UURL query arguments”来请求数据（请参阅httr软件包中的best-practices-api-packages文章）。

有许多请求可以发送到OWS服务。其中一个最基本的请求是getCapabilities，下面使用httr函数GET()和modify_url()演示了这个请求。下面的代码块演示了如何构建和发送查询，本例中是为了发现由联合国粮食及农业组织（FAO）运行的服务的功能：

1
2
3
4
5
6
7
library(httr)
base_url = "http://www.fao.org"
endpoint = "/figis/geoserver/wfs"
q = list(request = "GetCapabilities")
res = GET(url = modify_url(base_url, path = endpoint), query = q)
res$url
#> [1] "https://www.fao.org/figis/geoserver/wfs?request=GetCapabilities"

上面的代码块演示了如何使用GET()函数以编程方式构建API请求，该函数接受基本URL和查询参数列表，可以轻松扩展。请求的结果保存在res中，它是httr包中定义的response类的对象，是一个包含请求信息的列表，包括URL。通过执行browseURL(res$url)可以看到，结果也可以直接在浏览器中查看。提取请求内容的一种方式如下：

1
2
txt = content(res, "text")
xml = xml2::read_xml(txt)
1
2
3
4
5
xml
#> {xml_document} ...
#> [1] <ows:ServiceIdentification>\n  <ows:Title>GeoServer WFS...
#> [2] <ows:ServiceProvider>\n  <ows:ProviderName>UN-FAO Fishe...
#> ...

可以使用GetFeature请求和特定的typeName从WFS服务中下载数据（如下面的代码块所示）。

可用的名称取决于所访问的Web要素服务。可以使用Web技术以编程方式提取它们，也可以在浏览器中手动滚动浏览GetCapabilities输出的内容。

1
2
3
4
qf = list(request = "GetFeature", typeName = "area:FAO_AREAS")
file = tempfile(fileext = ".gml")
GET(url = base_url, path = endpoint, query = qf, write_disk(file))
fao_areas = read_sf(file)

请注意使用write_disk()来确保结果被写入磁盘而不是加载到内存中，从而可以使用sf导入这些数据。这个示例展示了如何使用httr获得对Web服务的底层访问权限，这对于理解Web服务的工作原理可能很有用。然而，对于许多日常任务来说，更高级的接口可能更合适，为此目的已经开发了许多R包和教程。ows4R包就是专门用于处理OWS服务的。

文件形式

地理数据集通常存储为文件或空间数据库中。文件格式可以存储矢量数据或栅格数据，而空间数据库如PostGIS可以同时存储两种类型的数据。如今，各种文件格式可能令人困惑，但自从20世纪60年代开始使用GIS软件（当时第一个广泛分发的空间分析程序SYMAP在哈佛大学创建）以来，已经进行了许多整合和标准化。

GDAL（应该发音为“goo-dal”，其中双“o”是对面向对象的引用），即地理数据抽象库，自2000年发布以来已经解决了许多与地理文件格式不兼容相关的问题。GDAL提供了一个统一且高性能的接口，用于读取和写入许多栅格和矢量数据格式。[^08-read-write-plot-4] 许多开放和专有的GIS程序，包括GRASS、ArcGIS和QGIS，在其图形用户界面后面使用GDAL来进行摄取和输出地理数据的相关格式。

[^08-read-write-plot-4]：正如我们在几何操作章中提到的，GDAL还包含一组实用函数，可以进行栅格镶嵌、重采样、裁剪、重投影等操作。

GDAL提供对200多种矢量和栅格数据格式的访问。下表呈现了一些常用的空间文件格式的基本信息。

(\#tab:formats)Selected spatial file formats.
Name Extension Info Type Model
ESRI Shapefile .shp (the main file) Popular format consisting of at least three files. No support for: files > 2GB; mixed types; names > 10 chars; cols > 255. Vector Partially open
GeoJSON .geojson Extends the JSON exchange format by including a subset of the simple feature representation; mostly used for storing coordinates in longitude and latitude; it is extended by the TopoJSON format Vector Open
KML .kml XML-based format for spatial visualization, developed for use with Google Earth. Zipped KML file forms the KMZ format. Vector Open
GPX .gpx XML schema created for exchange of GPS data. Vector Open
FlatGeobuf .fgb Single file format allowing for quick reading and writing of vector data. Has streaming capabilities. Vector Open
GeoTIFF .tif/.tiff Popular raster format. A TIFF file containing additional spatial metadata. Raster Open
Arc ASCII .asc Text format where the first six lines represent the raster header, followed by the raster cell values arranged in rows and columns. Raster Open
SQLite/SpatiaLite .sqlite Standalone relational database, SpatiaLite is the spatial extension of SQLite. Vector and raster Open
ESRI FileGDB .gdb Spatial and nonspatial objects created by ArcGIS. Allows: multiple feature classes; topology. Limited support from GDAL. Vector and raster Proprietary
GeoPackage .gpkg Lightweight database container based on SQLite allowing an easy and platform-independent exchange of geodata Vector and (very limited) raster Open

一个重要的发展是为了标准化和开源文件格式，即1994年成立了开放地理空间联盟（OGC）。除了定义简单要素数据模型，OGC还协调开发开放标准，例如在文件格式中使用的KML和GeoPackage。OGC认可的开放文件格式相对于专有格式具有几个优势：标准是公开的，确保透明性，并为用户提供了进一步开发和调整文件格式以适应其特定需求的可能性。

ESRI Shapefile是最流行的矢量数据交换格式；然而，它不是开放格式（尽管其规范是公开的）。它在1990年代初开发，并具有许多限制。首先，它是一个多文件格式，至少由三个文件组成。它只支持255列，列名限制为十个字符，文件大小限制为 2 GB。此外，ESRI Shapefile不支持所有可能的几何类型，例如，它无法区分多边形和多部分多边形。[^08-read-write-plot-5] 尽管存在这些限制，长时间以来都缺乏一个可行的替代方案。与此同时，GeoPackage出现了，似乎是ESRI Shapefile的一个更适合的替代候选。GeoPackage是一种用于交换地理空间信息的格式，也是OGC的标准。GeoPackage标准描述了如何在一个小型SQLite容器中存储地理空间信息的规则。因此，GeoPackage是一个轻量级的空间数据库容器，允许存储矢量和栅格数据，还可以存储非空间数据和扩展。除了GeoPackage，还有其他值得查看的地理空间数据交换格式。

GeoTIFF格式似乎是最显著的栅格数据格式。它允许将空间信息（如CRS）嵌入到TIFF文件中。与ESRI Shapefile类似，这个格式最早于1990年代开发，但作为开放格式。此外，GeoTIFF仍在不断扩展和改进。最近对GeoTIFF格式的最重要的一个新增是其称为COG（云优化的GeoTIFF）的变体。以COG形式保存的栅格对象可以托管在HTTP服务器上，因此其他人可以仅读取文件的部分内容，而无需下载整个文件。

还有许多其他的空间数据格式，由于书的篇幅限制，我们没有详细解释或在表@ref(tab:formats)中提及。此外，新的空间数据格式正在开发中（例如，GeoParquet和Zarr），并且现有的格式也在不断改进。如果您想了解更多关于其他格式的信息，我们建议您阅读GDAL关于矢量和栅格驱动程序的文档。此外，一些空间数据格式可以存储除矢量或栅格外的其他数据模型（类型）。这包括用于存储激光雷达点云的LAS和LAZ格式，以及用于存储多维数组的NetCDF和HDF格式。

空间数据通常也可以使用表格（非空间）文本格式进行存储，包括CSV文件或Excel电子表格。例如，这可以方便地与不使用GIS工具的人共享空间样本，或者与不接受空间数据格式的其他软件交换数据。然而，这种方法可能存在几个问题——对于存储比点更复杂的几何形状来说相当具有挑战性，并且不直接存储有关CRS的信息。

数据输入 (I)

执行诸如sf::read_sf()（我们用于加载矢量数据的主要函数）或terra::rast()（用于加载栅格数据的主要函数）等命令会悄无声息地触发一系列事件，从文件中读取数据。此外，有许多R包包含着各种各样的地理数据，或者提供了对不同数据源的简单访问。它们都会将数据加载到R中，更准确地说，将对象分配给存储在RAM中、从R会话的.GlobalEnv中访问的工作空间中。

矢量数据

空间矢量数据有多种文件格式。最常见的表示形式，如.geojson和.gpkg文件，可以直接使用sf函数read_sf()（或等效的st_read()）导入到R中，它在幕后使用了GDAL的矢量驱动程序。st_drivers()函数返回一个数据框，其中包含前两列的name和long_name，以及每个驱动程序的特性，这些特性对于GDAL（因此也对于sf）可用，包括写入数据的能力和在随后的列中存储栅格数据的能力，正如表@ref(tab:drivers)中对关键文件格式所示。以下命令显示了计算机GDAL安装中报告的前三个驱动程序（结果可能因安装的GDAL版本而异）以及它们特性的摘要。请注意，大多数驱动程序都可以写入数据（87个驱动程序中的51个），而仅有16种格式可以有效地表示栅格数据以及矢量数据（详细信息请参见?st_drivers()）：

1
2
3
sf_drivers = st_drivers()
head(sf_drivers, n = 3)
summary(sf_drivers[-c(1:2)])
(\#tab:drivers)Popular drivers/formats for reading/writing vector data.
name long_name write copy is_raster is_vector vsi
ESRI Shapefile ESRI Shapefile TRUE FALSE FALSE TRUE TRUE
GPX GPX TRUE FALSE FALSE TRUE TRUE
KML Keyhole Markup Language (KML) TRUE FALSE FALSE TRUE TRUE
GeoJSON GeoJSON TRUE FALSE FALSE TRUE TRUE
GPKG GeoPackage TRUE TRUE TRUE TRUE TRUE
FlatGeobuf FlatGeobuf TRUE FALSE FALSE TRUE TRUE

read_sf()的第一个参数是dsn，它应该是一个文本字符串或一个包含单个文本字符串的对象。文本字符串的内容在不同的驱动程序之间可能会有所不同。在大多数情况下，例如ESRI Shapefile（.shp）或GeoPackage格式（.gpkg），dsn将是一个文件名。read_sf()会根据文件扩展名猜测驱动程序，正如下面的示例中展示的针对.gpkg文件：

1
2
f = system.file("shapes/world.gpkg", package = "spData")
world = read_sf(f, quiet = TRUE)

对于某些驱动程序，dsn可以提供为文件夹名称、数据库的访问凭据或GeoJSON字符串表示形式（有关更多详细信息，请参见read_sf()帮助页面的示例）。

一些矢量驱动程序格式可以存储多个数据图层。默认情况下，read_sf()会自动读取在dsn中指定的文件的第一个图层；但是，您可以使用layer参数来指定任何其他图层。

read_sf()函数还允许使用两种可能的机制将文件的部分内容读入RAM。第一个与query参数相关，它允许使用OGR SQL查询文本指定要读取的数据部分。下面的示例从中提取了仅有坦桑尼亚的数据（图@ref(fig:readsfquery)：A）。通过指定我们要从"world"图层中选择所有列（SELECT *），其中name_long等于"Tanzania"：

1
tanzania = read_sf(f, query = 'SELECT * FROM world WHERE name_long = "Tanzania"')

如果您不知道可用列的名称，一个好的方法是仅读取数据的一行，可以使用'SELECT * FROM world WHERE FID = 1'。FID代表要素ID - 大多数情况下，它是行号；但是，其值取决于使用的文件格式。例如，在ESRI Shapefile中，FID从0开始，在其他一些文件格式中从1开始，甚至可以是任意值。

第二种机制使用wkt_filter参数。此参数需要一个表示我们要提取数据的研究区域的well-known text。让我们尝试一个小例子 - 我们希望从文件中读取与坦桑尼亚边界50,000米缓冲区相交的多边形。为了实现这一点，我们需要通过（a）创建缓冲区，（b）使用st_geometry()将sf缓冲区对象转换为sfc几何对象，以及（c）使用st_as_text()将几何对象转化为其well-known text表示形式来准备我们的"过滤器"：

1
2
3
tanzania_buf = st_buffer(tanzania, 50000)
tanzania_buf_geom = st_geometry(tanzania_buf)
tanzania_buf_wkt = st_as_text(tanzania_buf_geom)

现在，我们可以使用wkt_filter参数应用这个"过滤器"：

1
tanzania_neigh = read_sf(f, wkt_filter = tanzania_buf_wkt)

我们的结果如图：B所示，包含了坦桑尼亚以及其50公里缓冲区内的每个国家。


Reading a subset of the vector data using a query (A) and a wkt filter (B).

当然，某些选项是特定于某些驱动程序的^08-read-write-plot-6。例如，考虑坐标存储在电子表格格式（.csv）中的情况。要将此类文件读取为空间对象，我们自然必须指定表示坐标的列名（在下面的示例中为X和Y）。我们可以使用options参数来实现这一点。要了解可能的选项，请参阅相应的GDAL驱动程序描述的’Open Options’部分。对于逗号分隔值（csv）格式，请访问http://www.gdal.org/drv_csv.html。

^08-read-write-plot-6:可以在以下链接找到支持的矢量格式和选项列表： http://gdal.org/ogr_formats.html.

1
2
3
cycle_hire_txt = system.file("misc/cycle_hire_xy.csv", package = "spData")
cycle_hire_xy = read_sf(cycle_hire_txt,
  options = c("X_POSSIBLE_NAMES=X", "Y_POSSIBLE_NAMES=Y"))

除了描述’XY’坐标的列外，单个列还可以包含几何信息。Well-known text (WKT)、well-known binary (WKB)和GeoJSON格式就是这种情况的示例。例如，world_wkt.csv文件有一个名为WKT的列，表示世界各国的多边形。我们将再次使用options参数来指示这一点。

1
2
3
4
5
world_txt = system.file("misc/world_wkt.csv", package = "spData")
world_wkt = read_sf(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT")
# the same as
world_wkt2 = st_read(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT", 
                    quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE)

📌并非所有支持的矢量文件格式都存储有关其坐标参考系统的信息。在这些情况下，可以使用st_set_crs()函数添加缺失的信息。

作为最后一个示例，我们将展示read_sf()如何读取 KML 文件。KML 文件以 XML 格式存储地理信息——这是一种用于创建网页和以应用程序无关的方式传输数据的数据格式。在这里，我们访问来自网络的KML文件。该文件包含多个图层。st_layers()列出了所有可用的图层。我们选择第一个图层Placemarks，并通过read_sf()中的layer参数来指定：

1
2
3
4
5
6
7
8
9
10
11
12
13
u = "https://developers.google.com/kml/documentation/KML_Samples.kml"
download.file(u, "KML_Samples.kml")
st_layers("KML_Samples.kml")
#> Driver: KML 
#> Available layers:
#>              layer_name  geometry_type features fields crs_name
#> 1            Placemarks       3D Point        3      2   WGS 84
#> 2      Highlighted Icon       3D Point        1      2   WGS 84
#> 3                 Paths 3D Line String        6      2   WGS 84
#> 4         Google Campus     3D Polygon        4      2   WGS 84
#> 5      Extruded Polygon     3D Polygon        1      2   WGS 84
#> 6 Absolute and Relative     3D Polygon        4      2   WGS 84
kml = read_sf("KML_Samples.kml", layer = "Placemarks")

到目前为止，本节中所示的所有示例都使用了sf包进行地理数据导入。它快速且灵活，但值得看看其他特定文件格式的包。一个例子是geojsonsf包。一个基准测试表明，相对于读取.geojson文件，它的速度约快于sf包的10倍。

栅格数据

与矢量数据类似，栅格数据也有许多文件格式，其中一些支持多层文件。terra的rast()命令在提供仅有一层的文件时读取单个图层。

1
2
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
single_layer = rast(raster_filepath)

它也适用于您想读取多层文件的情况。

1
2
multilayer_filepath = system.file("raster/landsat.tif", package = "spDataLarge")
multilayer_rast = rast(multilayer_filepath)

之前的所有示例都是从存储在硬盘上的文件中读取空间信息。然而，GDAL 还允许直接从在线资源读取数据，如HTTP/HTTPS/FTP网络资源。我们需要做的唯一的事情就是在文件路径之前添加/vsicurl/前缀。让我们尝试连接到全球2000-2012年间500米分辨率的每月降雪概率数据。12月份的降雪概率以Cloud Optimized GeoTIFF (COG)文件的形式存储在以下网址中：\url{https://zenodo.org/record/5774954/files/clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif}。要读取在线文件，我们只需提供其URL，同时加上/vsicurl/前缀。

1
2
3
4
5
6
7
8
9
10
myurl = "/vsicurl/https://zenodo.org/record/5774954/files/clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif"
snow = rast(myurl)
snow
#> class       : SpatRaster 
#> dimensions  : 35849, 86400, 1  (nrow, ncol, nlyr)
#> resolution  : 0.00417, 0.00417  (x, y)
#> extent      : -180, 180, -62, 87.4  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif 
#> name        : clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0

由于输入数据是COG，实际上我们并没有将这个文件读入到我们的内存中，而是创建了一个连接，没有获取任何值。只有在应用任何基于值的操作（如crop()或extract()）时才会读取其值。这也使我们能够只读取数据的一小部分，而无需下载整个文件。例如，我们可以通过指定其坐标并应用extract()函数来获取雷克雅未克（Reykjavik）12月份的降雪概率（70%）：

1
2
3
4
5
rey = data.frame(lon = -21.94, lat = 64.15)
snow_rey = extract(snow, rey)
snow_rey
#>   ID clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0
#> 1  1                                                         70

通过这种方式，我们只下载了一个值，而不是整个庞大的GeoTIFF文件。

上面的例子只展示了一个简单（但有用）的情况，但还有更多可以探索的内容。/vsicurl/前缀不仅适用于栅格数据，还适用于矢量文件格式。只需在矢量文件的URL前添加该前缀，就可以直接从在线存储中使用read_sf()函数读取矢量。

值得注意的是，/vsicurl/并不是GDAL提供的唯一前缀，还有许多其他前缀，例如/vsizip/用于在不预先解压缩的情况下从ZIP归档中读取空间文件，或者/vsis3/用于即时读取在AWS S3存储桶中可用的文件。您可以在https://gdal.org/user/virtual_file_systems.html上了解更多信息。

数据输出 (O)

写入地理数据可以将数据从一种格式转换为另一种格式，并保存新创建的对象。根据数据类型（矢量或栅格）、对象类别（例如sf或SpatRaster）以及存储信息的类型和数量（例如对象大小、值范围），了解如何以最高效的方式存储空间文件非常重要。接下来的两个章节将演示如何实现这一点。

矢量数据

write_sf()是read_sf()的对应函数。它允许你将sf对象写入各种地理矢量文件格式，包括最常见的格式，如.geojson、.shp 和 .gpkg。根据文件名，write_sf()会自动决定使用哪个驱动程序。写入过程的速度也取决于驱动程序的性能。

1
write_sf(obj = world, dsn = "world.gpkg")

注意：如果你尝试再次写入相同的数据源，该函数会覆盖文件：

1
write_sf(obj = world, dsn = "world.gpkg")

与覆盖文件不同，我们可以使用append = TRUE将新图层添加到文件中，这在几种空间格式中都受到支持，包括GeoPackage。

1
write_sf(obj = world, dsn = "world_many_layers.gpkg", append = TRUE)

另一种选择是使用st_write()，因为它与write_sf()等效。然而，它有不同的默认设置——它不会覆盖文件（在尝试这样做时会返回错误），并显示所写文件格式和对象的简要摘要。

1
2
3
st_write(obj = world, dsn = "world2.gpkg")
#> Writing layer `world2' to data source `world2.gpkg' using driver `GPKG'
#> Writing 177 features with 10 fields and geometry type Multi Polygon.

layer_options参数可以用于许多不同的目的。其中之一是将空间数据写入文本文件。这可以通过在layer_options中指定GEOMETRY来实现。对于简单的点数据集，可以使用AS_XY（它会创建两列新列来存储坐标），或者对于更复杂的空间数据，可以使用AS_WKT（它会创建一个新列，其中包含空间对象的Well-Known Text表示）。

1
2
write_sf(cycle_hire_xy, "cycle_hire_xy.csv", layer_options = "GEOMETRY=AS_XY")
write_sf(world_wkt, "world_wkt.csv", layer_options = "GEOMETRY=AS_WKT")

栅格数据

writeRaster()函数用于将SpatRaster 对象保存到磁盘文件中。该函数需要输入有关输出数据类型和文件格式的信息，同时还接受特定于所选文件格式的GDAL选项（请参阅 ?writeRaster 获取更多详细信息）。

terra包在保存栅格数据时提供了七种数据类型：INT1U、INT2S、INT2U、INT4S、INT4U、FLT4S 和 FLT8S。[^08-read-write-plot-7] 数据类型决定了写入磁盘的栅格对象的位表示。要选择使用哪种数据类型取决于栅格对象的值范围。数据类型能够表示的值越多，文件在磁盘上占用的空间就越大。无符号整数（INT1U、INT2U、INT4U）适用于分类数据，而浮点数（FLT4S 和 FLT8S）通常用于表示连续数据。writeRaster() 默认使用 FLT4S。虽然这在大多数情况下是有效的，但如果要保存二进制或分类数据，输出文件的大小将变得不必要地大。因此，我们建议选择占用存储空间最少的数据类型，但仍能表示所有值（请使用 summary() 函数检查值范围）。

[^08-read-write-plot-7]: 不推荐使用 INT4U，因为 R 不支持32位无符号整数。

(\#tab:datatypes)Data types supported by the terra package.
Data type Minimum value Maximum value
INT1U 0 255
INT2S -32,767 32,767
INT2U 0 65,534
INT4S -2,147,483,647 2,147,483,647
INT4U 0 4,294,967,296
FLT4S -3.4e+38 3.4e+38
FLT8S -1.7e+308 1.7e+308

默认情况下，输出文件的格式由文件名决定。将文件命名为 *.tif 将创建一个 GeoTIFF 文件，如下所示：

1
writeRaster(single_layer, filename = "my_raster.tif", datatype = "INT2U")

一些栅格文件格式具有额外的选项，可以通过将GDAL参数提供给writeRaster() 函数的 options 参数来设置。在 terra 中，默认情况下，GeoTIFF 文件以 LZW 压缩方式写入，参数设置为 gdal = c("COMPRESS=LZW")。要更改或禁用压缩，我们需要修改这个参数。

1
2
writeRaster(x = single_layer, filename = "my_raster.tif",
            gdal = c("COMPRESS=NONE"), overwrite = TRUE)

此外，我们还可以使用filetype="COG"选项将我们的栅格对象保存为 COG（Cloud Optimized GeoTIFF，详见第@ref(file-formats) 节）。

1
2
writeRaster(x = single_layer, filename = "my_raster.tif",
            filetype = "COG", overwrite = TRUE)

输出可视化

R支持许多不同的静态和交互式图形格式。保存静态绘图的最通用方法是打开一个图形设备，创建绘图，然后关闭它，例如：

1
2
3
png(filename = "lifeExp.png", width = 500, height = 350)
plot(world["lifeExp"])
dev.off()

其他可用的图形设备包括 pdf()、bmp()、jpeg() 和 tiff()。您可以指定输出图的多个属性，包括宽度、高度和分辨率。

另外，一些图形包提供了自己的函数来保存图形输出。例如，tmap 包提供了tmap_save()函数。您可以通过指定对象名称和新文件的文件路径，将 tmap 对象保存为不同的图形格式或 HTML 文件。

1
2
3
library(tmap)
tmap_obj = tm_shape(world) + tm_polygons(col = "lifeExp")
tmap_save(tmap_obj, filename = "lifeExp_tmap.png")

另一方面，您可以使用mapview包中的mapshot()函数将创建的交互式地图保存为 HTML 文件或图像：

1
2
3
library(mapview)
mapview_obj = mapview(world, zcol = "lifeExp", legend = TRUE)
mapshot(mapview_obj, file = "my_interactive_map.html")

练习

E1. 列举并描述三种矢量、栅格和地理数据库格式的类型。

E2. 至少命名 sf 函数 read_sf() 和 st_read() 之间的两个区别。

  • read_sf() 是一个更简单、更直接的函数，用于快速读取矢量数据，而 st_read() 提供了更多的选项和控制，包括对图层、驱动程序和其他参数的选择。
  • read_sf() 返回一个简化的数据框（sf对象），而 st_read() 返回一个标准的数据框，并提供了一种选择性地读取几何数据的方法。

E3. 从 spData 包中读取 cycle_hire_xy.csv 文件作为一个空间对象（提示：它位于 misc 文件夹中）。加载对象的几何类型是什么？

E4. 使用 rnaturalearth 下载德国的边界，并创建一个名为 germany_borders 的新对象。将这个新对象写入GeoPackage格式的文件。

E5. 使用 geodata 包下载全球每月最低温度，空间分辨率为五分钟。提取六月的值，并将它们保存到名为 tmin_june.tif 的文件中（提示：使用 terra::subset()）。

E6. 创建德国边界的静态地图，并将其保存为PNG文件。

E7. 使用 cycle_hire_xy.csv 文件中的数据创建一个交互式地图。将此地图导出到名为 cycle_hire.html 的文件。

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