写到这里，三大秘诀和三大法宝亮过相了，基本操作介绍过了，扩展包也华丽登过场了，基本上，菜鸟们已经可以独立折腾了，这个系列的帖子也写不下去了。不过，第 05 篇提到过画地图的问题，就再多写一篇介绍一下空间数据吧。
开胃小菜: 这就是 R。这里没有做不到，只有想不到。（This is R. There is no if. Only how. by Simon Blomberg.)
空间数据，也就是地理信息系统 GIS 的数据，很多人都用昂贵的 ArcGIS 软件来处理。其实，免费的R 配上强大的扩展包，也能够处理很多 GIS 问题，有时甚至更灵活。老实说， dapeng 对 GIS 用得不多，这一篇也就不得不讲得浅显，见谅了。感兴趣的菜鸟可以看这里、这里以及这本书。
这里，我们一起在中国地图上，各省份区域用不同颜色来区分，颜色代表该省 PM10 的质量浓度。请先点击这里下载地图。这是个压缩包，下载后请解压缩到 _c:\R\data\chinamap_ 文件夹。这个文件夹下应该有三个文件，即 bou2_4p.dbf，bou2_4p.shp，bou2_4p.shx。
require(rgdal) # 用来处理 GIS 数据的扩展包。若是第一次使用，请自行安装。
## Loading required package: rgdal
## Warning: package 'rgdal' was built under R version 2.15.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 2.15.3
## rgdal: version: 0.8-5, (SVN revision 449) Geospatial Data Abstraction
## Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
## 1.9.2, released 2012/10/08 Path to GDAL shared files: d:/Program
## Files/R/R-2.15.2/library/rgdal/gdal GDAL does not use iconv for recoding
## strings. Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009,
## [PJ_VERSION: 470] Path to PROJ.4 shared files: d:/Program
## Files/R/R-2.15.2/library/rgdal/proj
require(plotrix) # 用来画图例的扩展包。自行安装
## Loading required package: plotrix
cm <- readOGR("C:\\R\\data\\chinamap", "bou2_4p") # 读入 GIS 数据
## OGR data source with driver: ESRI Shapefile
## Source: "C:\R\data\chinamap", layer: "bou2_4p"
## with 925 features and 7 fields
## Feature type: wkbPolygon with 2 dimensions
summary(cm)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 73.447 135.09
## y 6.319 53.56
## Is projected: NA
## proj4string : [NA]
## Data attributes:
## AREA PERIMETER BOU2_4M_ BOU2_4M_ID
## Min. : 0.00 Min. : 0.01 Min. : 2 Min. : 0
## 1st Qu.: 0.00 1st Qu.: 0.02 1st Qu.:233 1st Qu.: 588
## Median : 0.00 Median : 0.04 Median :464 Median :2235
## Mean : 1.04 Mean : 1.50 Mean :464 Mean :1856
## 3rd Qu.: 0.00 3rd Qu.: 0.09 3rd Qu.:695 3rd Qu.:3004
## Max. :175.59 Max. :129.93 Max. :926 Max. :3338
##
## ADCODE93 ADCODE99 NAME
## Min. : 0 Min. : 0 浙江省 :179
## 1st Qu.:330000 1st Qu.:330000 福建省 :168
## Median :350000 Median :350000 广东省 :154
## Mean :405762 Mean :405751 辽宁省 : 94
## 3rd Qu.:440000 3rd Qu.:440000 山东省 : 86
## Max. :810000 Max. :810000 (Other):243
## NA's : 1
levels(cm$NAME) # 省市名称。
## [1] "安徽省" "北京市" "福建省"
## [4] "甘肃省" "广东省" "广西壮族自治区"
## [7] "贵州省" "海南省" "河北省"
## [10] "河南省" "黑龙江省" "湖北省"
## [13] "湖南省" "吉林省" "江苏省"
## [16] "江西省" "辽宁省" "内蒙古自治区"
## [19] "宁夏回族自治区" "青海省" "山东省"
## [22] "山西省" "陕西省" "上海市"
## [25] "四川省" "台湾省" "天津市"
## [28] "西藏自治区" "香港特别行政区" "新疆维吾尔自治区"
## [31] "云南省" "浙江省" "重庆市"
pm <- c(100, 141, 80, 174, 99, 72, 104, 30, 175, 107, 121, 133, 135, 98, 120,
100, 135, 116, 132, 139, 149, 172, 136, 97, 118, NA, 133, 65, NA, 127, 86,
119, 147) # 按照上一条指令得到的省市名称的顺序，对应的省会城市 PM10 质量浓度。恕 dapeng 无法提供数据来源，就当瞎编的吧。
pm.col <- cm.colors(diff(c(30, 180)) + 1)[(pm - min(pm, na.rm = TRUE)) + 1] # 将 PM10 浓度用颜色代码表示。
cm$col <- cm$NAME # 添加一列颜色列，与省市名称相同。
levels(cm$col) <- pm.col # 更改颜色列的因子。
plot(cm, col = as.character(factor(cm$col))) # 作图
color.legend(120, 20, 123, 2, c(expression("180 " * mu * "g m"^"-3"), expression("120 " *
mu * "g m"^"-3"), expression(" 30 " * mu * "g m"^"-3")), rev(cm.colors(diff(c(30,
100)) + 1)), align = "rb", gradient = "y") # 添加图例。
axis(1)
axis(2)
box()
颜色越蓝，PM10 浓度越低。跟北京相比，上海虽然有死猪，但 PM10 浓度还是要低很多的。
练习 11.1 从 http://www.diva-gis.org/gdata 免费下载你感兴趣的 GIS 地图，画出来，并更改各区域的颜色。
( 连载中，待续 )