R語言在經(jīng)濟(jì)學(xué)GIS繪圖中的應(yīng)用

姚清仿

摘要：隨著不同學(xué)科的交叉融合發(fā)展，基于地理信息系統(tǒng)（GIS）的繪圖在經(jīng)濟(jì)學(xué)中發(fā)揮著越來越重要的作用。而具有強(qiáng)大可視化繪圖功能的R語言，通過編寫新代碼或調(diào)整已有代碼可輕松實(shí)現(xiàn)科研中數(shù)據(jù)呈現(xiàn)與圖形繪制的要求。文章以貴州省2017年9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率為例，使用3.5.1版本的R語言、R編輯器RStudio并載入相關(guān)的R-package繪制出有指北針、比例尺、不同類型圖例的地圖并提供相應(yīng)的腳本與說明，為經(jīng)濟(jì)學(xué)研究中直觀地表述研究結(jié)果提供了良好的范例，以方便從事經(jīng)濟(jì)學(xué)研究的人員直接更改使用。

關(guān)鍵詞： GIS; R語言; 經(jīng)濟(jì)學(xué)

【Abstract】?With the development of interdisciplinary integration， mapping based on Geographic Information System （GIS） plays an increasingly important role in economics. R language， which has powerful visual drawing function， can easily meet the requirements of data presentation and graphics drawing in scientific research by writing new codes or adjusting existing codes. Taking the GDP growth rate of nine prefecture-level administrative regions in Guizhou Province in 2017 as an example， this paper draws maps of different types of legends and provides corresponding scripts and explanations by using R language of version 3.5.1， RStudio editor and loading relevant R-package， which provides an intuitive description of the research results in economic research.Thus，people engaged in economic research can implement direct use from this example.

【Key words】 ?GIS; R language; economics

0 引 言

隨著不同學(xué)科的交叉融合發(fā)展，基于地理信息系統(tǒng)（GIS）的繪圖在經(jīng)濟(jì)學(xué)中發(fā)揮著越來越重要的作用。R語言同時(shí)具有統(tǒng)計(jì)分析、數(shù)據(jù)處理和可視化繪圖三大功能[1]。是一個(gè)自由、免費(fèi)、源代碼開放的軟件，是一種適合推廣應(yīng)用于制圖的優(yōu)秀工具[2]。此外，R語言與其它數(shù)繪圖工具不同，其在具體操作時(shí)可輕松地編寫代碼或者調(diào)整其它用戶的代碼，適應(yīng)特定圖形繪制要求[3]。雖然目前已有許多介紹R與GIS繪圖的資料[4-6]，但還沒有關(guān)于R語言就地圖的指北針、比例尺、圖例的連續(xù)或者離散進(jìn)行整合繪制的資料。此外，R語言版本和R的擴(kuò)展包間存在兼容問題，導(dǎo)致圖形繪制過程中一些函數(shù)或者功能的無法使用。并且即使是相同功能的R-package也需要不同參數(shù)進(jìn)行調(diào)用[7]，這成為經(jīng)濟(jì)學(xué)研究者的困難之處?；诖耍疚囊再F州省2017年9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率為例，使用3.5.1版本的R語言、R編輯器RStudio并載入相關(guān)的R-package包，繪制出有指北針、比例尺、圖例的地圖并提供相應(yīng)的腳本與說明，旨在為經(jīng)濟(jì)學(xué)研究中直觀地表述研究結(jié)果提供了良好的范例，具有一定的借鑒意義。

1 材料與方法

1.1 R語言與RStudio軟件的下載

3.5.1版本X64位的R語言安裝程序包從https：//cran.r-project.org/bin/windows/base/old/3.5.1/下載。由于RStudio集成了豐富的開發(fā)界面，提供了 方便的函數(shù)名識(shí)別及搜索功能，極大地方便了R編程計(jì)算[8]。RStudio安裝程序包可從https：//www.rstudio.com/products/rstudio/download/下載。

1.2 R-package的下載及其使用

R語言繪圖功能的實(shí)現(xiàn)很大程度上依賴于不同的R-package包，其中包含了許多用戶可以直接調(diào)用的函數(shù)[9]。本文需要下載的R-package軟件包內(nèi)容見表1。

R語言使用時(shí) R-package只需要安裝一次，但每次使用均需要重新加載[10]。下載和加載R-package包的代碼分別為：

Install，packages（“maptools”）#下載maptools包

library（maptools） #加載maptools包

