使用R畫地圖資料

skyme發表於2016-02-04
地圖

用R畫地圖資料

首先，從這裡下載中國地圖的GIS資料，這是一個壓縮包，完全解壓後包含三個檔案（bou2_4p.dbf、bou2_4p.shp和bou2_4p.shx），將這三個檔案解壓到同一個目錄下。

image

用R繪製地圖比較簡單。比如畫一下全國範圍的區域，可以用如下程式碼：

library(maptools)
mydat = readShapePoly("china-province-border-data.tar/china/bou2_4p.shp") #地圖包位置，根據自己的角標位置設定
plot(mydat)

生成的圖如下：

image

但是，可以看出這樣繪製的地圖的形狀有些扁平。這是因為，在繪圖的過程中，預設把經度和緯度作為普通資料，均勻平等對待，繪製在笛卡爾座標系上造成的。其實，地球的球面圖形如何對映到平面圖上，在地理學上是有一系列不同的專業演算法的。地圖不應該畫在普通的笛卡爾座標系上，而是要畫在地理學專業的座標系上。

也可以安裝maps和mapdata這兩個包，然後輸入下面的命令：

install.packages(mapdata)
install.packages("mapdata")
library(maps)
library(mapdata)
map("china")

生成的圖如下：

image

其中map()函式還可以加上很多引數，大致如下：

map(database = "world", regions = ".", exact = FALSE, boundary = TRUE,
  interior = TRUE, projection = "", parameters = NULL, orientation = NULL,
  fill = FALSE, col = 1, plot = TRUE, add = FALSE, namesonly = FALSE,
  xlim = NULL, ylim = NULL, wrap = FALSE, resolution = if (plot) 1 else 0,
  type = "l", bg = par("bg"), mar = c(4.1, 4.1, par("mar")[3], 0.1),
  myborder = 0.01, ...)

可以設定資料庫、地區、精確度、邊界等，在這裡就不一一詳述，具體的用法可以?map。

給地圖上色

在實際使用的過程中，我們往往會根據自己的需要對地圖中的某些省份著以特定的顏色，這時就可以通過調節plot命令中的fg引數來予以實現。

首先看看R繪製地圖的原理：

在繪製地圖時，每一個省市自治區或者島嶼都是用一個多邊形來表示的。之前的GIS資料，其實就是提供了每一個行政區其多邊形逐點的座標，然後R軟體通過順次連線這些座標，就繪製出了一個多邊形區域。在上面的資料中，一共包含了925個多邊形的資訊，之所以有這麼多是因為一些省份有很多小的附屬島嶼。在這925個多邊形中，每一個都對應一個唯一的ID，編號分別從1到925。

plot命令中的fg引數在本例中應該是一個長度為925的向量，其第i個分量的取值就代表了地圖中第i個多邊形的顏色。一個簡單的嘗試是執行下面這個命令看看效果：

plot(x,col=gray(924:0/924));

生成的圖如下所示：

image

於是自然就產生了一個問題：如何獲取某一個特定地區的ID，進而設定我們想要的顏色？事實上，在變數x中，就已經儲存了我們想要的資訊。在R中輸入“x[[2]]”或“x$att.data”，會得到一個925行7列的資料框，這其實是bou2_4p.dbf這個檔案中儲存的資訊，之前的read.shape()函式雖然讀取的是bou2_4p.shp檔案，但在預設情況下會把dbf檔案的資訊也放到變數之中。對於這個資料框，其行名就是每一個區域的ID編號，第一列和第二列分別是面積和周長，最後一列是該區域所屬的行政區名，其它的列應該也是一些編號性質的變數。於是，通過查詢相應的行政區對應的行名，就可以對fg引數進行賦值了。下面是我編的一個函式，用來生成所需的col向量：

getColor=function(mapdata,provname,provcol,othercol)
{
     f=function(x,y) ifelse(x %in% y,which(y==x),0);
     colIndex=sapply(mapdata@data$NAME,f,provname);
     col=c(othercol,provcol)[colIndex+1];
     return(col);
 }

其中mapdata是存放地圖資料的變數，在上面的例子中就是x，provname是需要改變顏色的地區的名稱，provcol是對應於provname的代表顏色的向量（名稱和數字均可），othercol是其它地區的顏色。舉例如下：

provname=c("北京市","天津市","上海市","重慶市");
 provcol=c("red","green","yellow","purple");
 plot(x,col=getColor(x,provname,provcol,"white"));

生成的圖資料如下：

image

生成人口資料分佈圖

利用類似的方法就可以根據自己的需要對不同的區域進行著色，下面再舉一例。從國家統計局獲取2007年我國各地區的人口資料，然後根據人口的多少對各省份進行著色。程式如下：

provname=c("北京市","天津市","河北省","山西省","內蒙古自治區",
        "遼寧省","吉林省","黑龍江省","上海市","江蘇省",
        "浙江省","安徽省","福建省","江西省","山東省",
        "河南省","湖北省","湖南省","廣東省",
        "廣西壯族自治區","海南省","重慶市","四川省","貴州省",
        "雲南省","西藏自治區","陝西省","甘肅省","青海省",
        "寧夏回族自治區","新疆維吾爾自治區","臺灣省",
        "香港特別行政區");
