星期四, 六月 28, 2012

用ggmap包进行地震数据的可视化

最近又发现了一个比较好玩的包ggmap。从名字上可以猜测出来,它的作用就是将ggplot2和map相结合。这样R语言用户能方便的获取各种静态地图数据,并在其基础上使用强大的ggplot绘图工具。ggmap包整合了四种地图资源,分别是Google、OpenStreetMaps、Stamen和Cloudmade。为了演示ggmap的作用,本例是从地震信息网获取最近一周的地震数据,得到其经纬度,然后以散点形式绘制在google地图上,另外也显示地震发生的密度估计。这个思路本质上和之前的一篇博文是一致的,但用ggmap包实现起来更为简单。

# 加载扩展包
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')+
                       opts(legend.position = "none")
 
更好玩的作法就是根据地震发生的日期生成不同的静态图,然后用animaiton包将其整合为一个gif动画。

# 为了生成动画,先准备好一个绘图函数
plotfunc <- function(x) {
    df <- subset(data,date <= x)
    df$lan <- as.numeric(df$lan)
    df$lon <- as.numeric(df$lon)
    p <- ggmap(get_googlemap(center = 'china', zoom=4,maptype='terrain'),,extent='device')+
        geom_point(data=df,aes(x=lon,y=lan),colour = 'red',alpha=0.7)
}
# 获取地震的日期
time <- sort(unique(data$date))
# 生成并保存动画
saveMovie(for( i in time) print(plotfunc(i)))

ggmap包中还有其它一些非常有用的函数。例如geocode函数可以根据地名字符串来查询经纬度,gglocator类似于基本包中的locator,它根据鼠标的点选来返回其坐标值。另外一个是mapdist函数,可以返回两点之间的地图距离和行驶时间。结合这几个函数可以直接在R中绘制地图,选择你的出发地和目标地,然后获得两地之间的距离。当然你还可以配合GPS等其它数据,创造出其它有意思的图形。