1.3 數(shù)據(jù)的讀取

省、市、縣的數(shù)據(jù)來自于國(guó)家基礎(chǔ)地理信息中心，也可以來自非官方的全球行政區(qū)域數(shù)據(jù)庫(kù)GADM2。當(dāng)數(shù)據(jù)存儲(chǔ)格式為.shp時(shí)，使用rgdal包中的函數(shù)readOGR（file.choose（））讀取。其中函數(shù)：

file.choose（）可以實(shí)現(xiàn)手動(dòng)點(diǎn)擊數(shù)據(jù)文件，從而將數(shù)據(jù)輸入R語言運(yùn)行的環(huán)境中。其執(zhí)行語句如下：

library（rgdal）#加載rgdal的包。

CHN3=readOGR（file.choose（））#讀取.shp格式的數(shù)據(jù)。

當(dāng)數(shù)據(jù)的存儲(chǔ)格式為.rds時(shí)，使用base包中的函數(shù)readRDS（file.choose（）），base屬于R語言中基礎(chǔ)包，其中的函數(shù)如file.choose（）、readOGR（）無需加載即可調(diào)用。其執(zhí)行語句如下：

CHN3=readOGR（file.choose（））#讀取.rds格式的數(shù)據(jù)。

貴州省2017年9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率的數(shù)據(jù)來自于貴州省統(tǒng)計(jì)局的《貴州統(tǒng)計(jì)年鑒2018》。通常需要將其于城市的ID對(duì)應(yīng)賦值之后，使用函數(shù)read.csv（）將GDP增長(zhǎng)率的數(shù)據(jù)讀入到R語言環(huán)境中。其執(zhí)行語句如下：

guizhou=subset（CHN3，CHN3@data[WT5HZ]$[WT5BZ]NAME_1=="Guizhou"）#獲取貴州省的地理數(shù)據(jù)。

Guizhou=data.frame（row.names（guizhou@data），guizhou@data）#將貴州省地級(jí)市的編號(hào)整合為數(shù)據(jù)框。由于讀入的地理信息的數(shù)據(jù)是SpatialPolygonsDataFrame類，需要通過@取得其中的數(shù)據(jù)。

CityID=Guizhou[，c（"row.names.Guizhou.data."，"NL_NAME_2"）]#獲取城市編號(hào)和名稱對(duì)應(yīng)的數(shù)據(jù)。

read.csv（file.choose（））#讀取存儲(chǔ)格式為.csv并含有城市編號(hào)的GDP 增長(zhǎng)率的數(shù)據(jù)。

1.4 繪圖主要步驟及其代碼實(shí)現(xiàn)

在讀取數(shù)據(jù)的基礎(chǔ)之上，需要將數(shù)據(jù)加工處理為繪圖的數(shù)據(jù)，即根據(jù)城市的編號(hào)將地理數(shù)據(jù)和GDP增長(zhǎng)率的數(shù)據(jù)的整合。首先，將 SpatialPolygonsDataFrame類的地理數(shù)據(jù)轉(zhuǎn)換為數(shù)據(jù)框類型。其次，將地理的信息數(shù)據(jù)和貴州省9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率的數(shù)據(jù)根據(jù)城市編碼合并。再次，通過計(jì)算每個(gè)城市中心的地理坐標(biāo)，可在每個(gè)地級(jí)行政區(qū)添加對(duì)應(yīng)的區(qū)域名稱。最后，調(diào)用函數(shù)ggplot（）即可繪制相應(yīng)的GIS圖。具體代碼如下：

gzmap=fortify（guizhou）#將地理信息的數(shù)據(jù)轉(zhuǎn)換為數(shù)據(jù)框的格式。

GZmap=merge（gzmap，GDP，by.x="id"，by.y="ID"，all.x=TRUE）#將GDP增長(zhǎng)率和數(shù)據(jù)框類型的地理數(shù)據(jù)按照城市編碼ID合并。

pos=data.frame（row.names（coordinates（guizhou）），coordinates（guizhou），guizhou@data[WT5HZ]$[WT5BZ]NL_NAME_2）# row.names（coordinates（guizhou）是各個(gè)行政區(qū)的編碼，coordinates（guizhou）是計(jì)算各個(gè)行政區(qū)的地理中心坐標(biāo)。guizhou@data[WT5HZ]$[WT5BZ]NL_NAME_2為各個(gè)行政區(qū)名稱，將三者整合為數(shù)據(jù)框。

