空間數據可視化筆記——simple features空間對象基礎

01-24

是不是感覺被封面圖和不明覺厲的題目給騙進來了哈哈哈，今天這篇是理論篇，沒有多少案例，而且還很長，所以靜不下心的小夥伴兒可以先收藏著，時間充裕了再看。

當今互聯網和大數據發展的如此迅猛，大量的運營與業務數據需要通過可視化呈現來給商業分析人員提供有價值的決策信息，而地理信息與空間數據可視化則是可視化分析中至關重要而且門檻較高的一類。

通常除了少數本身具備強大前端開發能力的大廠之外，很多中小型企業在內部預算資源有限的情況下，並不具備自建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語言中，傳統對於地理信息數據的支持主要是通過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對象可以通過簡單的幾個維度（xyzm）來刻畫現實世界中自然界幾乎所有的可見對象，可能是一座森林、一個城市、衛星雲圖等。它通常由兩部分組成：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:RR-3.4.1librarysfshape c.shp" using driver `ESRI Shapefile"Simple feature collection with 100 features and 14 fieldsgeometry type: MULTIPOLYGONdimension: XYbbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965epsg (SRID): 4267proj4string: +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-columnno longer have a plot method, andlack all of the other dedicated methods listed above for class sfnc_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 2plot(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 spnc2 <- st_as_sf(nc.sp) #sp to sf

sf對象就是以上提到的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語言可視化——關於ggplot所支持的數據地圖素材類型

左手用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 Pointfrom shapely.geometry import MultiPointfrom shapely.geometry import LineStringfrom shapely.geometry import MultiLineStringfrom shapely.geometry.polygon import LinearRingfrom shapely.geometry import Polygonfrom shapely.geometry import MultiPolygon

#########Point######

#Point

point = Point(0.0, 0.0)Point(point)

點對象的一些屬性信息：

point.areapoint.lengthpoint.boundspoint.xpoint.ypoint.geom_typepoint.distancepoint.coords[:]

#Collections of Points

points = MultiPoint([(0.0, 0.0), (1.0, 1.0)])MultiPoint(points)#點集合的屬性：points.areapoints.lengthpoints.bounds

###########Line###########

#LineStrings

line = LineString([(0, 0), (1, 1)])LineString(line)line.arealine.lengthline.boundsline.coordslist(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.arealines.lengthlines.boundslen(lines.geoms)list(lines.geoms)

############Polygon############

#Polygons

polygon = Polygon([(0, 0), (1, 1), (1, 0)])Polygon(polygon)polygon.areapolygon.lengthpolygon.boundslist(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.boundslen(polygons.geoms)len(polygons)

以上演示了在Python中構建基礎點、線、面以及點集合、線集合、面集合的構造方法。不要覺得這些東西很簡單，複雜的空間數據模型和框架都是由這些不起眼的點線面集合構成的，接下來可以查看下在geopandas導入shp或者json後，集合對象列的格式：

import geopandas as gp%matplotlib inlinechina_map=gp.GeoDataFrame.from_file("D:/R/mapdata/State/china.geojson",encoding="gb18030")

通過geojson或者shp導入的數據，其地理信息列（geometry）就是一組Polygon或者MutiPolygon的集合，這也是我們使用最多的地理信息空間數據格式，如果你想知道更多的關於R或者Python中地理信息可視化呈現的信息，可以參考這兩篇文章。

R語言可視化——關於ggplot所支持的數據地圖素材類型

左手用R右手Python系列12——空間數據可視化與數據地圖

以下幾個網址分別是R語言中的sf包主頁，sf項目簡介，Python中的geopandas主頁和python中shapely包對於空間幾何體構造的描述，相信這些內容對於你理解sf技術在空間地理信息方面的應用會有很大價值。

R——sf:

https://r-spatial.github.io/sf/https://r-spatial.github.io/sf/articles/sf1.htmlhttps://r-spatial.github.io/sf/articles/sf2.htmlhttps://r-spatial.github.io/sf/articles/sf3.htmlhttp://ggplot2.tidyverse.org/reference/ggsf.html

Python——geopandas:

http://geopandas.org/index.htmlhttp://toblerity.org/shapely/manual.html#linestrings

遺憾的是，以上所有資料都是全英文的，當然了想要學習新技術就要做好了跟英文鬥智斗勇的準備。

在線課程請點擊文末閱讀原文：

Hellobi Live | 9月12日 R語言可視化在商務場景中的應用

數據源文件請移步本人GitHub：

https://github.com/ljtyduyu/DataWarehouse/tree/master/File