帶南海九段線分位數地圖可視化（R語言版）

今天帶來一篇承諾蝦神的可視化博客。內容是使用R語言進行帶南海九段線分位數地圖可視化。蝦神的原博文地址以下（Python版）。

Python實現帶南海九段線分位數地圖完整可視化版本（附代碼及數據）

1999-2017年中國各省旅遊外匯收入分析及可視化（附代碼及數據）

數據及代碼github地址

1 數據下載

蝦神把代碼和數據放在了github上，沒接觸過github的人可能對如何下載不太熟悉，這裏也簡單介紹下兩種方式。若是想進一步瞭解git這種神器的能夠安裝git，而後按第一種方式下載（如下介紹默認已安裝git bash）。而第二種方式則徹底不用瞭解git方面的內容。

1 git clone

首先新建一個你放數據和分析的文件夾。而後先點擊瀏覽蝦神提供的github地址，在頁面中點擊clone or download，即跳出以下的頁面，複製https的地址。

而後右擊打開git bash。敲入命令行。

git clone https://github.com/allenlu2008/PythonDemo.git

而後敲擊回車，便可發現開始下載。固然因爲蝦神這個倉庫內容比較多，是一個比較漫長的下載過程。

2 直接下載（Download ZIP）

第二種方式依舊是點擊clone or download。這回是點擊Download ZIP。

接下來就進入傳統的瀏覽器點擊下載存儲文件的範疇了。後續就不闡述了。

2 R語言可視化

1 所需包的安裝

因爲R語言須要一些包的支持才能實現讀取空間數據和可視化。具體這塊前面有博客介紹過。這裏不詳述。

R語言讀取空間數據以及ArcGIS中OLS工具迴歸結果可視化R語言版

這裏主要提一下須要本次可視化須要安裝的幾個包：rgdal，ClassInt以及sp包。

2 繪圖

說到繪圖的話，其實R語言中有多套繪圖系統（我稱之爲系統是指他們的語法各有特點，使用起來有所差別），好比最基礎的plot系統，lattice系統，聞名遐邇的ggplot2系統。本篇博客最先是打算使用基於lattice的spplot繪圖來進行。後面發現利用spplot繪製這個圖並不簡單，但我也會先介紹使用spplot的畫法。最後我採用的是plot系統來繪製。本次繪圖函數樣例選擇2008年外匯收入數據來進行。此外後續出現的whsr表示外匯收入表格，chinau表示使用的中國省份地圖數據（關聯有外匯收入數據），jdx表示九段線地圖。

因爲本次是繪製分位數地圖，所以在讀取數據的基礎上須要先計算分位數。從蝦神Python代碼裏也發現有計算分位數進行分類的代碼。這裏也在R裏面寫了一下。這裏有兩種方式，一種是採用R包ClassInt來作分類。就如代碼的前兩行，一行代碼解決。此外也可使用這個包作熟悉的天然間斷點法分類。不過彷佛這個分類與蝦神的分位數地圖結果有點不一致，因此我也寫了一個本身手動計算的代碼。quantinle是R的分位數函數，輸入參數是給定一組向量（數據）以及須要的分位（例如25分位數就是0.25），返回值對應分位數的值。具體的計算步驟見代碼。不詳述。

##use ClassInt
q5 <- classIntervals(chinau$y2008, 5, style = 'quantile')

##self calculation
q25 <- quantile(whsr$y2008, probs = c(0.25))
q50 <- quantile(whsr$y2008, probs = c(0.5))
q75 <- quantile(whsr$y2008, probs = c(0.75))
iqr <- q75 - q25
u <- q75 + 1.5 * iqr
d <- q25 - 1.5 * iqr
u <- ifelse(u > max(whsr$y2008), max(whsr$y2008), u)
d <- ifelse(d < min(whsr$y2008), min(whsr$y2008), d)
box <- c(d, q25, q50, q75, u)

