注: 本文是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语言与其他空间计算系统中的潜在价值。
"Simple features" 究竟是什么?
我们可以把“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)等之间的兼容性较差。
- tidyverse 包不仅把操纵对象当做是一个数据框(然而sp 对象则是通过提供方法函数来实现),而且把对象视作一个长度相等的向量组成的列表,这一点儿sp包望尘莫及。
- 在使用ggplot2绘图时,先利用fortify函数将sp对象转化成数据框(该数据框里存放着每一个多边形构成点的信息),以此来尝试“简化”多边形对象,这样既不优雅,也不高效。
基本规范
数据类型
sf包的主要类型如下:
- “sf”: 一个数据框(或者tl_df):包含一到 多个空间几何对象列(通常由一组与数据框等长的列表组成)、一个用于标识当前空间几何对象列(sfc类)的属性(sf_column),
- "sfc": 一个由一组空间几何属性组成的列表列
- "sfg":一个空间几何列表列中的任一个元素(一个几何要素)
- "crs": 一个坐标参考系统(CRS),作为“sfc”对象的性质存储除了“sfg”对象之外,以上所有类型均可看做列表。“sfg”对象可看做是其他类型的子类,这些类主要有以下几种存储格式:
POINT:一个单点组成的数值型向量
MULTIPOINT:每行由多点组成的数值矩阵
LINESTRING:每行由多点组成的数值矩阵
POLYGON:多个数据矩阵(每行由多点组成)组成的列表(多边形边界内部可能嵌套若干个孔洞)
MULTILINRSTRING:多个数值矩阵(每行由多点组成)组成的列表
MULTIPLOYGON:POLYGON结构组成的list
GEOMETRYCOLLECTION: 一个或多个以上几何对象组成的结构。所有的几何对象都具有空值,表示几何对象的缺失(或者NA)。
函数与方法
CategorygoryFunctions
binary predicatesst_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_overlapsbinary operationsst_relate,st_distanceunary operationsst_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_infosettersst_set_agr,st_set_crsconstructorsst_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_colsin & outputst_read,st_read_db,st_write,st_write_db,read_sf,write_sf,st_drivers,st_layersplottingst_viewport,st_wrap_dateline,sf.colors表格 1:来自sf包的功能函数 ,按照函数类别顺序
以上所列函数中,部分函数同时作用于属性和几何对象,比如:aggregate 和 summarize 计算分组统计量。class methods
sfgas.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_zmsfc[, [<-, 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_sumsf[, [[<-, $<-, 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, unitecrs$, 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的开发、面世和普及与应用至关重要,我们对此表示感谢,同时匿名审稿人也给我们提供非常宝贵的意见。