写在前面

RgoogleMaps是SenseNetworks公司的Markus Loecher于2009年5月开始编写的一个R包. 此包的用途是把已有坐标数据作为一个个点自动打到从Google Maps上抓取的卫星地图上，好处是省去了人工计算抓取图片的参数和手动抓取、手动调整、打点的麻烦. 换句话说，现在我们只要在R中准备好经纬度数据再调用几个简单的函数就能画出一张 ( 至少从表面上看 ) 比较漂亮的图了 ( 其实也没这么简单 ).

话说敝人有去GIS选修课上睡过几次，奈何所学未有涉及地理专业，叙述方面难免缺失，只要大家能看懂就好了. 

准备工作

下载安装或读入以下R包：

sp

RgoogleMaps

rgdal （若需生成PNG格式图片）

ReadImages（若需生成JPEG格式图片，还要安装 libjpeg）

接着，前往以下网址申请一个Google Maps的API key. 至于API是什么，我想就不用解释了，不知道的童鞋可以去Google上百度一下.

http://code.google.com/apis/maps/signup.html

图1: Google Maps API 申请页面

阅毕许可协议，勾选 I have read and agree with the terms and conditions 即可获得一个Google Maps的API key. 将此API key（一长串的数字字母符号混合物）写入一个文本文档，保存为：

X:Documents and Settings你的用户名My Documents/API.key.txt

然后，将格式各异的原始数据处理为一个你我会读取的格式. 这里以我国地震监测部门四川台网所测定的2008年全年32640次地震震源的经纬度数据为例. 经过删减和规范，我们得到了一个体积与pattern俱佳的CSV文件，设其为data.csv.  点此下载

此包工作流程大致描述如下：

获得经纬度数据
↓

通过计算确定抓取图片所需参数

↓

访问Google服务器抓取和生成图片

↓

根据经纬度数据在图片上打点并输出

开工

markerdata<-read.csv("data.csv")

# 读取数据赋给对象markerdata. markerdata的作用有二：辅助确定抓取图片的参数；作为所打点的坐标.

boundingbox <- qbbox(lon = markerdata[,"longitudes"],

lat = markerdata[,"latitudes"])

# 通过对给出的经纬度数据进行简单的数学运算，确定要抓取图片的边界及中心经纬度数值，赋给对象boundingbox. 感兴趣的同学可以在R中直接输入qbbox看下作者的算法.

EQMap <- GetMap.bbox(boundingbox$lonR,
boundingbox$latR, destfile = "2008_SiChuan_EQs.png",
maptype="satellite", zoom=6)

# 显示出抓取图片的URI之后，即开始从Google的服务器上下载图片. 受Google的卫星图片数据本身的限制，图片的最大也是默认大小为640x640. 要是有某种高级卫星的数据就 …哈哈 .. 值得一提的是，如何设置缩放级别即参数zoom的值是一个值得斟酌的问题，zoom值太小时是大全景，打点时点都聚在一起看不清，zoom值太大了图像拉的太近待会打点时又无法打全.

图2: 抓取的原始图片

# p.s. 我做到这步总是有一个 Warning：

# In readGDAL(destfile, silent = TRUE) : GeoTransform values not available.

# 不过倒是能够正常完成抓取，不知何故，在此求高人解答.

result <- PlotOnStaticMap(EQMap,

lon=markerdata$longitudes, lat=markerdata$latitudes,

pch=20, cex = .1, col='red', verbose=0)

# 在已经读好的图片上打点. 其余的参数大家可以根据文档和自己的需要进行挖掘，如生成灰度图等. 参数verbose如果不为0，就会在绘图后给出若干行文字信息. 文档里清一色地解释为”level of verbosity”, 不知其具体意思，同求高人解答.

图3: 最终结果图

附：作者给出的示例

Coso地区位于美国的南加州，是美国第二大的地热发电厂所在地. geomapdata包中的cosomap.RData文件包含了该地地热区的断层等地质情况信息. 为了方便大家练手，此数据亦可在我这里下载.

首先，用R读入cosomap.RData. 接下来:

bb<-qbbox(lon=cosomap$POINTS$lon-360,

lat=cosomap$POINTS$lat)

# 注意这里需要对不满足Google Maps规定格式的原始经度数据进行简单换算.

MyMap<-GetMap.bbox(bb$lonR, bb$latR,

destfile = "Coso.png", maptype="satellite",zoom=11)

# 抓取图片，缩放级别为11，类型为默认的卫星地图. 除此之外，还有mobile, terrain, hybrid等地图类型.

tmp<-PlotOnStaticMap(MyMap,

lon=cosomap$POINTS$lon-360,lat=cosomap$POINTS$lat,

pch=20,cex = .5,col='red', verbose=0)

# 在图片上打点.

图4: Coso地区的断层分布

参考文献及网站

[1] Plotting on Google Static Maps in R
[2] Overlays on Google map tiles in R
[3] R for beginners 中文版
[4] 最近一周中国及周边版图地震情况
[5] 中国国家地震科学数据共享中心