獲取到box分位數分界線後，在R裏能夠經過ifelse給不一樣省份賦值（或者叫分類，就是按照分位數範圍將外匯收入劃分爲不一樣類別），這裏將分類設置爲plot字段（值爲1到6）。代碼以下。

chinau$plot <- ifelse(chinau$y2008 < box[1], 1, 
                      ifelse(chinau$y2008 < box[2], 2, 
                             ifelse(chinau$y2008 < box[3], 3, 
                                    ifelse(chinau$y2008 < box[4], 4, 
                                           ifelse(chinau$y2008 < box[5], 5, 6)))))

接着能夠設置色帶。這個顏色設置看我的喜愛，我是隨機挑了一種配色方案，使用colorRamPalette。

colspic <- colorRampPalette(c("#f38181", "#fce38a", "#eaffd0", "#95e1d3"))(length(box)+1)

以上就是主要的準備工做，接下來就進行可視化語句就好了。不過這是基於spplot系統的繪製地圖準備工做，若是是基於plot還會多一步，後面會說。這裏先關注如何用spplot來繪製地圖。先放代碼和結果再來解釋。

spplot(chinau, zcol = "plot", scales = list(draw = T), 
       main = "2008年中國各省國際旅遊外匯收入分位數專題圖\n數據來源:國家統計局", col.regions = colspic,
       at = seq(1, 7, 1))

等一等，爲何有X。

其實我後面畫的時候發現這個系統畫這個圖感受亂七八糟很麻煩，連畫九段線也有不少毛病，因此就只是先展現下spplot怎麼繪圖，後面有空再來水一篇博——啊呸，說錯了，再來寫一篇博客。前面那個結果圖是有問題的，僅僅是教學示意，必定要加上九段線，由於咱們中國一點都不能少！！！

先來說spplot的語法，首先第一個參數是繪圖的矢量或柵格數據（這裏是chinau），zcol表示要用來映射顏色的字段（選用了以前的plot字段，即參照分位數的分類），scales = list(draw = T)意思是要把經緯度標出來，若是是F就不標。main是填標題，at是映射字段的繪製數值。其餘參數詳情參見文檔或者等之後我來水——呸，我來介紹。

這裏也展現下ClassInt分類的結果。固然也是教學用，沒有畫九段線。

接下來進入plot系統的正經實現效果。首先就如蝦神博客效果，咱們知道本次實現的可視化自己就是一個拼圖操做，因此先須要用par函數來分割繪圖位置。par函數的參數十分豐富。此次主要用兩個，一個是fig，一個是new。代碼以下。

par(fig = c(0, 0.3, 0, 1), new = T)

fig的參數是一個向量，形式是(x1,x2,y1,y2)，這個向量是以畫布左下角爲原點，其實表示的就是接下來繪製的圖片佔據圖片位置和大小。按照這裏的代碼，接下來繪製的圖是佔畫布橫軸的前面30%，縱軸的所有。new是表示容許疊加新的圖，這樣才能把兩個圖畫在一個畫布上。解釋得有點繞口，可是若是你本身嘗試畫個圖就能理解，因此我就不繼續展開了。

接下來先將plot繪製地圖的關鍵作一些簡要介紹，這裏須要先在plot函數可視化前，作一個準備工做。即新建一個字段col，把顏色映射到這個字段上（字段值是顏色的十六進制），也是用ifelse的方法。作完以後就能夠開始作可視化。

因爲蝦神給的九段線數據是隻有九段線那塊區域，數據經緯度範圍與省份數據不相同，因此可視化的第一步首先要計算一個能包含這兩個數據的數據經緯度範圍。這裏就須要bbox了。

你說的是這個bbox？

仍是這個bbox？

或者是這個？

哦，對不起，串戲了。是這個。