names（pos）=c（"ID"，"Lon"，"Lat"，"cityname"）#將整合數(shù)據(jù)的各列重新命名，分別表示城市編碼、城市經(jīng)度、城市緯度、城市名稱。

2 圖例結(jié)果說明

2.1 圖例連續(xù)的GIS圖

為了更加直觀的展示2017年貴州省9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率，通過執(zhí)行相應(yīng)的R語言代碼繪制了其對(duì)應(yīng)的地理信息圖。具體代碼如下：

library（maptools）

library（sp）

library（ggplot2）

library（ggsn）

library（viridis）

GDP=read.csv（file.choose（））

names（GDP）[1]="ID"#將讀入的GDP增長(zhǎng)率的數(shù)據(jù)的第一列變量命名為ID。

CHN3<-readRDS（file.choose（））

guizhou=subset（CHN3，CHN3@data[WT5HZ]$[WT5BZ]NAME_1=="Guizhou"）

Guizhou=data.frame（row.names（guizhou@data），guizhou@data）

CityID=Guizhou[，c（"row.names.guizhou.data."，"NL_NAME_2"）]

pos=data.frame（row.names（coordinates（guizhou）），coordinates（guizhou），guizhou@data[WT5HZ]$[WT5BZ]NL_NAME_2）

names（pos）=c（"ID"，"Lon"，"Lat"，"cityname"）

gzmap=fortify（guizhou）

GZmap=merge（gzmap，GDP，by.x="id"，by.y="ID"，all.x=TRUE）

ggplot（GZmap，aes（long，lat））+#ggplot進(jìn)行圖形的繪制，將數(shù)據(jù)可視化。

scale_fill_gradient（low="green"，high="red"，name="GDP增長(zhǎng)率（%）"）+#繪制連續(xù)型的圖例，指定其顏色的變化范圍

geom_polygon（aes（group=group，fill=GDPincrease），colour="black"）+#設(shè)置區(qū)域的填充數(shù)據(jù)

coord_equal（）+#指定圖形的橫坐標(biāo)和縱坐標(biāo)的顯示比例

north（GZmap， location = "topleft"， symbol =9）+#繪制指北針

scalebar（GZmap，location="bottomright"，dist =100，height=0.025，dist_unit ="km"，st.size=2.5，st.dist=0.04，transform = TRUE， model = "GRS67"）+#繪制比例尺

geom_text（aes（x=pos[WT5HZ]$[WT5BZ]Lon，y=pos[WT5HZ]$[WT5BZ]Lat，label=cityname），data=pos）+

theme（

panel.grid = element_blank（），#設(shè)置網(wǎng)格線

panel.background = element_blank（），#設(shè)置圖形的背景

axis.text = element_blank（），#設(shè)置坐標(biāo)軸的信息

axis.ticks = element_blank（），#設(shè)置坐標(biāo)軸的刻度

axis.title = element_blank（），#設(shè)置坐標(biāo)軸的標(biāo)題

legend.position = "right"#指定圖例的位置。

）

在以上代碼中，指北針函數(shù)的參數(shù)symbol表示指北針的不同類型，可以使用函數(shù)northSymbols（）查看并選擇相應(yīng)的類型。比例尺的函數(shù)中，dist用于設(shè)置圖上1單位代表的實(shí)際距離的數(shù)值。st.size表示的是數(shù)值的字體大小，st.dist表示的是數(shù)值與比例圖標(biāo)的距離，Model表示比例尺圖標(biāo)的顯示類型，執(zhí)行函數(shù)help（scalebar）便可查找相應(yīng)的Model類型。

執(zhí)行以上代碼，即可得到含有貴州省9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率的GIS圖。具體效果如圖1所示。結(jié)果顯示，黔東南苗族侗族自治州的GDP增長(zhǎng)率小于7%，顯示為白色。其中顏色最深的是黔西南布依苗族自治州，其GDP增長(zhǎng)率組最高約為12.5%。其它的行政區(qū)的顏色越深，GDP增長(zhǎng)率相應(yīng)的就越高。

2.2 圖例離散的GIS圖

在經(jīng)濟(jì)學(xué)研究中，有時(shí)需要將數(shù)據(jù)進(jìn)行分類展示。為此，需要使用圖例離散函數(shù)scale_fill_viridis（），該函數(shù)中的參數(shù)direction表示離散分類的排序方式。Option表示其中顏色的顯示方式，有A-E5種選擇方式。alpha設(shè)置顏色的透明度，begin、end分別設(shè)置顏色的色調(diào)。具體代碼如下：

