空间数据可视化笔记——simple features空间对象基础
是不是感觉被封面图和不明觉厉的题目给骗进来了哈哈哈，今天这篇是理论篇，没有多少案例，而且还很长，所以静不下心的小伙伴儿可以先收藏着，时间充裕了再看。
当今互联网和大数据发展的如此迅猛，大量的运营与业务数据需要通过可视化呈现来给商业分析人员提供有价值的决策信息，而地理信息与空间数据可视化则是可视化分析中至关重要而且门槛较高的一类。
通常除了少数本身具备强大前端开发能力的大厂之外，很多中小型企业在内部预算资源有限的情况下，并不具备自建BI和完整可视化框架的能力。需要借助第三方提供的开源可视化平台或者商用BI可视化工具进行大批量业务数据的地理信息维度数据可视化。
第三方的开源在线可视化平台和商用BI工具意味着要么你需要付出高昂的版权费用（盗版的就不提了，总不能在老外面前丢人吧），要么需要付出高昂的员工培训成本。
最重要的是，无论是在线开源开始化框架（底层的如百度Echarts、阿里的dataV和二次开发的完整可视化产品，如永洪BI、诸葛IO、SmartBI、BDP等等）还是商用BI工具（如Tableau、PowerBI等），在可视化理念设计还是可视化呈现风格以及适用的场景上，均已经做过预设，固定好的，你只能选择接受和遵从。
这就意味着团队分析人员如果没有技术开发人员的配合或者服务商的独家服务支持，可能在软件服务使用期内一直只能接受产品中限定的所有规则，想要自定义自己的风格、扩展性能或者自定义新的可视化框架几无可能。
但是R语言和Python等免费开源软件在近两年随着大数据热潮迅猛发展，其内部生态系统对空间数据可视化的支持越来越完善，而空间数据可视化的前沿应用技术和高性能存储方案随着硬件条件的改善都不再成问题。（当然学习R和Python的过程本身就需要付出高昂的时间成本）R和Python这些特性意味着只要掌握基础语法和核心框架使用技巧，你完全可以避过前端，定制自己的可视化平台。R与Python均提供有可以定制可视化应用的成熟框架。
今天这一篇跟大家分享空间数据可视化应用的前沿基础理念，以R语言为主，最后会贯穿一下Python中的简单实现。因为今天讲到的几何对象地理信息数据结构应用比较广泛，不仅在R和Python中有着重要应用，在PowerBI和Tableau甚至很多开源图表库中都有着很普遍的应用，意义比较重大。
----------------
R
----------------
在R语言中，传统对于地理信息数据的支持主要是通过sp包、maptools包和maps包和ggplot2包中的geom_ploygon()或者geom_line()来实现的。
maptools包用于地理信息数据的导入导出（I/O）,支持较多是shp格式、数据框等。sp用于构造地理信息数据以及进行各种需求的计算和转换等。maps包和ggplot2包用于对地理信息数据按照其自身的投影信息和地理属性进行映射和视觉信息号编码。
但是以上技术组合maptools+sp+ggplot2（maps）面临着很大缺陷，这些地理信息数据结构存储上是分割的，地理信息边界数据和地理信息属性数据是通过列表组合的，且不说将业务数据合并，在实际应用时，需要同时设定ID进行属性信息和地理信息合并，而业务数据的纳入则需要二次合并。这样导致花费在处理与转换地理信息数据上的时间甚至会超过作图的时间。
而今天我要分享的内容就是是空间地理可视化前言应用的新方法，将地理信息数据浓缩成单个列表，每一个单独的地理信息对象都被压缩成数据框中的单个记录，这样无需ID，我们的整个空间地理信息数据框就完美的容纳了属性信息和地理信息，地理信息对象作为一个特殊空间的地理信息字段，其每一个记录都是一个压缩的地理信息几何体，可能是点、线、面，也可能是点族、面族、线族甚至以上对象的混合体。
首先来一段官方对于sf的解释：
What is a feature?
A feature is thought of as a thing, or an object in the real world, such as a building or a tree. As is the case with objects, they often consist of other objects. This is the case with features too: a set of features can form a single feature. A forest stand can be a feature, a forest can be a feature, a city can be a feature. A satellite image pixel can be a feature, a complete image can be a feature too.
Features have a geometry describing where on Earth the feature is located, and they have attributes, which describe other properties. The geometry of a tree can be the delineation of its crown, of its stem, or the point indicating its centre. Other properties may include its height, color, diameter at breast height at a particular date, and so on.
The standard says: “ A simple feature is defined by the OpenGIS Abstract specification to have both spatial and non-spatial attributes. Spatial attributes are geometry valued, and simple features are based on 2D geometry with linear interpolation between vertices. ” We will see soon that the same standard will extend its coverage beyond 2D and beyond linear interpolation. Here, we take simple features as the data structures and operations described in the standard.
英文原文写的太优美，不敢随便翻译，大致意思就是说，sf对象可以通过简单的几个维度（x\y\z\m）来刻画现实世界中自然界几乎所有的可见对象，可能是一座森林、一个城市、卫星云图等。它通常由两部分组成：geometry 和 attributes，前者用于描述对象在自然界的所处的位置（主要由空间地理信息坐标构成），后者用于描述它自身的各维度属性（属性信息可能有很多，取决于我们想要了解该对象的哪方面信息）。
以下是sf对象在R语言中的组织形式：
As attributes are typically stored in data.frame objects (or the very similar tbl_df), we will also store feature geometries in a data.frame column. Since geometries are not single-valued, they are put in a list-column, a list of length equal to the number of records in the data.frame, with each list element holding the simple feature geometry of that feature. The three classes used to represent simple features are:
sf, the table (data.frame) with feature attributes and feature geometries, which contains
sfc, the list-column with the geometries for each feature (record), which is composed of
sfg, the feature geometry of an individual simple feature.
以上是R语言中sf包（也即该项技术在R语言中应用的扩展包）的官方文档所引用的技术资料对于simple features特性的解释。
sf对象所能容纳的控件对象主要有以下几种：
除了这些常用的之外，还有很多极具理论价值但是应用 不多的几何类型：
如果你觉得我说的有点抽象，也可参考下面这张图的示意图解释：
图中绿色线条代表一个Simple feature对象，也即一个地理几何对象的所有属性和信息，蓝色线条代表一个Simple feature geometry（sfg）对象，即单个地理几何对象的地理信息部分（内部主要包含边界点信息，可能有多组，也可能 只有一组，以列表格式存储）。红色线条区域代表所有记录的地理信息属性列，是一个有与数据框等长的列表组成的列，英文表示为Simple feature geometry list-colum（sfc）。
而在传统的maptools导入之后，同样的shp文件数据，属性信息和空间地理信息数据分开的。转换为sf对象之后，整体来看，数据结构呈现更加友好、清晰易读。
以上简单揭示了sf的定义、结构和特点，接下来我们深入到sf的sfc列内部，探索它的基础元素生成过程。
sfc列的每一个单独的元素都是一个地理空间对象集合，可能是单个点、线或者面的集合，也有可能是多个点、线、面的集合。那么这些点、线、面的对象时如何组成的呢，sf包中提供了全套的应用函数和方法来处理sf对象。
sf对象的基本操作和属性方法：
library("sf")
nc <- st_read(system.file("shape/nc.shp", package="sf"))
#nc.shp是sf包内置案例文件，因而加载时无需指定绝对路径。
Reading layer `nc' from data source `D:\R\R-3.4.1\library\sf\shape\nc.shp' using driver `ESRI Shapefile'
Simple feature collection with 100 features and 14 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID): 4267
proj4string: +proj=longlat +datum=NAD27 +no_defs在使用sf包导入shp格式数据的时候，软件输出了该数据集所包含的的一些格式描述信息。主要是行列信息、几何对象类型、维度类型、边界框信息、投影编码信息。
class(nc) #类型
[1] "sf" "data.frame"
attr(nc, "sf_column") #几何对象属性列
"geometry"
methods(class = "sf") #自带方法
[1] $<- [ [[<- aggregate cbind coerce
[7] initialize merge plot print rbind show
[13] slotsFromS3 st_agr st_agr<- st_as_sf st_bbox st_boundary
[19] st_buffer st_cast st_centroid st_convex_hull st_coordinates st_crs
[25] st_crs<- st_difference st_geometry st_geometry<- st_intersection st_is
[31] st_line_merge st_make_valid st_point_on_surface st_polygonize st_precision st_segmentize
[37] st_set_precision st_simplify st_split st_sym_difference st_transform st_triangulate
[43] st_union st_voronoi st_zm
see '?methods' for accessing help and source code
nc.no_sf <- as.data.frame(nc)
class(nc.no_sf)
## [1] "data.frame"
虽然带有sf属性的数据框可以强制转化为普通数据库，但是转换之后将会丢失掉所有的sf自带属性。
However, such objects:
no longer register which column is the geometry list-column
no longer have a plot method, and
lack all of the other dedicated methods listed above for class sf
nc_geom <- st_geometry(nc) #提取几何对象列
sf对象的基本元素操作：
point:
x <- st_point(c(1,2));x;class(x);str(x);plot(x)
POINT (1 2)
[1] "XY" "POINT" "sfg"
Classes 'XY', 'POINT', 'sfg' num [1:2] 1 2
plot(x)
multipoint:
p <- rbind(c(3.2,4), c(3,4.6), c(3.8,4.4), c(3.5,3.8), c(3.4,3.6), c(3.9,4.5))
mp<-st_multipoint(p);mp;class(mp);str(mp);plot(mp)
linestring:
s1 <- rbind(c(0,3),c(0,4),c(1,5),c(2,5))
ls <- st_linestring(s1);ls;class(ls);str(ls);plot(ls)
multilinestring:
s2 <- rbind(c(0.2,3), c(0.2,4), c(1,4.8), c(2,4.8))
s3 <- rbind(c(0,4.4), c(0.6,5))
mls <- st_multilinestring(list(s1,s2,s3));mls;class(mls);str(mls);plot(mls)
polygon:
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2));pol;class(pol);str(pol);plot(pol)
multipolygon:
p1 <- rbind(c(0,0), c(1,0), c(3,2), c(2,4), c(1,4), c(0,0))
p2 <- rbind(c(1,1), c(1,2), c(2,2), c(1,1))
pol <-st_polygon(list(p1,p2))
p3 <- rbind(c(3,0), c(4,0), c(4,1), c(3,1), c(3,0))
p4 <- rbind(c(3.3,0.3), c(3.8,0.3), c(3.8,0.8), c(3.3,0.8), c(3.3,0.3))[5:1,]
p5 <- rbind(c(3,3), c(4,2), c(4,3), c(3,3))
mpol <- st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5)));mpol;class(mpol);str(mpol);plot(mpol)
geometrycollection:
gc <- st_geometrycollection(list(mp, mpol, ls))
gc;class(gc);str(gc);plot(gc)
以上是最为常见的七种空间几何对象列举，分别对应点、线、面、多点、多线、多面、点线面集合。这些空间几何对象都可以封装在一个单独的list中，同时与地理信息属性或者其他任何物理空间事物进行匹配组成一个sf对象，进而更完美的呈现空间对象特征。
编码方式：WKT/WKB
x <- st_linestring(matrix(10:1,5))
st_as_text(x)
st_as_binary(x)
[1] 01 02 00 00 00 05 00 00 00 00 00 00 00 00 00 24 40 00 00 00 00 00 00 14 40 00 00 00 00 00 00 22 40 00 00 00 00 00 00 10 40
[42] 00 00 00 00 00 00 20 40 00 00 00 00 00 00 08 40 00 00 00 00 00 00 1c 40 00 00 00 00 00 00 00 40 00 00 00 00 00 00 18 40 00
[83] 00 00 00 00 00 f0 3f因为计算机内存和硬盘在信息存储机制上的差异，在存储资源成本高昂的情况下，对于大批量数据而言（特别是空间地理信息数据），需要进行二进制转码。
WKB/WKB与原生R对象之间的转换：
st_as_sfc("LINESTRING(10 5, 9 4, 8 3, 7 2, 6 1)")[[1]]
LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)
st_as_sfc(structure(list(st_as_binary(x)), class = "WKB"))[[1]]
LINESTRING (10 5, 9 4, 8 3, 7 2, 6 1)Reading and writing（sf对象读写）
setwd("E:/数据可视化/R/R语言学习笔记/数据可视化/ggplot2/地图可视化/sf")
st_write(nc,"nc.shp")
write_sf(nc,"nc.shp") #静默写入：
nc_data<-st_read("nc.shp",stringsAsFactors=FALSE,quiet=TRUE)
nc_data<-read_sf("nc.shp",stringsAsFactors=FALSE) #静默写出以上代码演示了sf对象可以由shp格式地理信息文件读取，同时也可写出为shp格式的地理信息文件。同时需要提醒大家的是，sf对象除了支持读入shp文件之外，也支持json格式地理信息文件，这里又多了一条获取空间地理信息数据的途径。
Coordinate reference systems and transformations（坐标系统转换）
st_crs(nc)<-"+proj=longlat +datum=NAD27 +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat"
#读入数据没有任何crs属性时，需要直接设置crs属性，可以赋给一个projstring字符串，也可以设置WSG84空间投影代码。
nc.web_mercator <- st_transform(nc, 3857) 当自带的投影CRS信息需要转换市，使用以上函数进行转换。
Conversion, including to and from sp（sf对象和sp对象之间的转换）
nc.sp <- as(nc, "Spatial") #sf to sp
nc2 <- st_as_sf(nc.sp) #sp to sfsf对象就是以上提到的Simple features空间地理信息数据框，而sp对象就是我们使用传统的sp包和maptools包导入的 spatial PolygonsDataFrame。
这些又给我们提供了一个很好的导入传统sp对象的方法，你可以选择先将一个空间地理信息文件导入成sf对象，然后再转换成sp对象，之后提供给plot函数、ploygon函数或者ggplot函数使用。仅需一个简单的as函数即可。这一点很有必要说明，因为你现在继续运行老代码的时候，特别是使用maptools导入shp格式数据，已经开始出现警示，并且忠告我们readShapePoly函数即将被遗弃，并强烈建议使用rgdal包中的readORG函数或者sf包中st_read函数。
如果你还是继续停留在复制黏贴代码，不想跟进新技术的话，那么过不了多久，你的老代码都要失效了。
对于sf对象的应用，R语言系统中的plot系统、grid系统和ggplot2系统都提供原生的支持，特别是ggplot2的开发版（开发版板一般都托管在GitHub上，正式版里面暂时还没有提供sf的接口）已经提供了了sf的接口，看官方的最新文档你会看到 多了一个geom_sf()几何对象函数，这就意味着ggplot2为这项sf新技术单独写了一套优化方案，今后的空间数据可视化再也不会是geom_ploygon()一家独大的天下了，势必会出现geom_ploygon和geom_sf双雄争霸，而且sf代表着新技术，很有可能后来居上，让我们拭目以待吧，我以后的案例也都会倾向于提供geom_sf的应用案例。
如果你想要了解更多geom_sf()的用法，请参见这两篇文章，同时可以直接去官方源文档查看介绍。
左手用R右手Python系列12——空间数据可视化与数据地图
-----------------------
Python:
-----------------------
对于sf对象的理论和属性信息的阐述，下面就不再花费时间了，这里我只解释sf技术在Python中的实际应用。Python中的空间地理信息数据可视化主要依赖geopandas，关于这一点，前一篇文章已经有过介绍了，geopandas中主要有两种数据对象，GeoDataFrame和GeoSeries,其中GeoSeries列便是存储着空间地理信息数据的列表集合对象（geometry），其理念与R中的sf对象是一致的。
左手用R右手Python系列12——空间数据可视化与数据地图
china_map=gp.GeoDataFrame.from_file("D:/***/bou2_4p.shp", encoding = 'gb18030')
china_map=gp.GeoDataFrame.from_file("D:/***/china.geojson", encoding = 'gb18030')Python中的GeoDataFrame也同时提供了对于shp、json格式空间地理信息数据源的支持。
在Python中，基础的点线面几何对象主要是通过shapely包来进行支持的。它提供了如同R语言中的sf一样的地理信息数据格式，先将独立几何对象的空间信息进行压缩封装在一个独立的空间几何对象中，然后用这些独立空间几何对象组成空间几何对象集，也即一列由列表组构成的GeoSeries，同时也可为这些独立对象配备属性值信息，最终形成的GeoDataFrame，就是和R语言中的sf(simple features)对象一致的，含有地理空间信息集合的数据框。
按照以上对于R语言部分的描述，同样只取7种最为常用的空间几何对象进行讲解。
- Points / Multi-Points
- Lines / Multi-Lines
- Polygons / Multi-Polygons
#导入shapely库中所有以上对象的函数：
from shapely.geometry import Point
from shapely.geometry import MultiPoint
from shapely.geometry import LineString
from shapely.geometry import MultiLineString
from shapely.geometry.polygon import LinearRing
from shapely.geometry import Polygon
from shapely.geometry import MultiPolygon#########Point######
#Point
point = Point(0.0, 0.0)
Point(point)点对象的一些属性信息：
point.area
point.length
point.bounds
point.x
point.y
point.geom_type
point.distance
point.coords[:]#Collections of Points
points = MultiPoint([(0.0, 0.0), (1.0, 1.0)])
MultiPoint(points)
#点集合的属性：
points.area
points.length
points.bounds
###########Line###########
#LineStrings
line = LineString([(0, 0), (1, 1)])
LineString(line)
line.area
line.length
line.bounds
line.coords
list(line.coords)
line.coords[:]
#线的各种属性信息#Collections of Lines
coords = [((0, 0), (1, 1)), ((-1, 0), (1, 0))]
lines = MultiLineString(coords)
MultiLineString(lines)
MultiLineString(lines.geoms)
lines.area
lines.length
lines.bounds
len(lines.geoms)
list(lines.geoms)
############Polygon############
#Polygons
polygon = Polygon([(0, 0), (1, 1), (1, 0)])
Polygon(polygon)
polygon.area
polygon.length
polygon.bounds
list(polygon.exterior.coords)
list(polygon.interiors)#Collections of Polygons
coords = [(0, 0), (1, 1), (1, 0)]
r = LinearRing(coords)
s = Polygon(r)
t = Polygon(s.buffer(1.0).exterior,[r])
polygons = MultiPolygon([polygon,s,t])
MultiPolygon(polygons)
len(polygons.geoms)
polygons.bounds
len(polygons.geoms)
len(polygons)
以上演示了在Python中构建基础点、线、面以及点集合、线集合、面集合的构造方法。不要觉得这些东西很简单，复杂的空间数据模型和框架都是由这些不起眼的点线面集合构成的，接下来可以查看下在geopandas导入shp或者json后，集合对象列的格式：
import geopandas as gp
%matplotlib inline