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

Python實現(xiàn)帶南海九段線分位數(shù)地圖完整可視化版本（附代碼及數(shù)據(jù)）

1999-2017年中國各省旅游外匯收入分析及可視化（附代碼及數(shù)據(jù)）

數(shù)據(jù)及代碼github地址

1 數(shù)據(jù)下載

蝦神把代碼和數(shù)據(jù)放在了github上，沒接觸過github的人可能對如何下載不太熟悉，這里也簡單介紹下兩種方式。如果想進(jìn)一步了解git這種神器的可以安裝git，然后按第一種方式下載（以下介紹默認(rèn)已安裝git bash）。而第二種方式則完全不用了解git方面的內(nèi)容。

1 git clone

首先新建一個你放數(shù)據(jù)和分析的文件夾。然后先點擊瀏覽蝦神提供的github地址，在頁面中點擊clone or download，即跳出如下的頁面，復(fù)制https的地址。

image

然后右擊打開git bash。敲入命令行。

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

然后敲擊回車，即可發(fā)現(xiàn)開始下載。當(dāng)然由于蝦神這個倉庫內(nèi)容比較多，是一個比較漫長的下載過程。

image

2 直接下載（Download ZIP）

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

image

接下來就進(jìn)入傳統(tǒng)的瀏覽器點擊下載存儲文件的范疇了。后續(xù)就不闡述了。

2 R語言可視化

1 所需包的安裝

由于R語言需要一些包的支持才能實現(xiàn)讀取空間數(shù)據(jù)和可視化。具體這塊前面有博客介紹過。這里不詳述。

R語言讀取空間數(shù)據(jù)以及ArcGIS中OLS工具回歸結(jié)果可視化R語言版

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

2 繪圖

說到繪圖的話，其實R語言中有多套繪圖系統(tǒng)（我稱之為系統(tǒng)是指他們的語法各有特色，使用起來有所差異），比如最基礎(chǔ)的plot系統(tǒng)，lattice系統(tǒng)，聞名遐邇的ggplot2系統(tǒng)。本篇博客最早是打算使用基于lattice的spplot繪圖來進(jìn)行。后面發(fā)現(xiàn)利用spplot繪制這個圖并不簡單，但我也會先介紹使用spplot的畫法。最后我采用的是plot系統(tǒng)來繪制。本次繪圖函數(shù)樣例選擇2008年外匯收入數(shù)據(jù)來進(jìn)行。此外后續(xù)出現(xiàn)的whsr表示外匯收入表格，chinau表示使用的中國省份地圖數(shù)據(jù)（關(guān)聯(lián)有外匯收入數(shù)據(jù)），jdx表示九段線地圖。

由于本次是繪制分位數(shù)地圖，因此在讀取數(shù)據(jù)的基礎(chǔ)上需要先計算分位數(shù)。從蝦神Python代碼里也發(fā)現(xiàn)有計算分位數(shù)進(jìn)行分類的代碼。這里也在R里面寫了一下。這里有兩種方式，一種是采用R包ClassInt來做分類。就如代碼的前兩行，一行代碼解決。此外也可以使用這個包做熟悉的自然間斷點法分類。不過似乎這個分類與蝦神的分位數(shù)地圖結(jié)果有點不一致，所以我也寫了一個自己手動計算的代碼。quantinle是R的分位數(shù)函數(shù)，輸入?yún)?shù)是給定一組向量（數(shù)據(jù)）以及需要的分位（例如25分位數(shù)就是0.25），返回值對應(yīng)分位數(shù)的值。具體的計算步驟見代碼。不詳述。

##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分位數(shù)分界線后，在R里可以通過ifelse給不同省份賦值（或者叫分類，就是按照分位數(shù)范圍將外匯收入劃分為不同類別），這里將分類設(shè)置為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)))))

接著可以設(shè)置色帶。這個顏色設(shè)置看個人喜好，我是隨機挑了一種配色方案，使用colorRamPalette。

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

以上就是主要的準(zhǔn)備工作，接下來就進(jìn)行可視化語句就行了。不過這是基于spplot系統(tǒng)的繪制地圖準(zhǔn)備工作，如果是基于plot還會多一步，后面會說。這里先關(guān)注如何用spplot來繪制地圖。先放代碼和結(jié)果再來解釋。

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

image

等一等，為什么有X。

image

其實我后面畫的時候發(fā)現(xiàn)這個系統(tǒng)畫這個圖感覺亂七八糟很麻煩，連畫九段線也有很多毛病，所以就只是先展示下spplot怎么繪圖，后面有空再來水一篇博——啊呸，說錯了，再來寫一篇博客。前面那個結(jié)果圖是有問題的，僅僅是教學(xué)示意，一定要加上九段線，因為我們中國一點都不能少?。?！

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

這里也展示下ClassInt分類的結(jié)果。當(dāng)然也是教學(xué)用，沒有畫九段線。

image

接下來進(jìn)入plot系統(tǒng)的正經(jīng)實現(xiàn)效果。首先就如蝦神博客效果，我們知道本次實現(xiàn)的可視化本身就是一個拼圖操作，所以先需要用par函數(shù)來分割繪圖位置。par函數(shù)的參數(shù)十分豐富。這次主要用兩個，一個是fig，一個是new。代碼如下。

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

