如何用R繪制地圖

本文主要包含三種繪制地圖的方法：繪制基礎地圖、基于空間數據格式(shapefile)繪制地圖以及如何調用百度地圖和谷歌地圖的數據來繪制地圖。

基礎地圖

方法

從map()包里獲取地圖數據，用geom_polygon()(可以用顏色填充)或者geom_path()(不能填充)繪制。

#install.packages(“ggplot2″)

#install.packages(“maps”)

library(ggplot2)

library(maps) # 為了獲取數據

##

## # ATTENTION: maps v3.0 has an updated ‘world’ map. #

## # Many country borders and names have changed since 1990. #

## # Type ‘?world’ or ‘news(package=”maps”)’. See README_v3. #

# 美國地圖數據

states_map <- map_data(“state”)

head(states_map)

## long lat group order region subregion

## 1 -87.46201 30.38968 1 1 alabama

## 2 -87.48493 30.37249 1 2 alabama

## 3 -87.52503 30.37249 1 3 alabama

## 4 -87.53076 30.33239 1 4 alabama

## 5 -87.57087 30.32665 1 5 alabama

## 6 -87.58806 30.32665 1 6 alabama

# geom_polygon()

ggplot(states_map, aes(x=long,y=lat,group=group)) +

geom_polygon(fill=”white”,colour=”black”) +

labs(title = “USA Map”)

# 中國地圖

library(mapdata)

map(“china”, col = “red4″, ylim = c(18,54), panel.first = grid())

title(“China Map”)

# 世界地圖數據

world_map <- map_data(“world”)

head(world_map)

## long lat group order region subregion

## 1 -69.89912 12.45200 1 1 Aruba

## 2 -69.89571 12.42300 1 2 Aruba

## 3 -69.94219 12.43853 1 3 Aruba

## 4 -70.00415 12.50049 1 4 Aruba

## 5 -70.06612 12.54697 1 5 Aruba

## 6 -70.05088 12.59707 1 6 Aruba

#sort(unique(world_map$region))

# 繪制指定區域的地圖數據

# 繪制歐洲足球五大聯賽所在地

euro <- map_data(“world”, region = c(“UK”,”France”, “Spain”,”Germany”, “Italy”))

ggplot(euro, aes(x=long, y = lat, group=group,fill=region)) +

geom_polygon(colour=”black”) +

scale_fill_brewer(palette = “Set2″) +

scale_y_continuous(limits=c(40,60)) +

scale_x_continuous(limits=c(-25,25)) +

labs(title = ” Euorpe’s Big Five Football Leagues”)

繪制等值區域圖

當我們創建一個地圖后，如果根據變量值對不同區域填充不同的顏色呢?

方法

很簡單，只要把變量值和地圖數據合并在一起，然后把一個變量映射到fill上就可以了。

head(USArrests) # 1973年的數據

## Murder Assault UrbanPop Rape

## Alabama 13.2 236 58 21.2

## Alaska 10.0 263 48 44.5

## Arizona 8.1 294 80 31.0

## Arkansas 8.8 190 50 19.5

## California 9.0 276 91 40.6

## Colorado 7.9 204 78 38.7

crimes <- data.frame(state= tolower(rownames(USArrests)), USArrests)

# 合并數據集

crime_map <- merge(states_map,crimes,by.x=”region”,by.y = “state”)

# head(crime_map)

library(plyr) # 加載數據清洗軟件包

##

## Attaching package: ‘plyr’

##

## The following object is masked from ‘package:maps':

##

## ozone

# 按照 group, order排序

crime_map <- arrange(crime_map,group,order)

# head(crime_map)

ggplot(crime_map, aes(x=long,y=lat, group = group, fill = Assault)) +

geom_polygon(colour = “black”) +

coord_map(“polyconic”) +

labs(title = “USA Map”)

# 更改配色

ggplot(crimes, aes(map_id = state, fill = Assault)) +

geom_map(map = states_map, colour = “black”) +