這個函數能夠把空間數據的座標範圍返回給咱們。這樣子就能夠經過簡單的計算獲得包含兩個數據的經緯度範圍。固然我把全部的代碼都連在了一行，因此看起來比較複雜。這個函數首先先繪製九段線，同時，繪製的數據範圍是包含兩個數據的。xlim和ylim是表示繪製的數據範圍，須要的是一個二元向量（如c(20,30），即表示繪製經（緯）度的範圍是20到30）。main跟spplot同樣是標題。axes表示不須要畫任何座標軸。這樣子結果如圖。

plot(jdx, 
     xlim = c(min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])), max(c(bbox(chinau)[1,2], bbox(jdx)[1,2]))),
     ylim = c(min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])), max(c(bbox(chinau)[2,2], bbox(jdx)[2,2]))), 
     main = '2008年中國各省國際旅遊外匯收入分位數專題圖\n數據來源:國家統計局',
     axes = F)

如今先把省份的數據繪製上去就能獲得基本的地圖。col是顏色映射的字段，add表示能夠把新的圖繪製到原來的圖上，至關於PS的疊加一個新圖層。這樣咱們就看到了基本的圖形。

plot(chinau, col = chinau$col, add = T)

上面的圖很單調，須要格網作點綴。這個時候就須要另一個函數axis了。axis的第一個參數是side，取值爲1（下）,2（左）,3（上）,4（右），意思就是畫上下左右的四條座標軸。at表示要繪製的刻度值範圍，tck是表示刻度值的長度，0的話是不標刻度，1的話是網格線，默認爲-0.01，col是顏色。這裏再次用到了bbox，這是爲了劃分網格，這裏採用的是10經度爲劃分標準。固然我還加上了四條軸，具體的代碼後面會分享。這裏就不詳述了，主要闡述函數用法。結果如圖，愈來愈接近蝦神效果圖了。

axis(1, at = seq(round(min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))-10), 
                 round(max(c(bbox(chinau)[1,2], bbox(jdx)[1,2]))+10),
                 10), 
     tck = 1, col = 'grey')

地圖目前相較於蝦神的圖（標註不考慮），最關鍵的差距就是圖例了。這裏使用polygon和text來繪製圖例。polygon是繪製多邊形，首先須要給出組成x和y的多邊形座標，border是線的顏色，col是多邊形填充顏色。text是根據座標標註文字。x,y是點座標，字符串是標註內容。固然生成多邊形座標以及標註文字座標依舊使用了bbox，這裏能夠根據具體地圖替換主要改變前面的係數:1.1還有後面的加數靈活繪製。效果如圖。

polygon(x = c(rep(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))), 7), 
              seq(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))), 
                  round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+6, 
                  1),
              rep(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+6, 7),
              seq(round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+5, 
                  round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1]))), 
                  -1)),
        y = c(seq(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))), 
                  round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])))+3, 
                  0.5), 
              rep(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])))+3, 7),
              seq(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1])))+2.5, 
                  round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))), 
                  -0.5),
              rep(round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))), 7)), 
        border = 'black',
        col = 'white')
text(x = round(1.1*min(c(bbox(chinau)[1,1], bbox(jdx)[1,1])))+10, 
     y = (round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))) + 
                  round(1.1*min(c(bbox(chinau)[2,1], bbox(jdx)[2,1]))) + 3)/2, "無數據")

後續就一步步添加圖例，再加上箱線圖就能夠獲得蝦神的效果圖了。箱線圖繪製是R語言中比較簡單的，使用boxplot便可，這裏不贅述。

固然咱們也能夠作一版，原理同上面同樣，也是設置par，把九段線部分重繪製一遍。結果如圖。

打完收工，至於動圖，因爲R語言作動圖會牽扯到另一個軟件，咱們就等下一篇博客再水——呸，再介紹。全部的代碼後面會公開給你們。請稍安勿躁。本文全部使用數據來自於蝦神Github，解釋權由他說了算。最後一句——bbox牛逼。