注: 本文是R语言sf包的核心开发者和维护者——来自德国明斯特大学的地理信息学教授:
Edzer Pebesma
的一篇关于sf包的简介,发表于2018年7月的R语言期刊,主要讲述了sf的定位、功能、开发现状及现存问题和今后展望,sf包是一个非常了不起的工具,在R语言中引入了空间数量分析领域通用的标准规范(simple feature),结合tidyverse工具箱组合,R语言中处理、转化与绘制地理空间数据的复杂度降了一个数量级。
by Edzer Pebesma
摘要 Simple features是一种在计算机中编码矢量空间数据(点、线、面等)的标准化方法。sf包在R语言中引入了simple features对象,它基本具备和sp、rgeos、rgdal一样的矢量空间数据处理能力。本文主要描述此包的基本功能,其在R语言诸多扩展生态系统中的地位,以及在连接R语言与其他空间计算系统中的潜在价值。
我们可以把“Features”看做是一种包含特定空间位置或范围的“事物”或对象;Featrue gemetry 是指某种空间几何特征(位置或范围),或者你也可以直接把其当成是一个点、点集合、线、线集合、多边形、多边形集合,甚至是以上多种对象的结合。简单来说,simple features 就是线集合、多边形集合的特征(这些线集合或者多边形集合是由很多点连接的直线段构成的)。Features 还包含有其他典型属性(如时间、颜色、名称、质量),我们称之为特征属性。并非所有的空间现象都可以被轻易称为“事物”或“对象”。某些连续性现象:比如水温或者海拔等诸如此类的,我们最好把其看成是来自于连续空间(或时间)上的映射函数(Scheider et al., 2016),这些也经常被作为栅格数据而非矢量数据(点集合、线集合、多边形集合)。
Simple feature 对象 (Herring,2011)是呈现、编码矢量空间数据的国际通用标准,主要用于呈现点、线、多边形等几何对象(IOS,2004)。它被广泛用于各种空间数据库系统(Herring,2010),GeoJSON(Butler et al.,2016),GeoSPL(Perry and Herring,2012),以及开源地理信息工具库:如GDAL(Warmerdam,2008),GEOS(GEOS Development Team,2017) 和liblwgeom(一个PostGIS的附加模块,Obe and Hsu(2015))。
sf 包(Pebesma,2018)是R语言中一个读取、写入、操纵、计算simple features对象的工具包。它几乎重塑了sp包中所有矢量空间数据(点集合、线集合、多边形集合等)处理功能(Pebesma and Bivand,2005;Bivand et al.,2013),rgdal ( Bivand et al.,2017) 和 rgeos ( Bivand and Rundel,2017)。由于sp包具有约400个直接的反向依赖和几千个间接依赖,所以很有必要重写一个包来替代它。
首先,在sp包的开发期间,simple features标准还尚未出现,ESRI shapefile那时在矢量空间数据的存储和转换上来处于统治地位。但是由于ESRI shapefile缺乏清晰开放的标准,其本身混乱、繁多的配置文件及其在呈现空间数据上的诸多缺陷,给sp包造成了不利影响,比如在呈现多边形集合上的孔洞时,盲目的使用封闭外边界来标记孔洞。这种方式严重影响图形绘制,阻碍其与其他同类型工具库之间的兼容性。
simple feature 格式 标准目前已经被广泛采纳,但是sp包仍然习以为常的将矢量空间数据强制转化为R的内部对象。这也意味着你无法灵活的进行双向数据操作。比如:导入数据、操纵数据、导出数据后才能得到同样的几何对象。但对于sf而言,这根本就是不是问题。
另一个重要原因是R语言在读写空间数据(GDAL)以及操纵空间几何对象(GEOS)时重度依赖的外部扩展库均以对simple feature标准给予了强有力的支持。
除此之外,sp和当前比较流行的数据操纵工具箱:如tidyverse(Wickham et., 2017)和 ggplot2(wickham,2016)等之间的兼容性较差。
sf包的主要类型如下:
POINT:一个单点组成的数值型向量 MULTIPOINT:每行由多点组成的数值矩阵 LINESTRING:每行由多点组成的数值矩阵 POLYGON:多个数据矩阵(每行由多点组成)组成的列表(多边形边界内部可能嵌套若干个孔洞) MULTILINRSTRING:多个数值矩阵(每行由多点组成)组成的列表 MULTIPLOYGON:POLYGON结构组成的list GEOMETRYCOLLECTION: 一个或多个以上几何对象组成的结构。
所有的几何对象都具有空值,表示几何对象的缺失(或者NA)。
Categorygory
Functions
binary predicates | st_contains,st_contains_properly,st_convered_by,st_covers,st_crosses,st_disjoint,st_equals,st_equals_exact,st_intersects,st_is_within_distance,st_within,st_touches,st_overlaps |
binary operations | st_relate,st_distance |
unary operations | st_dimension,st_area,st_length,st_is_longlat,st_is_simple,st_is_valid,st_jitter,st_geohash,st_geometry_type,st_sample,st_line_sample,st_join,st_interpolate_aw,st_make_grid,st_graticule,st_extSoftVersion,raw_ToHex,st_proj_info |
setters | st_set_agr,st_set_crs |
constructors | st_sfc,st_sf,st_as_sf,st_as_sfc,st_point,st_multipoint,st_linestring,st_multilinestring,st_polygon,st_multipolygon,st_geometrycollection,st_combine,st_bind_cols |
in & output | st_read,st_read_db,st_write,st_write_db,read_sf,write_sf,st_drivers,st_layers |
plotting | st_viewport,st_wrap_dateline,sf.colors |
表格 1: 来自 sf 包的功能函数 ,按照函数类别顺序
以上所列函数中,部分函数同时作用于属性和几何对象,比如:aggregate 和 summarize 计算分组统计量。
class
methods
sfg | as.matrix, c, coerce, format, head, Ops, plot, print, st_as_binary, st_as_grob, st_as_text, st_transform, st_coordinates, st_geometry, st_boundary, st_buffer, st_centroid, st_convex_hull, st_difference, st_intersection, st_line_merge, st_make_valid, st_node, st_point_on_surface, st_polygonize, st_segmentize, st_simplify, st_split, st_sym_difference, st_triangulate, st_union, st_voronoi, st_cast, st_collection_extract, st_is, st_zm |
sfc | [, [<-, as.data.frame, c, coerce, format, Ops, print, rep, st_as_binary, st_as_text, st_bbox, st_coordinates, st_crs, st_crs<-, st_geometry, st_precision, st_set_precision, str, summary, st_transform, st_boundary, st_buffer,st_centroid, st_convex_hull, st_difference, st_intersection, st_line_merge, st_make_valid, st_node, st_point_on_surface, st_polygonize, st_segmentize, st_simplify, st_split, st_sym_difference, st_triangulate, st_union, st_voronoi, st_cast, st_collection_extract, st_is, st_zm, obj_sum, type_sum |
sf | [, [[<-, $<-, aggregate, cbind, coerce, merge, plot, print, rbind, st_agr, st_agr<-, st_bbox, st_coordinates, st_crs, st_crs<-, st_geometry, st_geometry<-, st_precision, st_set_precision, st_transform, st_boundary, st_buffer, st_centroid, st_convex_hull, st_difference, st_intersection, st_line_merge, st_make_valid, st_node, st_point_on_surface, st_polygonize, st_segmentize, st_simplify, st_split, st_sym_difference, st_triangulate, st_union, st_voronoi, st_cast, st_collection_extract, st_is, st_zm, anti_join, arrange, distinct, filter, full_join, gather, group_by, inner_join, left_join, nest, mutate, rename, right_join, sample_frac, sample_n, select, semi_join, separate, slice, spread, summarise, transmute, ungroup, unite |
crs | $, is.na, Ops, print, st_as_text, st_crs |
表格 2: sf 包中的类方法
st_interpolate_aw函数可以针对几何对象进行基于面积加权的差值计算(Do et al., 2015)。st_join可以基于空间类型连接成对的表格。
sf包的一般方法已经展示在上面表格2中了,其中很多方法主要服务于矢量空间数据的创建、抽取、转换,当然也有很函数属于不经常用到的低频函数。在可能的情况下,方法作用于一个几何对象(sfg)、一个几何对象集合(sfc)或一个带有属性的集合对象集合(sf),同时返回一个相同类的对象。
序列化
simple featrue对象定义了两个序列化标准:常见文本表现形式(WKT)和二进制文本表现形式(WKB)。常见文本表现形式是日常打印时默认的输出格式,sfc列可以利用st_as_sfc函数从WKT格式的字符串向量中直接读取。
> library(sf)
Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3
> (pt <- st_point(c(0,1)))
POINT (0 1)
> (pol <- st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))))
POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))
> st_as_sfc("POINT(0 1)") # returns sfc:
Geometry set for 1 feature
geometry type: POINT
dimension: XY
bbox: xmin: 0 ymin: 1 xmax: 0 ymax: 1
epsg (SRID): NA
proj4string: NA
POINT (0 1)
R native simple feature geometries can be written to WKB usingst_as_binary:
> st_as_binary(st_point(c(0,1)))[1] 01 01 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 f0 3f
> st_as_binary(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))))
[1] 01 03 00 00 00 01 00 00 00 05 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
[26] 00 00 00 00 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 00 00 00 00 00 00 00
[51] 00 f0 3f 00 00 00 00 00 00 f0 3f 00 00 00 00 00 00 00 00 00 00 00 00 00 00
[76] f0 3f 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
同理,二进制编码的几何对象同样可以使用sf_as_sfc函数来进行还原。
在sf包中,所有与底层库GDAL、GEOS和liblwgeom之间的通信,以及与空间数据库之间进行的空间几何对象读写操作,均使用c++编写的二进制序列化和反序列化。这样可以使得代码高效、稳健,对于所有可能的几何对象类型,都使用统一的接口进行操纵。
球面几何对象
GEOS库提供了很多用于处理二维空间的运算函数。对于未做投影处理的地理空间数据,提供的坐标通常是经纬度,表征的是球面上的点,而非投影后的平面。sf包允许针对此类数据进行所有几何操作,但在操作过程中,GEOS包会弹出提示信息。而像st_area
、st_length
、st_distance
、st_is_within_distance
和st_segmentze
之类的专门用于进行球面运算操作的函数,则主要依赖于底层库——lwgeom(Pebesma)。相对于geosphere包(Hijmans, 2016a)而言,这个包的主要优势在于,它支持所有simple feature
对象的距离计算,但geosphere包只能计算点与点之间的距离。st_sample
函数目前的功能已经有所改变,现在主要服务于球面区域点采样后的坐标计算。
如果有一个全面的用于处理球面几何对象的函数系统,那当时再好不过的事情了。目前看来具备这种潜质的候选包主要包括s2(Rubak和Ooms, 2017)(Rubak和Ooms, 2017)、liblwgeom (PostGIS的一部分)、CGAL (Fabri和Pion, 2009)和boost.Geometry。
操纵工具
在sf的开发过程中,为了使新的数据结构与tidyverse能够较好的协同工作,投入了大量的精力。这方面主要体现在给dplyr包提供了很多用于操纵sf对象的方法(见表2)以及给ggplot2提供了用于快速针对sf对象进行绘图映射的geom系列函数。
tidyverse工具箱规定了四个主要原则,如下所示:
1. 返回基础数据结构 我们仅使用R语言中最基础的数据结构,并且完全支持两种序列化标准(WKT、WKB)
2. 使用管道函数来精简代码过程。 这样来设计函数和方法可以非常容易的适用于管道式工作流程。像st_crs
之类的替换函数建议用st_set_crs
这种形式,看起来更加优雅。
3. 拥抱函数式编程。 保持函数类型安全,支持空几何体和空列表,并通过提供缩放和移动多边形选项来创造性地完成了重载操作。
>
pol * 2 + pt
POLYGON ((0 1, 2 1, 2 3, 0 3, 0 1))
st_join之类的空间连接函数,允许用户传递一个与st_intersects兼容的连接函数,从而使应用于连接的空间谓词完全可定制。
4. 人性化设计。 基于过去10年编写和维护sp包的经验,我们一直保持sf包的简介和易学。尽可能抽象通用方法,来保持底层代码占用较少的内存。所有函数和方法均以st_
(对于“spacetime”,遵循PostGIS约定)开头,以保持它们的统一性和可识别性,并可提供tab键来实现函数的模糊提醒。
绘图
图1(左)显示了具有多个属性的“sf”对象的默认图:没有提供颜色参数,默认颜色取决于变量是数值(上)还是因子(下)。图1如下:
图1: 左图:带有两个属性的sf对象的默认图;右图:带有颜色键、坐标轴和经纬度的单个属性的绘图。
图2: 使用ggplot2::geom_sf生成的图,现在弯曲的经纬网遵循固定比例的的经纬度线。
> library(sf)
> nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))`
> plot(nc[, c(9,5)])
When we plot a single attribute, a color key is default (unlesskey.pos=NULL). The following command
> plot(nc[, 9], key.pos = 1, axes = TRUE, graticule = TRUE)
在图一右侧图中增加经纬网和坐标轴。
Figure 2 shows a plot generated byggplot2(version 2.2.1 or later):
> library(ggplot2)
> library(tidyr)
> library(dplyr)
> nc2 <- nc %>% st_transform(32119) %>% select(SID74, SID79, geom) %>%+ gather(VAR, SID, -geom)
> ggplot() + geom_sf(data = nc2, aes(fill = SID)) + facet_wrap( ~ VAR, ncol = 1)
图3:sf依赖的其他R语言工具包及外部扩展系统
对于某些用户来说,开始sf包的学习就如同翻开一本新书,同时合上一本旧书。但是和那本旧书相比,这本新书的内容、页码完全不同。目前还不知道,那些R语言中数百个使用了sp包提供的类和方法的包,是否会、以及何时会将修改为依赖sf包的类和方法。
最常听到的问题是在这本新书中栅格数据在哪里:sp为网格数据提供了简单的类,栅格(Hijmans, 2016b)提供了大量的类和密集的方法来使用它们,并与sp向量类紧密集成。当前版本的栅格数据是通过将sf对象转换为(较小的一组)sp对象,从而使其可以兼容其中的一小部分函数。在撰写本文时,我们只能说,这是一个高度活跃、探索和发展中的领域,我们很乐意向感兴趣的读者指出,这一讨论的中大家关注的主流趋势在向何处发展。除了栅格数据之外,时间序列类的空间特征(例如监测站的观测数据)很难映射成sf对象:要么必须将时间切片放入列中,要么添加一个时间列,并为每个观测重复空间几何特征。栅格数据、空间时间序列和栅格时间序列是该项目今后探索的重点领域。
该软件包的一个新功能是能够检索空间度量,并使用显式度量单元设置距离参数等 (Pebesma et al., 2016):
> st_area(st_transform(nc[1, ], 2264)) # NC state plane, US foot
12244955726 US_survey_foot^2
> st_crs(2264)$units
[1] "us-ft"
> st_area(st_transform(nc[1, ], 2264)) %>% units::set_units(km^2) # convert:
1137.598 km^2
这可能会让人困惑,但却有可能防止一系列的科学错误。
与其他计算系统的连接和可伸缩性
在许多情况下,使用R分析空间数据从导入数据开始,或者从文件或数据库导出数据结束。这种能力主要由可见文本(WKT)和二进制文本(WKB)序列化提供,它们是sf
标准的一部分,并且受到sf
的支持。与GDAL、GEOS和liblwgeom库的通信都使用WKB方法。GDAL目前有93种不同的空间向量数据连接驱动程序(文件格式、数据库、web服务)。图3显示了sf包和其他R包和系统库的依赖关系。之所以将sf包构构筑于这些系统上,主要因为这些系统是由R语言外部致力于空间数据探索的研究机构和社会组织使用和维护的,反映了这些组织在关于空间数据研究上达成的默契和共识。
除了使用GDAL之外,sf还可以直接读写空间数据库。目前主要通过RPostgreSQL来与PostGIS一起工作,当然,使用RPostgres以及DBI来读写空间数据库的功能仍然进一步开发完善中。初步研究表明,使用dbplyr框架可以在R中处理大量耗费内存的空间数据库。这不仅消除了R的内存限制,而且还从这些数据库的持久空间索引中获益。
对于平面数据,sf动态地为空间二进制谓词构建空间索引(st_intersects
,st_containsetc
.) 以及二进制形式 (st_intersection,st_difference 等)。一篇关于在sf中设置空间索引的博文 描述了如何使用索引操纵大内存的空间数据集。对于球面数据,还需要研究liblwgeom或s2提供的索引。
总结与延伸阅读
我们引入了一个新包 —— sf
,在R语言中操纵simple feature对象,并且成为最前沿的用于部分替代sp
包家族的潜力股。它为R语言中空间矢量数据的处理提供了新的基础类,已经得到了广泛的关注和应用。在实现sf过程中,维护了几个经过良好验证的概念(几何对象与属性的分离),为sf创建了新的连接(dplyr、ggplot2、空间数据库),并探讨了新的概念(单位、空间索引等)。
为了进一步了解sf的全部功能及其原理,读者可以参考软件包附带的六篇小品文。
致谢
如果没有Roger Bivand的前期工作和付出,sf
包的面世几乎是不可能的。该包的贡献者包括 Ian Cook, Tim Keitt, Michael Sumner, Robin Lovelace, HadleyWickham, Jeroen Ooms, and Etienne Racin。同时也非常感谢在Githubs上面给sf项目提供问题反馈的所有贡献者。特别鸣谢Dirk Eddelbuettel开发的Rcpp(Eddelbuettel et al., 2011; Eddelbuettel,2013)。
来自R语言联盟的支持对于sf的开发、面世和普及与应用至关重要,我们对此表示感谢,同时匿名审稿人也给我们提供非常宝贵的意见。
参考文献
参考期刊原文
原文链接
https://km.sankuai.com/preview?downloadURL=https%3A%2F%2Fkm.sankuai.com%2Fapi%2Ffile%2Fcdn%2F164977341%2F196731151%3FcontentType%3D0%26isNewContent%3Dfalse&attachmentId=196731151&name=RJ-2018-009.pdf
翻译水平有限,如有错误敬请指正!