pop=c(1633,1115,6943,3393,2405,4298,2730,3824,1858,7625,
        5060,6118,3581,4368,9367,9360,5699,6355,9449,
        4768,845,2816,8127,3762,4514,284,3748,2617,
        552,610,2095,2296,693);
provcol=rgb(red=1-pop/max(pop)/2,green=1-pop/max(pop)/2,blue=0);
plot(x,col=getColor(x,provname,provcol,"white"),xlab="",ylab="");

生成的圖如下：

image

其中顏色越深的地方代表人口數越多，反之為人口數越少。

畫部分省地圖

在繪製地圖的過程中，還有一個比較有用的引數是recs，它是一個由多邊形ID組成的向量，表示在地圖中只畫出這些ID所代表的區域。利用這個引數，就可以畫出某一部分的地圖，例如下面的例子是我國中部六省的地圖：

getID=function(mapdata,provname)
{
    index=mapdata$att.data$NAME %in% provname;
    ids=rownames(mapdata$att.data[index,]);
    return(as.numeric(ids));
}
midchina=c("河南省","山西省","湖北省","安徽省","湖南省","江西省");
plot(x, col = getColor(x, midchina, rep("green", 6),
    "white"), border = "white", xlab = "", ylab = "")

生成的圖：

image

用R畫中國地圖並標註城市位置

畫出的圖上仍然可以用points()函式和text()函式加上點和文字，而maptools包中還提供了一個pointLabel()函式，用來解決文字標籤的重疊問題。

par(mar=rep(0,4))
dat = read.csv(text = "城市,jd,wd
    北 京,116.4666667,39.9
    上 海,121.4833333,31.23333333
    天 津,117.1833333,39.15
    重 慶,106.5333333,29.53333333
    哈爾濱,126.6833333,45.75
    長 春,125.3166667,43.86666667
    沈 陽,123.4,41.83333333
    呼和浩特,111.8,40.81666667
    石家莊,114.4666667,38.03333333
    太 原,112.5666667,37.86666667
    濟 南,117,36.63333333
    鄭 州,113.7,34.8
    西 安,108.9,34.26666667
    蘭 州,103.8166667,36.05
    銀 川,106.2666667,38.33333333
    西 寧,101.75,36.63333333
    烏魯木齊,87.6,43.8
    合 肥,117.3,31.85
    南 京,118.8333333,32.03333333
    杭 州,120.15,30.23333333
    長 沙,113,28.18333333
    南 昌,115.8666667,28.68333333
    武 漢,114.35,30.61666667
    成 都,104.0833333,30.65
    貴 陽,106.7,26.58333333
    福 州,119.3,26.08333333
    臺 北,121.5166667,25.05
    廣 州,113.25,23.13333333
    海 口,110.3333333,20.03333333
    南 寧,108.3333333,22.8
    昆 明,102.6833333,25
    拉 薩,91.16666667,29.66666667
    香 港,114.1666667,22.3
    澳門,113.5,22.2")
library(maps)
library(mapdata)
map("china", col = "darkgray", ylim = c(18, 54), panel.first = grid())
points(dat$jd, dat$wd, pch = 19, col = rgb(0, 0, 0, 0.5))
text(dat$jd, dat$wd, dat[, 1], cex = 0.9, col = rgb(0,
    0, 0, 0.7), pos = c(2, 4, 4, 4, 3, 4, 2, 3, 4, 2, 4, 2, 2,
    4, 3, 2, 1, 3, 1, 1, 2, 3, 2, 2, 1, 2, 4, 3, 1, 2, 2, 4, 4, 2))
axis(1, lwd = 0); axis(2, lwd = 0); axis(3, lwd = 0); axis(4, lwd = 0)

生成的圖：

image

ggplot2繪製地圖

以中國地圖為例，下載最新的ArcGIS向量地圖資料，這種地圖資料包含了很多資訊，這是畫地圖的基礎資料。下載地址：

http://download.csdn.net/detail/lgstarzkhl/9427677

用以下程式碼進行地圖繪製：

library(maptools)
library(ggplot2)
library(plyr)
#讀取地圖檔案
china_map<-readShapePoly("C:/Users/feng/Desktop/chinaprovinceborderdata_tar_gz/china-province-border-data.tar/Lambert/省級行政區.shp")
#提取用於繪圖的地圖資料
x<-china_map@data
xs<-data.frame(x,id.1=seq(0:33)-1)
#將地圖資料轉換為資料框
china_map1<-fortify(china_map)
#新增一個id.1欄位，用於和上面的xs（各省市資料）糅合，合併
china_map1$id.1<-china_map1$id
#去掉china_map1中的id欄位，避免在糅合資料的時候，出現兩個相同欄位id和id.1，保證只用id.1來糅合
china_map2<-china_map1[,-7]
#糅合地圖資料
china_mapdata<-join(china_map2, xs, type = "full")
#繪製地圖
ggplot(china_mapdata, aes(x = long, y = lat, group = group,fill=NAME))+
geom_polygon( )+
geom_path(colour = "grey40")+
scale_fill_manual(values=colours(),guide=FALSE)

生成的圖：

image

總結

使用R的地圖擴充套件包可以畫出各種樣式的地圖，具體的展現形態及結合方式還有待進一步挖掘。

相關文章