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

以前我一直覺得Python的繪圖工具與R語言ggplot2比起來，不夠優雅，這也是我一直堅定的選擇使用R+ggplot2深入的學習數據可視化的原因，ggplot2在坐標系的整合與兼容性和擴展性上確實技高一籌，所以ggplot2成了可視化的巨無霸，成了可視化界的微信，不僅自身生態日趨完善，而且還有眾多的開發者為其開發輔助功能包（你可以理解為依附於微信的小程序）。

最近偶然在學習Python可視化的過程中，了解到了geopandas，確實第一眼看著很眼熟，或許你第一眼就能把它與pandas聯繫起來。的確，它跟pandas有著千絲萬縷的聯繫，並且繼承了pandas的諸多高頻函數。而geo是什麼鬼呢？

geo是地理信息系統的簡寫，geopandas是Python中用於處理空間地理信息數據的後起之秀（為什麼是後起之秀呢，因為有個叫basemap的包，據說很難用，我還沒有深入了解過呢）。

今天要講解的主角是R語言中的sf包和Python中的geopandas庫。

為什麼今天把geopandas和R語言空間數據可視化寫在一起，因為他們很巧合的用到了相同的地理信息處理技術，無論是數據源的支持上、還是空間數據的結構存儲還是投影設置上都是如此。

我以前寫過大量ggplot2空間數據可視化的文章，但是那些大多是基於shp數據源，和ggplot2中的geom_ploygon或者gemo_map函數來製作的。

對，你沒看錯，真的有這種操作~

shp數據源導入所依賴的maptools包已經快要被遺棄了，而geom_ploygon、gemo_map函數所支持的數據結構複雜難懂，無論是對於初學者還是老手們都是一種挑戰，需要做大量的數據合併、轉換、匹配，在前期的數據處理上花費的時間和代碼量已經遠遠超過了可視化的代碼量。

好在新技術總是不斷地出現，數據源上json格式的數據為我們提供了更為便捷、高效、低廉的空間數據信息，而sf包則可以使用直觀易懂的Simple Features數據結構來從新規整地圖數據源，使得過去需要分別準備地理邊界屬性信息和地理邊界經緯點信息來呈現地理空間信息數據結構。

巧合的是，python中的geopandas用了同樣的 技術來簡化空間數據可視化的複雜度，其核心理念也是通過壓縮單個地理多邊形為一個Simple Features，使得所有的地理多邊形與其屬性信息嚴格對齊，行政一個呈現友好的帶有地理信息數據的數據框。

也許以上描述過於抽象，因為涉及到到的內容比較深入，我實在是不知道該如何把這些內容將的通俗易懂，接下來會使用圖片輔助演示。

R中製作地圖傳統的方式是使用geom_ploygon+maptools+shp數據

library(ggplot2)library(plyr)library(maptools)

#數據導入：

china_map<-readShapePoly("D:/R/rstudy/CHN_adm/bou2_4p.shp")

Warning message:

#已經開始提出警告了！

(替代方案，使用rgdal中的readORG函數或者sf包中的st_read函數)

use rgdal::readOGR or sf::st_read

china_map1<-fortify(china_map)

#從SP（空間數據對象）中剝離地理多邊形邊界點信息和多邊形屬性信息

x<-china_map@data xs<-data.frame(id=row.names(x),x) china_map_data <- join(china_map1, xs, type = "full")

#導入業務數據

province_city <- read.csv("D:/R/rstudy/Province/chinaprovincecity.csv") mydata <- read.csv("D:/R/rstudy/Province/geshengzhibiao.csv")china_data <- join(china_map_data, mydata, type="full")

#可視化代碼：

ggplot(china_data,aes(long,lat))+ geom_polygon(aes(group=group),fill="white",colour="grey60")+ geom_point(data =province_city,aes(x = jd,y = wd),colour="red")+ coord_map("polyconic")+ theme_void()

過程何其辛苦！

為什麼使用maptools+geom_ploygon技術組合這麼辛苦呢，問題出在數據源上，如果你想要詳細了解maptools導入的空間信息結果以及goem_ploygon根據什麼規則映射地圖信息，請看這一篇。