scale_fill_gradient(low=”#FFFFFF”, high = “#BB4444″) +

expand_limits(x = states_map$long, y = states_map$lat)

對于犯罪率這個指標，從上圖可以看出采用連續取值的方法無法很好地反映出信息，這時采用離散取值反而更容易解釋。

# 離散顏色標度

qa <- quantile(crimes$Assault, c(0,0.2,0.4,0.6,0.8,1.0))

qa

## 0% 20% 40% 60% 80% 100%

## 45.0 98.8 135.0 188.8 254.2 337.0

# 新增一個分位數類別變量

crimes$Assault_q <- cut(crimes$Assault, qa, labels = c(“0-20%”, “20-40%”,”40-60%”,

“60-80%”, “80-100%”),

include.lowest = TRUE)

states <- ddply(states_map, .(region),summarise, lat = mean(lat,na.rm = TRUE),

long = mean(long,na.rm = TRUE))

crimes <- merge(crimes, states, by.x = “state”, by.y = “region”)

# 繪制離散分類地圖

p <- ggplot(crimes, aes(map_id = state, fill = Assault_q)) +

geom_map(map = states_map, colour = “black”) +

scale_fill_brewer(palette = “Set2″) +

expand_limits(x = states_map$long, y =states_map$lat) +

coord_map(“polyconic”) +

labs(fill=”Assault Rate\nPercentile”, title = “USA Map”)

p

# 加入州名對應的標簽

p + geom_text(aes(x=long,y=lat,label=state),size=3,colour=”black”) +

theme_bw() +

xlab(“long”) + ylab(“lat”)

# 如果你想去掉網格線和坐標框，那么接著往下翻!

# 創建空白背景地圖

theme_clean <- function(base_size=12){

require(grid)

theme_grey(base_size)

theme(

axis.title = element_blank(),

axis.text = element_blank(),

panel.background = element_blank(),

panel.grid = element_blank(),

axis.ticks.length = unit(0, “cm”),

axis.ticks.margin = unit(0, “cm”),

panel.margin = unit(0, “lines”),

plot.margin = unit(c(0,0,0,0), “lines”),

complete = TRUE

)

}

p + theme_clean()

## Loading required package: grid

基于空間數據格式(shapefile)創建地圖

ESRI公司的Shapefile文件是描述空間數據的幾何和屬性特征的矢量數據結構的一種格式。 一個Shapefile文件最少包括三個文件：主文件(.shp):存儲地理要素的幾何圖形的文件; 索引文件(.shx):存儲圖形要素與屬性信息索引的文件; dBASE表文件(*.dbf):存儲要素信息屬性的dBase表文件。

除此之外還有可選的文件包括：空間參考文件(.prj), 幾何體的空間索引文件(.sbn 和 .sbx), 只讀的Shapefiles的幾何體的空間索引文件(.fbn 和 .fbx), 列表中活動字段的屬性索引(.ain 和 .aih), 可讀寫Shapefile文件的地理編碼索引(.ixs), 可讀寫Shapefile文件的地理編碼索引(.mxs), dbf文件的屬性索引(.atx), 以XML格式保存元數據(.shp.xml), 用于描述.dbf文件的代碼頁,指明其使用的字符編碼的描述文件(*.cpg)。

需要注意的是，主文件是一個直接存取，變長記錄的文件，其中每個記錄描述一個實體的數據，我們稱之為稱為shape。

下面這個網站中可以下載全球各個國家完整的shapefile格式數據： shapefile數據下載網站

備注：慎用中國地圖數據!

方法

利用maptools()包中的readShapePoly()載入空間數據文件，用fortify()把數據轉化成數據框的格式，然后畫圖。

library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

setwd(“~/Desktop/dataset/map”)

# 載入空間數據并轉化成數據框

china_shp <- readShapePoly(“china/bou2_4p.shp”)

# 返回一個 SpatialPolygonsDataFrame 對象

# str(china_shp)

china_map <- fortify(china_shp)

## Regions defined for each Polygons

# 繪制中國地圖

p <- ggplot(china_map, aes(x = long, y = lat, group = group)) +

geom_path() +

labs(title = “China Map”)

# 繪制空白背景的地圖

p + theme_clean()

這里我們只利用了shapefile中最基本的邊界信息，還沒有對地理信息數據進行更進一步的分析。我們還可以將不同格式的地理數據整合起來，例如如何在上面的地圖上繪制出我國的鐵路、水系分布等內容。

調用百度地圖和谷歌地圖的數據

接下來，我將介紹如何從百度地圖和谷歌地圖中獲取心儀的地圖數據信息。

百度地圖

library(devtools)

#install_github(“badbye/baidumap”)

library(baidumap)

# 隨便輸入幾個經緯度坐標

lon = matrix(c(117.93780, 24.55730, 117.93291, 24.57745, 117.23530, 24.64210,

117.05890, 24.74860), byrow=T, ncol=2)

# 將經緯度坐標轉換成真實地理信息

location = getLocation(lon, formatted = T)

location

## lon=117.9378;lat=24.5573 lon=117.93291;lat=24.57745

## “福建省廈門市海滄區坂南路” “福建省廈門市海滄區大溪路”

## lon=117.2353;lat=24.6421 lon=117.0589;lat=24.7486

## “福建省漳州市南靖縣” “福建省漳州市南靖縣X607″

# 獲取廈門大學經緯度坐標，返回json格式文件

getCoordinate(‘廈門大學’) # json

## 廈門大學

## “{\”status\”:0,\”result\”:{\”location\”:{\”lng\”:118.10229694322,\”lat\”:24.442898974406},\”precise\”:0,\”confidence\”:30,\”level\”:\”商圈\”}}”

ad <- getCoordinate(‘廈門大學’, formatted = TRUE)

names(ad) <- NULL

# 繪制地圖

# 自己修改了一些參數，并將修改后的package掛在github上，所以我選擇從github上安裝ggmap包。

# install_github(“fibears/ggmap”)

library(ggmap)

p <- getBaiduMap(“廈門市思明區”,zoom = 12)

## Map from URL : http://api.map.baidu.com/staticimage?width=400&height=400￠er=118.13453488213,24.468728076403&zoom=12&scale=2

ggmap(p) +

geom_point(aes(x=ad[1], y =ad[2]))

當當當，奇跡發生了——廈門大學位于海中央!這是因為谷歌地圖和百度地圖經緯度坐標存在一定的偏差， 而getBaiduMap()軟件包的作者利用了ggmap()包中坐標轉換的思想，這導致最終繪制出來 的圖形不準確。因此，我們還是考慮轉入谷歌陣營。

谷歌地圖

首先介紹一個概念：火星坐標系統

火星坐標系統是一種國家保密插件，也叫做加密插件或者加偏或者SM模組，其實就是對真實坐標系統進行人為的加偏處理，按照特殊的算法，將真實的坐標加密成虛假的坐標，而這個加偏并不是線性的加偏，所以各地的偏移情況都會有所不同。而加密后的坐標也常被人稱為火星坐標系統。

所有的電子地圖、導航設備，都需要加入國家保密插件。第一步，地圖公司測繪地圖，測繪完成后，送到國家測繪局，將真實坐標的電子地圖，加密成“火星坐標”，這樣的地圖才是可以出版和發布的，然后才可以讓GPS公司處理。第二步，所有的GPS公司，只要需要汽車導航的，需要用到導航電子地圖的，都需要在軟件中加入國家保密算法，將COM口讀出來的真實的坐標信號，加密轉換成國家要求的保密的坐標。這樣，GPS導航儀和導航電子地圖就可以完全匹配，GPS也就可以正常工作了。

由于谷歌地圖被GFW屏蔽了，所以想調用其API需要翻墻。

主要有以下幾個步驟：

安裝SSLedge

在RStudio中更改proxy

首先，我們先簡單認識下http和https：

http是普通超文本協議，其信息室明文傳送，而https就是安全超文本傳輸協議，需要證書和提供安全連接，https是嵌套了SSL加密的http連接，其內容會由SSL先加密，然后再傳送。

為了更方便地使用網絡，我將只使用https代理，對于Http類型的網站使用直接連接的方式。

# 查看信息

Sys.getenv()

# 兩種方式設置proxy

# 利用Sys.setenv()

Sys.setenv(https_proxy=”https://user:password@ip:port”)

# 修改.Renviron 文檔

接下來我們來看看如何調用谷歌地圖的API來繪圖。

# 可以直接在cran中下載package

# install_github(“fibears/ggmap”) # 自己修改了部分鏈接代碼，所以我選擇從github下載

# library(ggmap)

setwd(“~/Desktop/dataset/others”)

# 獲取坐標及地圖數據

ad1 <- as.numeric(geocode(“福建省廈門市思明南路422″,source = “google”))

## Information from URL : https://maps.googleapis.com/maps/api/geocode/json?address=%E7%A6%8F%E5%BB%BA%E7%9C%81%E5%8E%A6%E9%97%A8%E5%B8%82%E6%80%9D%E6%98%8E%E5%8D%97%E8%B7%AF422&sensor=false

xmu <- get_map(“廈門市思明區”,zoom = 13, maptype = “roadmap”)

## Map from URL : https://maps.googleapis.com/maps/api/staticmap?center=%E5%8E%A6%E9%97%A8%E5%B8%82%E6%80%9D%E6%98%8E%E5%8C%BA&zoom=13&size=640×640&scale=2&maptype=roadmap&language=en-EN&sensor=false

## Information from URL : https://maps.googleapis.com/maps/api/geocode/json?address=%E5%8E%A6%E9%97%A8%E5%B8%82%E6%80%9D%E6%98%8E%E5%8C%BA&sensor=false

ggmap(xmu, extent = “normal”) +

geom_point(aes(x=ad1[1], y =ad1[2]))

需要注意的是，利用geocode函數檢索經緯度數據時，最好選擇使用道路數據，這樣可以提高檢索的準確率。

最后引用肖凱大神博客中的一個案例：本例是從地震信息網獲取最近一周的地震數據，得到其經緯度，然后以散點形式繪制在google地圖上，同時也顯示地震發生的密度估計。

# 加載擴展包

# install.packages(“animation”)

# install.packages(“XML”)

library(ggmap)

library(animation)

library(XML)

# 從網頁上抓取數據，并進行清理

webpage <-‘http://data.earthquake.cn/datashare/globeEarthquake_csn.html’

tables <- readHTMLTable(webpage,stringsAsFactors = FALSE)

raw <- tables[[6]]

data <- raw[-1,c(‘V1′,’V3′,’V4′)]

names(data) <- c(‘date’,’lan’,’lon’)

data$lan <- as.numeric(data$lan)

data$lon <- as.numeric(data$lon)

data$date <- as.Date(data$date, “%Y-%m-%d”)

# 用ggmap包從google讀取地圖數據，并將之前的數據標注在地圖上。

ggmap(get_googlemap(center = ‘china’, zoom=4,maptype=’terrain’),extent=’device’) +

geom_point(data=data,aes(x=lon,y=lan),colour = ‘red’,alpha=0.7) +

stat_density2d(aes(x=lon,y=lan,fill=..level..,alpha=..level..),

size=2,bins=4,data=data,geom=’polygon’)+

theme(legend.position = “none”)

## Map from URL : https://maps.googleapis.com/maps/api/staticmap?center=china&zoom=4&size=640×640&scale=2&maptype=terrain&sensor=false

## Information from URL : https://maps.googleapis.com/maps/api/geocode/json?address=china&sensor=false

## Warning: Removed 47 rows containing non-finite values (stat_density2d).

## Warning: Removed 47 rows containing missing values (geom_point).