42 条评论:

  1. 看到这个,顿时觉得用ggplot什么直接加载GIS数据画地图爆弱了...拜一下。

    回复删除
  2. raw[-1,c('V1','V3','V4')] 不明白

    回复删除
  3. 应该是data <- tables(-1,c('V1','V3','V4'))

    回复删除
    回复
    1. raw是前面从原始网页里头抓下来的第六个表格,也就是有数据的可用的表格。我发现最近好些匿名的留言都被系统自动放到垃圾评论里头去了。各位路过的见谅啊。

      删除
  4. get_googlemap这个方法是哪个包的?怎么说找不到这个方法?

    回复删除
    回复
    1. 需要安装加载ggmap包
      install.packages('ggmap')
      library(ggmap)

      删除
  5. 此评论已被作者删除。

    回复删除
  6. 为什么说:错误于nchar(els[[1]]) : 多字节字符串8有错

    回复删除
    回复
    1. 不知道你在哪一步出的问题,看一下那个els变量是怎么回事。

      删除
    2. Error: unexpected 'in' in "help(Error in",我也出现差不多的问题,请博主指点一下

      删除
  7. 超级感谢!!!不过
    data <- raw[-1,c('V1','V3','V4')]这一步读不上数据啊

    回复删除
  8. 还有一个问题哈~~ legend怎么加在第一个图上的??

    回复删除
    回复
    1. 我刚试了,数据应该可以读进来的,可能是你那边网络问题?legend是通过控制图例坐标到右下角加上的

      删除
  9. 安装xml包发生错误,请问如何解决?

    回复删除
  10. 在运行“tables <- readHTMLTable(webpage,stringsAsFactors = FALSE)”出错,
    Error in nchar(els[[1]]) : invalid multibyte string 8
    博主,怎么办呢?

    回复删除
    回复
    1. 这部分是要连网操作的,因为是下载读取网络上的内容嘛,是否你没有连网?

      删除
    2. 我连了网还是不行

      删除
    3. 同样的问题
      联网也是不行的

      删除
  11. 老师您好,不知道您用的都是什么版本的程序,我运行的时候总是提示有问题

    回复删除
    回复
    1. 用的2.15,现在的3.0是可能有问题,所以我一直没升级。

      删除
    2. 老师,我用的ggmap_2.3,animation2.2,XML3.96,按照老师上面的代码读入的时候首先显示的错误是'opts' is deprecated. Use 'theme' instead. (Deprecated; last used in version 0.9.1),
      然后警告信息:
      1: Removed 931 rows containing non-finite values (stat_density2d).
      2: Removed 931 rows containing missing values (geom_point).
      有931条信息被删了,图上只有四个点,不知道原因是什么,、
      还有就是我不太理解as.numeric的作用是什么,为什么数值的跟原来差那么大
      刚接触R,不是很懂,还望老师能够给予解答

      删除
    3. 我目前的环境是x86_64-pc-linux-gnu (64-bit),R2.15.3,ggmap2.2,XML_3.95。刚刚尝试是可以画出静态图的,动画的没试。as.numeric是为了将字符串转为数值型。

      删除
    4. 上来跟老师您汇报一声,我成功试出来了,最后还是用到了Ubuntu,在Windows底下不行,老师能解释一下其中的原因吗,上次问您说as.numeric什么作用,其实从字面上就能了解是什么意思,只是在window下用这个函数之后的结果是跟字符串显示的数字是不一样的,不知道是为什么

      删除
  12. 您好,最近被R的可视化震撼到了,网上找了个例子,但是一直实现不出来,不知道您有什么办法,先谢谢老师了。http://spatialanalysis.co.uk/2012/06/mapping-worlds-biggest-airlines/

    回复删除
    回复
    1. 它那个NASA地图没下载过来,不过你参考这个例子也是可以的http://www.stanford.edu/~cengel/cgi-bin/anthrospace/great-circles-on-a-recentered-worldmap-in-ggplot,关键是画出弧线,就是靠geosphere包实现的,然后转数据框格式时需要一个函数。
      fortify.SpatialLinesDataFrame <- function(model, data, ...) {
      ldply(model@lines, fortify)
      }

      删除
  13. 请问Geocode和Mapdist如何实现中文地址解析呢?

    回复删除
  14. 你好,很高心看到你的作品。我现在的问题是不想让出来的地图中有那个等高线。怎么办?谢谢,我试了半天也没解决。

    回复删除
  15. 我 解决了刚问的那个问题。谢谢你啊

    回复删除
  16. get_googlemap()中有language = "en-EN"参数,但改为language = "zh-CN"后还是英文的,可是输出信息中有:
    Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=china&zoom=4&size=%20640x640&maptype=terrain&language=%20zh-CN&sensor=false
    Google Maps API Terms of Service : http://developers.google.com/maps/terms
    Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=china&sensor=false
    Google Maps API Terms of Service : http://developers.google.com/maps/terms
    Warning message:
    Removed 4 rows containing missing values (geom_point).

    如果直接打开上面的URL,能看到中文的地图.我的问题是还要哪些设置才能取到中文地图?

    回复删除
    回复
    1. 我也没有成功过。或许要改系统中的语言环境设置?

      删除
    2. > map <- get_googlemap('taichung city hall', markers = df, path = df, scale = 2,key='AIzaSyDyGhYwN3VMjxPER-DVnAF3TVsdMvG8P6w')
      錯誤在download.file(url, destfile = paste0("ggmapFileDrawer/", destfile), :
      無法開啟目標檔案 'ggmapFileDrawer/.png' ,原因是 'No such file or directory'

      為什麼會這樣?!

      删除
  17. 参照Git的修正,主要是语言选项缺了个sep

    回复删除
  18. 老师您好,我在使用saveMovie()函数的时候,出现了个错误:Error in data.frame(ll.lat = ll[1], ll.lon = ll[2], ur.lat = ur[1], ur.lon = ur[2]) :
    arguments imply differing number of rows: 0, 1
    动态图片一直没有显示出来,这个错误怎么解决啊?

    回复删除
  19. 老师你好,谢谢您的例子,让我这个R语言初学者受益匪浅。

    这个地震数据的动态图是叠加在前一天数据上生成的, 如果我只是想生成每一天的数据图,然后使用这些每天独立的静态图合成一个动态图反映点数据的分布变化,应该怎样做呢?

    思考了很久没有想出来,在此求教了!得到解答 不胜感激~

    回复删除
    回复
    1. 嘻嘻 我知道了 df <- subset (all, date < x) 这个 subset 改一下就可以了, 真是太笨了

      删除
  20. g <- ggmap(get_googlemap(center = 'china', zoom=4,maptype='terrain'), extent='device')
    plotfunc <- function(x) {
    df <- subset(data,date <= x)
    df$lan <- as.numeric(df$lan)
    df$lon <- as.numeric(df$lon)
    p <- geom_point(data=df,aes(x=lon,y=lan),colour = 'red',alpha=0.7)
    }

    回复删除
  21. 您好,请问,这个URL现在不可用,因为不懂地震,去地震网上找上找了半天,不知道该用哪个。可以再发一个可用的链接吗?谢谢。

    回复删除