数据科学中的R和Python: 基于OpenStreetMap的地理信息绘图

星期六, 四月 14, 2012

基于OpenStreetMap的地理信息绘图

开放街道地图（OpenStreetMap，简称OSM）是一个网上地图协作计划，目标是创造一个内容自由且能让所有人编辑的世界地图。OSM可以根据用户的手持GPS装置、航空摄影照片、卫星影像、其他自由内容甚至单靠用户对目标区域的知识绘制。OSM网站的灵感来自维基百科等网站。这可从该网地图页的“编辑”按钮及其完整修订历史获知。经注册的用户可上载GPS路径及使用内置的编辑程式编辑数据。2010年海地大地震中，OpenStreetMap极大地提高了地面搜救小队的工作效率。这个地图以GeoEye等公司提供的最新卫星照片为基础，又加入了很多新资料，工作人员和志愿者可以随时用随身的电脑或手机在地图上即时标注救护站、帐篷和倒塌的大桥。2012年4月，继苹果和Foursquare相继放弃使用谷歌地图后，维基百科也放弃使用谷歌地图，转向使用OpenStreetMap。

本例结合OpenStreetMap包、XLM包和ggplot2包来介绍如何整合地图资料和其它数据。我们的任务是将中国城区人口最多的前十位城市标注在地图上。基本步骤如下：

从wikipedia上抓取城市人口资料；
利用Google API获得十个城市的经纬度；
根据上面两类数据渲染得到地图如下(点击可看大图)，圆的大小显示了城市居民数的多寡。

R代码如下：

library(XML)
library(OpenStreetMap)
library(ggplot2)
# 从wiki获取中国前十大城市人口数据
url <-'http://zh.wikipedia.org/wiki/%E4%B8%AD%E5%9B%BD%E5%9F%8E%E5%B8%82%E5%88%97%E8%A1%A8'
tables <- readHTMLTable(url)
raw <- tables[[2]]
data <- raw[1:10,c(2,4,5)]
names(data) <- c('cname','name','pop')
data$name <- as.character(data$name)
# 将人口数从字符型转为数值型
data$pop <- as.numeric(gsub(',','',data$pop))
 
# 建立函数，输入参数为城市名称，返回经纬坐标
# 先利用google API得到xml形式结果，再用XML包函数提取经纬度
geoencode <- function(name) {
  requestUrl<-paste(
   "http://maps.googleapis.com/maps/api/geocode/xml?address=",
   name,"&sensor=false", sep="")
  xmlResult<-xmlTreeParse(requestUrl,isURL=TRUE)
  root <- xmlRoot(xmlResult)
  latitude <-xmlValue(root[['result']][['geometry']][['location']][['lat']])
  longitude <-xmlValue(root[['result']][['geometry']][['location']][['lng']])
  return(list(latitude=latitude,longitude=longitude))
}
# 得到十个城市的经纬度，结果为矩阵形式，整合为数据框
tempdata <- sapply(data$name,geoencode)
data <- cbind(data,data.frame(t(tempdata)))
data$latitude <-as.numeric(unlist(data$latitude))
data$longitude <-as.numeric(unlist(data$longitude))
# 将经纬度转为open street map绘图所需的输入参数
data <- cbind(data,projectMercator(data$latitude,data$longitude))
# 获取从东经100到125，北纬21到41之间的地图资料
map <- openmap(c(41,100),c(21,125),zoom=5,type = "osm")
# 利用ggplot2绘图包增加城市散点，点的大小表示人口
autoplot(map)+geom_point(data=data,aes(x=x,y=y
  ,size=pop),shape=16,colour='gold',alpha = 0.6)+
  geom_point(data=data,aes(x=x,y=y
  ,size=pop),shape=1,colour='red4')+
  scale_size_continuous(range= c(8,16))+
  opts(axis.line=theme_blank(),axis.text.x=theme_blank(),
       axis.text.y=theme_blank(),axis.ticks=theme_blank(),
       axis.title.x=theme_blank(),
       axis.title.y=theme_blank()) +
  geom_text(data=data,colour='black',
            aes(x=x,y=y,label=cname),hjust=1.4,vjust=0)
# 保存图片
ggsave("plot1.png", dpi=600)