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

姚清仿

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

關(guān)鍵詞： GIS; R語言; 經(jīng)濟學(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)濟學(xué)中發(fā)揮著越來越重要的作用。R語言同時具有統(tǒng)計分析、數(shù)據(jù)處理和可視化繪圖三大功能[1]。是一個自由、免費、源代碼開放的軟件，是一種適合推廣應(yīng)用于制圖的優(yōu)秀工具[2]。此外，R語言與其它數(shù)繪圖工具不同，其在具體操作時可輕松地編寫代碼或者調(diào)整其它用戶的代碼，適應(yīng)特定圖形繪制要求[3]。雖然目前已有許多介紹R與GIS繪圖的資料[4-6]，但還沒有關(guān)于R語言就地圖的指北針、比例尺、圖例的連續(xù)或者離散進行整合繪制的資料。此外，R語言版本和R的擴展包間存在兼容問題，導(dǎo)致圖形繪制過程中一些函數(shù)或者功能的無法使用。并且即使是相同功能的R-package也需要不同參數(shù)進行調(diào)用[7]，這成為經(jīng)濟學(xué)研究者的困難之處。基于此，本文以貴州省2017年9個地級行政區(qū)的GDP增長率為例，使用3.5.1版本的R語言、R編輯器RStudio并載入相關(guān)的R-package包，繪制出有指北針、比例尺、圖例的地圖并提供相應(yīng)的腳本與說明，旨在為經(jīng)濟學(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ù)名識別及搜索功能，極大地方便了R編程計算[8]。RStudio安裝程序包可從https：//www.rstudio.com/products/rstudio/download/下載。

1.2 R-package的下載及其使用

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

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

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

library（maptools） #加載maptools包

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

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

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

library（rgdal）#加載rgdal的包。

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

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

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

貴州省2017年9個地級行政區(qū)的GDP增長率的數(shù)據(jù)來自于貴州省統(tǒng)計局的《貴州統(tǒng)計年鑒2018》。通常需要將其于城市的ID對應(yīng)賦值之后，使用函數(shù)read.csv（）將GDP增長率的數(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）#將貴州省地級市的編號整合為數(shù)據(jù)框。由于讀入的地理信息的數(shù)據(jù)是SpatialPolygonsDataFrame類，需要通過@取得其中的數(shù)據(jù)。

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

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

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

在讀取數(shù)據(jù)的基礎(chǔ)之上，需要將數(shù)據(jù)加工處理為繪圖的數(shù)據(jù)，即根據(jù)城市的編號將地理數(shù)據(jù)和GDP增長率的數(shù)據(jù)的整合。首先，將 SpatialPolygonsDataFrame類的地理數(shù)據(jù)轉(zhuǎn)換為數(shù)據(jù)框類型。其次，將地理的信息數(shù)據(jù)和貴州省9個地級行政區(qū)的GDP增長率的數(shù)據(jù)根據(jù)城市編碼合并。再次，通過計算每個城市中心的地理坐標(biāo)，可在每個地級行政區(qū)添加對應(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增長率和數(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）是各個行政區(qū)的編碼，coordinates（guizhou）是計算各個行政區(qū)的地理中心坐標(biāo)。guizhou@data[WT5HZ]$[WT5BZ]NL_NAME_2為各個行政區(qū)名稱，將三者整合為數(shù)據(jù)框。

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

2 圖例結(jié)果說明

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

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

library（maptools）

library（sp）

library（ggplot2）

library（ggsn）

library（viridis）

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

names（GDP）[1]="ID"#將讀入的GDP增長率的數(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進行圖形的繪制，將數(shù)據(jù)可視化。

scale_fill_gradient（low="green"，high="red"，name="GDP增長率（%）"）+#繪制連續(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ù)值。st.size表示的是數(shù)值的字體大小，st.dist表示的是數(shù)值與比例圖標(biāo)的距離，Model表示比例尺圖標(biāo)的顯示類型，執(zhí)行函數(shù)help（scalebar）便可查找相應(yīng)的Model類型。

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

2.2 圖例離散的GIS圖

在經(jīng)濟學(xué)研究中，有時需要將數(shù)據(jù)進行分類展示。為此，需要使用圖例離散函數(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增長率（%）"，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頁）

axis.ticks = element_blank（），

axis.title = element_blank（），

legend.position = "right"

）

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

3 結(jié)束語

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

參考文獻

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

[2] R語言在地球化學(xué)數(shù)據(jù)趨勢面分析中的應(yīng)用 [J]. 數(shù)學(xué)的實踐與認(rèn)識， 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é)報：自然科學(xué)版， 2019， 2019（143）： 124-35.

[8] 王斌會. 多元統(tǒng)計分析及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.