ggplot（GZmap，aes（long，lat））+

scale_fill_viridis（discrete = TRUE，labels=c（"[5.2，11.4]"，"（11.4，12.1]"，"（12.1，12.5]"） ，name="GDP增長(zhǎng)率（%）"，alpha = 0.4， begin = 0.15， end = 0.3，direction = -1， option = "E"）+#設(shè)置離散圖例。

geom_polygon（aes（group=group，fill=cut_number（GDPincrease， 3）），colour="black"）+

coord_equal（）+

north（GZmap， location = "topleft"， symbol =9）+

scalebar（GZmap，location="bottomright"，dist =100，height=0.025， dist_unit ="km"，st.size=2.5，st.dist=0.04，transform = TRUE， model = "GRS67"）+

geom_text（aes（x=pos[WT5HZ]$[WT5BZ]Lon，y=pos[WT5HZ]$[WT5BZ]Lat，label=cityname），data=pos）+

theme（

panel.grid = element_blank（），

panel.background = element_blank（），

axis.text = element_blank（），

[JY]（下轉(zhuǎn)第353頁(yè)）

axis.ticks = element_blank（），

axis.title = element_blank（），

legend.position = "right"

）

執(zhí)行以上代碼，即可得到含有貴州省9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率圖例離散的GIS圖。研究指示：2017年貴州省9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率被分為了三梯隊(duì)，其中，安順市黔西南布依族苗族自治州的GDP增長(zhǎng)率位于12.1%～12.5%之間屬于第一梯隊(duì)。銅仁地區(qū)、貴陽(yáng)市、六盤水市、黔東南苗族侗族自治州的GDP增長(zhǎng)率位于5.2%～11.4%之間，屬于第三梯隊(duì)。其它地級(jí)行政區(qū)的GDP增長(zhǎng)率則位于11.4%～12.1%之間為第二梯隊(duì)。

3 結(jié)束語

本文基于經(jīng)濟(jì)學(xué)相關(guān)領(lǐng)域的背景，以貴州省2017年9個(gè)地級(jí)行政區(qū)的GDP增長(zhǎng)率為例。使用3.5.1版本的R語言、R編輯器RStudio并載入相關(guān)的R-package包，從數(shù)據(jù)的讀取、繪圖的步驟進(jìn)行闡述。從而繪制出有指北針、比例尺，不同類型圖例的地圖，并提供相應(yīng)的腳本以方便從事經(jīng)濟(jì)學(xué)研究的人員直接更改使用，本文為經(jīng)濟(jì)學(xué)研究中直觀地表述研究結(jié)果提供了良好的范例。今后針對(duì)同類型數(shù)據(jù)分析時(shí)，僅需要更改相應(yīng)的數(shù)據(jù)與文本，即可已完成所有的繪圖工作，適合推廣使用。

參考文獻(xiàn)

[1] ADLER J. R in a nutshell： A desktop quick reference [M]. O'Reilly Media， Inc.， 2010.

[2] R語言在地球化學(xué)數(shù)據(jù)趨勢(shì)面分析中的應(yīng)用 [J]. 數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)， 2019， 49（1）： 210-217.

[3] CARSON M A， BASILIKO N. Approaches to R education in Canadian universities [J]. Fresearch， 2016， 5（2802）.

[4] BIVAND R， PEBESMA E， GóMEZ RUBIO V. Applied Spatial Data Analysis with R [M]. 2013： 59-82.

[5] LOVELACE R， CHESHIRE J. Spatial Data Visualisation with R [M]. 2015： 20.

[6] BRUNSDON C A C， LEX. An Introduction to R for Spatial Analysis and Mapping [M]. Sage， 2015.

[7] R語言在生物科學(xué)研究繪圖中的應(yīng)用 [J]. 華東師范大學(xué)學(xué)報(bào)：自然科學(xué)版， 2019， 2019（143）： 124-35.

[8] 王斌會(huì). 多元統(tǒng)計(jì)分析及R語言建模 [M]. 暨南大學(xué)出版社， 2010.

[9] IHAKA R， GENTLEMAN R. R： A Language for Data Analysis and Graphics [J]. Journal of Computational and Graphical Statistics， 1996， （5）：299-314.

[10]CHAPMAN C， FEIT M D. An Overview of the R Language [M]. Springer International Publishing， 2015.