一篇小短文助你打開數據可視化的任督二脈！

我能告訴你的是，geom_ploygon製作地圖的時候，剝離了地理信息邊界點數據和多邊形屬性信息，所以你需要同時兼顧、處理兩個包含空間信息的數據框，如果是對不同區域進行等值線映射，你還需要對這兩個數據框進行合併操作。

而sf包則使用了新的、更為優雅簡潔的空間信息呈現技術——Simple Features

以上便是使用shp+maptools+geom_ploygon技術的核心數據結構概況，接下來我們會跟大家講解新技術組合下所支持的空間數據結構。

sf包則也是同時支持shp數據源和json數據源

library("sf")library("ggplot2")china_map<-st_read("D:/R/rstudy/CHN_adm/bou2_4p.shp",stringsAsFactors=FALSE,quiet=TRUE)Encoding(china_map$NAME)<-"GBK"china_map<-st_read("D:/R/mapdata/State/china.geojson",stringsAsFactors=FALSE,quiet=TRUE)

####轉換編碼

st_crs(china_map)$epsg<-4267st_crs(china_map)$proj4string<-"+proj=longlat +datum=NAD27 +no_defs"china_map<-st_transform(china_map,3395)

#####合併數據

mydata<-data.frame(NAME=unique(china_map$NAME),value=runif(34,1,100))china_map<-merge(china_map,mydata,by="NAME")

ggplot() + geom_sf(data=china_map,aes(fill=value))

所以使用sf提供的新技術，製作數據地圖通常僅需以上幾步。任務量大大縮減。

而Python的geopandas包則也提供了相同的空間數據結構處理技術。

import pandas as pdimport numpy as npimport geopandas as gpimport matplotlib.pyplot as pltprovince_city = pd.read_csv("D:/R/rstudy/Province/chinaprovincecity.csv", encoding = "gb18030") china_map=gp.GeoDataFrame.from_file("D:/R/rstudy/CHN_adm/bou2_4p.shp", encoding = "gb18030")china_map=gp.GeoDataFrame.from_file("D:/R/mapdata/State/china.geojson", encoding = "gb18030")

geopandas包同時支持導入shp素材和json素材，導入之後得數據結構與R語言中得sf導入之後得結構是一致得，地理多邊形邊界點信息都被壓縮成了一個非常整齊的列表存儲，列表內每一個單獨的子項目都代表著一個多邊形。

實際上導入之後，你可以看到它的結構是一種特殊的帶有地理信息列的數據框。

geopandas.geodataframe.GeoDataFrame

這種格式數據框繼承了大多數pandas普通數據框的函數及屬性，可以直接針對其使用plot函數繪圖。

china_map.plot(figsize=(20,12))

其內部數據結構與sf包的數據結構如出一轍。

我們可以給其指定一個數值變數，使得映射出來的地圖各區塊根據數值單獨填色。

china_map.plot(column="AREA",figsize=(20,12),cmap="Greens")

如果你想要在此圖層上添加另外一個散點圖層，則需設定兩個具有同樣投影信息的GeoDataFrame對象。

使用剛才導入的province_city數據，將其合併進china_map中去。

china_map=china_map.merge(province_city,left_on="NAME", right_on="province", how="left")from shapely.geometry import Pointchina_map["center"]=gp.GeoSeries([Point(x, y) for x, y in zip(china_map["jd"], china_map["wd"])])

china_map_ploygon=china_map[["AREA","NAME","geometry","class"]]china_map_ploygon=china_map_ploygonchina_map_point=china_map[["AREA","NAME","center","class"]]china_map_point=china_map_point.set_geometry("center")

###設置相同的投影：

china_map_ploygon.crs={"init": "epsg:3395"}china_map_point.crs ={"init": "epsg:3395"}

###地圖可視化過程：

base=china_map_ploygon.plot(column="AREA", edgecolor="black",figsize=(20,12),cmap="Greens")china_map_point.plot(ax=base,marker="o",color="red",markersize=5)

最後讓我們再次看一下R語言中的sf數據對象和Pyhton中的geodatafame對象的對比。

數據源請移步GitHub：

https://github.com/ljtyduyu/DataWarehouse