fig的參數(shù)是一個向量，形式是(x1,x2,y1,y2)，這個向量是以畫布左下角為原點，其實表示的就是接下來繪制的圖片占據(jù)圖片位置和大小。按照這里的代碼，接下來繪制的圖是占畫布橫軸的前面30%，縱軸的全部。new是表示允許疊加新的圖，這樣才能把兩個圖畫在一個畫布上。解釋得有點繞口，但是如果你自己嘗試畫個圖就能理解，所以我就不繼續(xù)展開了。

接下來先將plot繪制地圖的關(guān)鍵做一些簡要介紹，這里需要先在plot函數(shù)可視化前，做一個準(zhǔn)備工作。即新建一個字段col，把顏色映射到這個字段上（字段值是顏色的十六進(jìn)制），也是用ifelse的方法。做完之后就可以開始做可視化。

由于蝦神給的九段線數(shù)據(jù)是只有九段線那塊區(qū)域，數(shù)據(jù)經(jīng)緯度范圍與省份數(shù)據(jù)不相同，所以可視化的第一步首先要計算一個能包含這兩個數(shù)據(jù)的數(shù)據(jù)經(jīng)緯度范圍。這里就需要bbox了。

你說的是這個bbox？

image

還是這個bbox？

image

或者是這個？

image

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

image

這個函數(shù)可以把空間數(shù)據(jù)的坐標(biāo)范圍返回給我們。這樣子就可以通過簡單的計算得到包含兩個數(shù)據(jù)的經(jīng)緯度范圍。當(dāng)然我把所有的代碼都連在了一行，所以看起來比較復(fù)雜。這個函數(shù)首先先繪制九段線，同時，繪制的數(shù)據(jù)范圍是包含兩個數(shù)據(jù)的。xlim和ylim是表示繪制的數(shù)據(jù)范圍，需要的是一個二元向量（如c(20,30），即表示繪制經(jīng)（緯）度的范圍是20到30）。main跟spplot一樣是標(biāo)題。axes表示不需要畫任何坐標(biāo)軸。這樣子結(jié)果如圖。

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年中國各省國際旅游外匯收入分位數(shù)專題圖\n數(shù)據(jù)來源:國家統(tǒng)計局',
     axes = F)

image

現(xiàn)在先把省份的數(shù)據(jù)繪制上去就能得到基本的地圖。col是顏色映射的字段，add表示可以把新的圖繪制到原來的圖上，相當(dāng)于PS的疊加一個新圖層。這樣我們就看到了基本的圖形。

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

image

上面的圖很單調(diào)，需要格網(wǎng)做點綴。這個時候就需要另外一個函數(shù)axis了。axis的第一個參數(shù)是side，取值為1（下）,2（左）,3（上）,4（右），意思就是畫上下左右的四條坐標(biāo)軸。at表示要繪制的刻度值范圍，tck是表示刻度值的長度，0的話是不標(biāo)刻度，1的話是網(wǎng)格線，默認(rèn)為-0.01，col是顏色。這里再次用到了bbox，這是為了劃分網(wǎng)格，這里采用的是10經(jīng)度為劃分標(biāo)準(zhǔn)。當(dāng)然我還加上了四條軸，具體的代碼后面會分享。這里就不詳述了，主要闡述函數(shù)用法。結(jié)果如圖，越來越接近蝦神效果圖了。

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')

image

地圖目前相較于蝦神的圖（標(biāo)注不考慮），最關(guān)鍵的差距就是圖例了。這里使用polygon和text來繪制圖例。polygon是繪制多邊形，首先需要給出組成x和y的多邊形坐標(biāo)，border是線的顏色，col是多邊形填充顏色。text是根據(jù)坐標(biāo)標(biāo)注文字。x,y是點坐標(biāo)，字符串是標(biāo)注內(nèi)容。當(dāng)然生成多邊形坐標(biāo)以及標(biāo)注文字坐標(biāo)依舊使用了bbox，這里可以根據(jù)具體地圖替換主要改變前面的系數(shù):1.1還有后面的加數(shù)靈活繪制。效果如圖。

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, "無數(shù)據(jù)")

image

后續(xù)就一步步添加圖例，再加上箱線圖就可以得到蝦神的效果圖了。箱線圖繪制是R語言中比較簡單的，使用boxplot即可。

image

當(dāng)然我們也可以做一版有小窗的，原理同上面一樣，也是設(shè)置par，把九段線部分重繪制一遍。結(jié)果如圖。

image

最后蝦神要求的，箱線圖+點繪圖。這里也描述下。由于只針對單年數(shù)據(jù)所有省份箱線圖，先設(shè)置一個類別字段讓所有省份數(shù)據(jù)都相同，設(shè)置為1，然后以這個字段來繪制boxplot箱線圖。加上點的思路就是以設(shè)置的類別字段作為x，外匯收入為y，顏色與地圖相同。代碼如下。詳細(xì)意義就不闡述了。詳情可參見boxplot的文檔。效果最后如圖。

chinau$boxu <- 1

boxplot(y2008 ~ boxu, data = as.data.frame(chinau))
points(factor(chinau$boxu), chinau$y2008, col = chinau$col, pch = 16)

image

打完收工，至于動圖，由于R語言做動圖會牽扯到另外一個軟件，我們就等下一篇博客再水——呸，再介紹。所有的代碼后面會公開給大家。請稍安勿躁。本文所有使用數(shù)據(jù)來自于蝦神Github，解釋權(quán)由他說了算。最后一句——bbox牛逼。