R 菜鸟入门篇 第11篇 空间数据

写到这里,三大秘诀和三大法宝亮过相了,基本操作介绍过了,扩展包也华丽登过场了,基本上,菜鸟们已经可以独立折腾了,这个系列的帖子也写不下去了。不过,第 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()

plot of chunk unnamed-chunk-1

颜色越蓝,PM10 浓度越低。跟北京相比,上海虽然有死猪,但 PM10 浓度还是要低很多的。

练习 11.1 从 http://www.diva-gis.org/gdata 免费下载你感兴趣的 GIS 地图,画出来,并更改各区域的颜色。

( 连载中,待续 )

R 菜鸟入门篇 第11篇 空间数据》上有26条评论

  1. monomiao

    大鹏,你太能写了!
    你老板让你新加的数据处理做完了?

    回复
    1. dapeng 文章作者

      没有……就是因为论文写不下去,我那糟糕的英文写得我恶心,不得不写点中文来中和一下……

      回复
  2. 卡卡

    很有意思,不过有些参数和函数的细节没搞明白,只能照葫芦画瓢了

    回复
  3. ellie

    我点击下载地图,提示以下信息:This item might not exist or is no longer availableThis item might have been deleted, expired, or you might not have permission to view it. Contact the owner of this item for more information.

    回复
  4. ellie

    这个包rgdal安装不上,说sp的版本是1.0.8 要求>=1.0.9,但是sp包已经安装上了,需要卸载了,重新装吗?

    回复
  5. ellie

    已经可以了,原因是sp需要升级。需要重新关闭rstudio,然后安装sp 1.0.9,这样就可以安装上了

    回复
    1. dapeng 文章作者

      地图的事儿之前也有人提出来过下载不了,可是我这里测试都是可以下载的,奇怪。

      回复
  6. PureCotton

    大鹏 我的Rstudio貌似无法读中文 数据里如果有中文的话读出来都是乱码

    回复
    1. dapeng 文章作者

      编码的问题。Rstudio打开文件后,菜单栏File下选择Reopen with Encoding,找到合适的编码方式打开即可。

      回复
      1. PureCotton

        大鹏,找到你说的选项了,但是不知为何是灰色的,无法选择阿

        回复
        1. dapeng 文章作者

          重装一下试试?或者用别的文本编辑器转换一下编码。

        2. PureCotton

          大鹏 这个GIS的dbf文件在excel下打开中文列也显示不出来是为啥呢?平时我的excel显示中文没有问题啊

  7. 南南北

    求问那个cm.colors 怎么定义不同
    cm.colors(diff(c(30, 180)) + 1)[(pm – min(pm, na.rm = TRUE)) + 1]

    我加了自己的数据 出来没有差别

    回复
    1. 大鹏 文章作者

      是不是浓度范围的原因?贴来完整的代码看看。

      回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注

To create code blocks or other preformatted text, indent by four spaces:

    This will be displayed in a monospaced font. The first four 
    spaces will be stripped off, but all other whitespace
    will be preserved.
    
    Markdown is turned off in code blocks:
     [This is not a link](http://example.com)

To create not a block, but an inline code span, use backticks:

Here is some inline `code`.

For more help see http://daringfireball.net/projects/markdown/syntax